Visualizing unlimited Cloud Optimized GeoTiffs(COG) with single FlatGeobuf(FGB)
Partial set of CONUS 30m DEM's, 514 Cloud Optimized GeoTiffs, sourced from FGB file, zoom level 8+
...or just STAC killer. That's what I've been calling it.
So, here it is. Let's create an FGB file with unlimited features. Each feature has 2 attributes:
We do this with our old friend, gdaltindex. We pass a whole slew of COG's and it returns them as an FGB:
gdaltindex -lyr_name dem30 -tileindex file dem30m.fgb /path/to/cogs/*.tif
This leads to a very interesting advantage. FGB is a cloud native or streaming format. Using the built in indexing we can get byte ranges of data out of the file, without reading the entire file, whether it's on a web server or cheap cloud storage like S3. This doesn't require any intermediate servers or services. No API's, just a web browser and the data. Using our highly efficient streaming format, we can get a bounding box of just the COG's we need, from an otherwise unlimited number of COG's!
See where I'm going with this?
By specifying a bounding box, we can get a list of COG's that intersect that bounding box. In the example, that bounding box is a web browser view port that's running LeafletJS. I'm also using Georaster-for-leaflet and flatgeobuf-geojson with Leaflet. This also works with Open Layers
Time to throw a bit of code at you. The following is the core piece of JavaScript that runs in your browser each time you pan/move/zoom your map:
// get browser viewport bbox, spatial extent
let b = o.map.getBounds();
// read only the features intersecting the bbox from the FGB file
for await (let feature of flatgeobuf.deserialize(o.opts.fgbUrl, {minX: b.getWest() ,maxX: b.getEast() ,minY: b.getSouth() ,maxY: b.getNorth()})) {
// skip COG if already loaded
if(feature.properties.file in o.loadedCogs) continue;
// add COG layer to map
parseGeoraster(o.opts.cogUrl + feature.properties.file, {readOnDemand: true}).then(georaster => {
o.georasterOpts.georaster = georaster;
o.loadedCogs[feature.properties.file] = new GeoRasterLayer(o.georasterOpts);
this.loadedCogs[feature.properties.file].addTo(o.map);
});
}
That's the basic premise. But so much more can be done. Here's a feature from our FGB file:
OGRFeature(dem30):1 file (String) = n30w082.tif POLYGON ((-82.00166666667 30.0016666666641,-80.9983333333359 30.0016666666641,-80.9983333333359 28.99833333333,-82.00166666667 28.99833333333,-82.00166666667 30.0016666666641))
The feature has only 2 attributes. That's all we need to make it work. So much more information can be added, such as meta data. Perhaps we added some pre-calculated data as well. A feature may look something like this:
OGRFeature(dem30):1 file (String) = n30w082.tif group (String) = CONUS desc (String) = 30 meter raster minElev (real) = 123.456 maxElev (real) = 654.321 avgElev (real) = 777.777 POLYGON ((-82.00166666667 30.0016666666641,-80.9983333333359 30.0016666666641,-80.9983333333359 28.99833333333,-82.00166666667 28.99833333333,-82.00166666667 30.0016666666641))
Once we add more attributes to the feature, we are able to filter in a whole new way. Initially, we want to get our COG's within a specific area. With the bounding box we've done that. In the browser using JavaScript, we can filter by any attribute. On the server side, being that this is just OGR, we can use SQL queries in our GDAL ogr* commands. Further, you could add JSON with any free from type of data as an attribute.
Further, we have put only COG's in the FGB. There's no reason we couldn't put other file types that have a known spatial extent. JSON, CSV, GPKG, whatever. Don't expect those formats to perform like COG. As long as your browser JavaScript has a handler for the format, it will work.
If you like the structure that STAC brings, but don't want to interact with an API, the strict structure or the complexity, this might be an alternative.
Apr 10 2024 - Scott Parks