gti_example.md

raster from multiple sources, disparate grid specs

This FlatGeoBuf file contains two polygon footprints, one is a northern hemisphere sea-ice concentration GeoTIFF image, the other southern (the values are scaled from 0-100% and stored with a palette, but we're just using the raw integer value here)

https://github.com/mdsumner/gti/raw/refs/heads/main/inst/extdata/tmerc_seaice.fgb

The source GeoTIFFs are in NSIDC south Polar Stereographic, and north Polar Stereographic as is standard. We use GTI format to register these as disparate data sources within a latent raster frame for rendering at a later time. We could store all the daily sea ice files this way, with an index on field 'date', so we could coalesce to a particular date (or set of dates, for examples like Sentinel scenes).

The south file is 316x332, and the north is 304x448.

Each GeoTIFF is stored in the 'location' field of the FlatGeoBuf, using the '/vsicurl/' protocol for streaming data from a url for GDAL.

The FlatGeoBuf itself sets Transverse Mercator (on longitude_0=147), and has metadata layer items that specify the resolution of the latent raster (RESX/RESY = 25000), as per the GTI specification.

Open as a raster with xarray, this is a "vector" file but we're declaring it a raster using the "GTI:" prefix, and we've set the minimal layer details for it to identify as a raster layer. At 25km resolution we get a grid of 487x1263.

import xarray
ds = xarray.open_dataset("GTI:https://github.com/mdsumner/gti/raw/refs/heads/main/inst/extdata/tmerc_seaice.fgb", engine = "rasterio")

ds.rio.crs
CRS.from_wkt('PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",147],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

<xarray.Dataset> Size: 2MB
Dimensions:      (band: 1, x: 487, y: 1263)
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 4kB -6.183e+06 -6.158e+06 ... 5.942e+06 5.967e+06
  * y            (y) float64 10kB 1.594e+07 1.591e+07 ... -1.559e+07 -1.561e+07
    spatial_ref  int64 8B ...
Data variables:
    band_data    (band, y, x) float32 2MB ...

The creation script is at this repo: https://github.com/mdsumner/gti

This plot sets the very highest value (land,missing) and zero concentrations to NA:

The two sources used are these, south then north:

"/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/10_Oct/S_20241017_concentration_v3.0.tif"

"/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/north/daily/geotiff/2024/10_Oct/N_20241017_concentration_v3.0.tif"

I use docker:

https://github.com/mdsumner/gdal-builds/pkgs/container/gdal-builds

rocker-gdal-dev-python

It doesn't have to all be so bleeding edge, but it's easier to keep up that way. GTI is GDAL 3.9, so it's very new

ok gotcha thanks, I'll give that a shot. It might be that the newest rasterio is using an older GDAL and overriding my python GDAL 3.9 install? Not sure.

Maybe, I went through packages to make sure they use my GDAL, PROJ, and GEOS but it's all bit chaotic now, sometimes I have to force reinstall pyproj and haven't got back to clean up yet

I think Pangeo docker has very recent GDAL so that might be a better way

Thanks for the tip, in a fresh conda environment with only xarray and rioxarray from the conda-forge channel this works.

Note that the warper doesn't need GTI you can just give it multiple sources, that's how the API works (I'm not really across how rasterio uses it, yet)