GTI driver: recognize STAC GeoParquet catalogs

[rouault]

@cholmes / @TomAugspurger

A toot from @mdsumner (https://mastodon.social/@[email protected]/113053401068184789) made me aware of STAC GeoParquet, and I've quickly realized this was very close to be handled by the recent GDAL Tile Index driver (https://gdal.org/en/latest/drivers/raster/gti.html), which is a generalization of the STACIT one, able to understand any OGR vector layer.

With this PR:

$ gdalinfo GTI:/vsicurl/https://github.com/stac-utils/stac-geoparquet/raw/main/tests/data/naip.parquet
Driver: GTI/GDAL Raster Tile Index
Files: none associated
Size is 23688, 14787
Coordinate System is:
PROJCRS["NAD83 / UTM zone 14N",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["UTM zone 14N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-99,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["North America - between 102°W and 96°W - onshore and offshore. Canada - Manitoba; Nunavut; Saskatchewan. United States (USA) - Iowa; Kansas; Minnesota; Nebraska; North Dakota; Oklahoma; South Dakota; Texas."],
        BBOX[25.83,-102,84,-96]],
    ID["EPSG",26914]]
Data axis to CRS axis mapping: 1,2
Origin = (408231.000000000000000,3873862.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  408231.000, 3873862.000) (100d 0'20.49"W, 35d 0'11.59"N)
Lower Left  (  408231.000, 3859075.000) (100d 0'14.63"W, 34d52'11.63"N)
Upper Right (  431919.000, 3873862.000) ( 99d44'46.03"W, 35d 0'18.33"N)
Lower Right (  431919.000, 3859075.000) ( 99d44'41.68"W, 34d52'18.34"N)
Center      (  420075.000, 3866468.500) ( 99d52'30.71"W, 34d56'15.22"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Description = Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Description = Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Description = Blue
Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined
  Description = NIR (near-infrared)

[mdsumner]

Nice, I love the convergence happening here. Some of it's really obvious but seen from very different lenses, and it's hard to see those without diving into other communities.

[coveralls]

coverage: 69.358% (+0.005%) from 69.353%
when pulling 47dbb5d on rouault:gti_stac_geoparquet
into 526b28c on OSGeo:master.

[cholmes]

Love it! Wasn't even aware of the GDAL Tile Index driver - sounds awesome.

[TomAugspurger]

Very neat!

[MathewNWSH]

Quick question: will it be possible to process STAC GeoParquet items provided via the GTI driver using https://gdal.org/en/latest/drivers/raster/vrt.html#default-pixel-functions?

[rouault]

Quick question: will it be possible to process STAC GeoParquet items provided via the GTI driver using https://gdal.org/en/latest/drivers/raster/vrt.html#default-pixel-functions?

no, the GTI driver only uses the "raw" items

[rouault]

no, the GTI driver only uses the "raw" items

but, what you could actually do is to "wrap" each item in a VRT using a pixel function and reference thoses VRT in the GTI catalog. Cool trick: the content of a VRT file can be used as a valid GDAL connection name, so you don't actually need to create a physical VRT file.

Cf

$ gdalinfo '<VRTDataset rasterXSize="20" rasterYSize="20"><VRTRasterBand dataType="Byte" band="1"><SimpleSource><SourceFilename relativeToVRT="0">/home/even/gdal/gdal/build_cmake/autotest/gcore/data/byte.tif</SourceFilename></SimpleSource></VRTRasterBand></VRTDataset>'
Driver: VRT/Virtual Raster
Files: /home/even/gdal/gdal/build_cmake/autotest/gcore/data/byte.tif
Size is 20, 20
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,   20.0)
Upper Right (   20.0,    0.0)
Lower Right (   20.0,   20.0)
Center      (   10.0,   10.0)
Band 1 Block=20x20 Type=Byte, ColorInterp=Undefined