Portfolio for Scott Parks: VRTMate

GDAL virtual raster (.vrt) is a popular way of managing a collection of rasters as a single vitrual raster. VRT's contain no data, only the path to the data, spatial information and processing information. VRT's can be unlimited in size. Our demo .vrt will be the USGS 10 meter DEM, USGS_Seamless_DEM_13.vrt, that contains 1,430 Cloud Optimized Geotiff's (COG). All remote, no local data. For displaying individual rasters currently in the viewport (bbox) on an interactive map, we want the path to the data and spatial information for all rasters in the .vrt.

To see a full working interactive web map of how VRT, COG and FGB cloud-native data work together, view the Serverless section in the notebook.

The JavaScript here is actually running. The code and results are displayed below. You can also open your browser inspector (Ctrl+Chift+I) to see output.

Documentation for VRTMate:

VRTMate does not exhaustively process a .vrt the way GDAL does. We are only retreiving file paths, overall and individual raster spatial information contained in the .vrt.

Class creation options:
   VRTMate(text vrt)
      vrt: URL to remote .vrt, use http(s)://
   
Methods:
   readVRT()
      Loads meta data for each raster. Returns .vrt overall meta data on success.
      
   getSize()
      returns the size of VRTMate in memory
      
   rastersByBBox(Object { bbox: {xmin: float, ymin: float, xmax: float, ymax: float}, bands: int[] })
      Returns meta data of all rasters in .vrt intersecting bounding box
      Object { 
         bbox: { xmin: -125, xmax: -120, ymin: 45, ymax: 49 }
         ,bands: [3,11] // Optional, default 1. Array of bands to from multi-band .vrt
      }
   
   search(object obj, text searchKey)
      obj: Any JS object, but particularly for a single raster from array of rasters returned by rastersByBBox or property 'rasters'
      ,searchKey: object key to find in obj tree
      
Properties (after call to readVRT()):
   info   : object, overall .vrt meta data.
   rasters: array of objects. Meta data about each raster in .vrt.

A note about raster urls/file paths in a .vrt. For web client use we need a path to all rasters that begin with http(s)://. In the .vrt a raster url/path may be a path that begins /vsicurl/... or /vsis3/... or maybe a local file system path. It's important to have a 'prefix' defined based on whatever is in the .vrt. In addition to a prefix, you may need to get just part of the raster path. Maybe the basename or the directory above the basename. You'll see 'basenameSlice' in the config. This is passed to array.slice, to parse our final raster path. Use any method you wish.

Viewing the data as a single raster from the command line:

gdal raster info --of=text https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt

# ...or

gdalinfo https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt

We'll be using our VRTMate JavaScript class to parse and build meta data for each raster in the .vrt. This config is not used directly by the class, for organization only. We'll only use the 'url' attribute to create the new instance:

let cfg = {
   url: "https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt"
   ,prefix: 'https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/'
   ,basenameSlice: -2 
}

Now, parse our .vrt and get a new instance to work with:

var vrt = new VRTMate(cfg.url);
let info = await vrt.readVRT();

Now, display spatial data that encompasses all rasters in the .vrt (example, not actual):

console.log('VRT meta: ', info); 

{
  "vrt": "https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt"
  "geoTransform": {
    "xmin": -180.0005555560936,
    "xres": 9.259259266022042e-05,
    "xoff": 0,
    "ymax": 72.0005555562949,
    "yoff": 0,
    "yres": -9.25925926599094e-05
  },
  "SRS": "\"EPSG\",\"4269\"",
  "EPSG": 4269,
  "bbox": {
    "xmin": -180.0005555560936,
    "xmax": 180.00000026239937,
    "ymin": -15.000555618067892,
    "ymax": 72.0005555562949
  },
  "width": 3888006,
  "height": 939612,
  "rasterCnt": 1430,
  "center": [
    28.499999969113503,
    -0.00027764684710973597
  ]
}

Using vrt.getSize() we can get the size of the vrt in memory:

console.log('VRT in memory size:', vrt.getSize());

1412220

Using vrt.search() we can get the overall .vrt bounding box:

console.log('BBOX vrt.search: ', vrt.search(vrt.info, 'bbox'));

BBOX vrt.search:  
Array [ {…} ]
   0: Object { key: "bbox", value: {…} }
      key: "bbox"
      value: Object { xmin: -180.0005555560936, xmax: 180.00000026239937, ymin: -15.000555618067892, … }
         xmax: 180.00000026239937
         xmin: -180.0005555560936
         ymax: 72.0005555562949
         ymin: -15.000555618067892

Supply bounding box and optionally raster bands. In a web client you would dynamically create the bounding box based on the map viewport. A combined zoom and bounding box strategy ultimately should decide the number of rasters to return. The more rasters and larger the bbox, the greater number of rasters returned.

Iterate through the results display spatial/path result for each individual raster in the .vrt. It's important to notice how configuration option 'prefix' and 'basenameSlice' work together here, compared to the 'vsicurl' in the JSON. Note the actual url after the JSON. When adding prefix/basenameSlice to your config, you'll have to examine the source .vrt to figure this out:

let rasters = vrt.rastersByBBox({ xmin: -125, xmax: -120, ymin: 45, ymax: 49 }, 3); 

for(key in rasters) {
   let url = 'basenameSlice' in cfg ? cfg.prefix + rasters[key].raster.split('/').slice(cfg.basenameSlice).join('/') : rasters[key].raster;
   console.log(rasters[key]);
   console.log(url);
}

Rasters meeting bbox criteria (example, not actual):

{
  "vrtSeq": 793,
  "bbox": {
    "xmin": -120.00055551227075,
    "xmax": -118.99944440042844,
    "ymax": 45.00055553666532,
    "ymin": 43.99944442482638,
    "center": {
      "y": 44.49999998074585,
      "x": -119.4999999563496
    }
  },
  "geoWidth": 1.001111111842306,
  "geoHeight": 1.001111111838938,
  "pixelWidth": 10812,
  "pixelHeight": 10811,
  "sort": 43120,
  "raster": "/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n45w120/USGS_13_n45w120.tif"
}

"https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n45w120/USGS_13_n45w120.tif" 

{
  "vrtSeq": 794,
  "bbox": {
    "xmin": -121.00055551300113,
    "xmax": -119.99944440115883,
    "ymax": 45.00055553666532,
    "ymin": 43.99944442482638,
    "center": {
      "y": 44.49999998074585,
      "x": -120.49999995707998
    }
  },
  "geoWidth": 1.001111111842306,
  "geoHeight": 1.001111111838938,
  "pixelWidth": 10812,
  "pixelHeight": 10811,
  "sort": 43121,
  "raster": "/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n45w121/USGS_13_n45w121.tif"
}

"https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n45w121/USGS_13_n45w121.tif"

{
  "vrtSeq": 795,
  "bbox": {
    "xmin": -122.00055551373151,
    "xmax": -120.9994444018892,
    "ymax": 45.00055553666532,
    "ymin": 43.99944442482638,
    "center": {
      "y": 44.49999998074585,
      "x": -121.49999995781036
    }
  },
  "geoWidth": 1.001111111842306,
  "geoHeight": 1.001111111838938,
  "pixelWidth": 10812,
  "pixelHeight": 10811,
  "sort": 43122,
  "raster": "/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n45w122/USGS_13_n45w122.tif"
}

"https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n45w122/USGS_13_n45w122.tif"

That's it!