Merge pull request #11514 from rouault/fix_11512

GTI STACGeoParquet: add support for proj:code, proj:wkt2, proj:projjson
OSGeo · Dec 19, 2024 · 989e655 · 989e655
2 parents fbc0ed6 + 27ff763
commit 989e655
Show file tree

Hide file tree

Showing 7 changed files with 268 additions and 15 deletions.
diff --git a/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_code.geojson b/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_code.geojson
@@ -0,0 +1,7 @@
+{
+"type": "FeatureCollection",
+"name": "sentinel2_stac_geoparquet",
+"features": [
+{ "type": "Feature", "properties": { "type": "Feature", "stac_version": "1.1.0", "assets.rededge3.href": "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/12/S/VD/2023/12/S2A_T12SVD_20231213T181818_L2A/B07.tif", "assets.rededge3.title": "Red Edge 3 - 20m", "assets.rededge3.type": "image/tiff; application=geotiff; profile=cloud-optimized", "assets.rededge3.roles": [ "data", "reflectance" ], "assets.rededge3.eo:bands": [ { "name": "B07", "common_name": "rededge", "center_wavelength": 0.783, "full_width_half_max": 0.028 } ], "assets.rededge3.gsd": 20.0, "assets.rededge3.proj:shape": [ 5490, 5490 ], "assets.rededge3.proj:transform": [ 20, 0, 399960, 0, -20, 3900000 ], "assets.rededge3.raster:bands": [ { "nodata": 0, "data_type": "uint16", "spatial_resolution": 20, "scale": 0.0001, "offset": -0.1 } ], "assets.rededge3.file:checksum": "1220f2623a57da9b09fe8253c2dc9a6377227041e86d192715c6d753638cb0977d99", "assets.rededge3.file:size": 53583888, "proj:code": "EPSG:32612" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -112.099469246340007, 35.238074627405702 ], [ -112.086441301193005, 34.249018117524301 ], [ -111.171637815685003, 34.254252630787903 ], [ -110.892875963752005, 35.136386493702297 ], [ -110.892735914686995, 35.243021163215097 ], [ -112.099469246340007, 35.238074627405702 ] ] ] } }
+]
+}
diff --git a/...ata/gti/sentinel2_stac_geoparquet.geojson → ...ntinel2_stac_geoparquet_proj_epsg.geojson b/...ata/gti/sentinel2_stac_geoparquet.geojson → ...ntinel2_stac_geoparquet_proj_epsg.geojson
@@ -2,6 +2,6 @@
 "type": "FeatureCollection",
 "name": "sentinel2_stac_geoparquet",
 "features": [
-{ "type": "Feature", "properties": { "type": "Feature", "stac_version": "1.1.0", "assets.rededge3.href": "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/12/S/VD/2023/12/S2A_T12SVD_20231213T181818_L2A/B07.tif", "assets.rededge3.title": "Red Edge 3 - 20m", "assets.rededge3.type": "image/tiff; application=geotiff; profile=cloud-optimized", "assets.rededge3.roles": [ "data", "reflectance" ], "assets.rededge3.eo:bands": [ { "name": "B07", "common_name": "rededge", "center_wavelength": 0.783, "full_width_half_max": 0.028 } ], "assets.rededge3.gsd": 20.0, "assets.rededge3.proj:shape": [ 5490, 5490 ], "assets.rededge3.proj:transform": [ 20, 0, 399960, 0, -20, 3900000 ], "assets.rededge3.raster:bands": [ { "nodata": 0, "data_type": "uint16", "spatial_resolution": 20, "scale": 0.0001, "offset": -0.1 } ], "assets.rededge3.file:checksum": "1220f2623a57da9b09fe8253c2dc9a6377227041e86d192715c6d753638cb0977d99", "assets.rededge3.file:size": 53583888, "proj:epsg": 32612, }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -112.099469246340007, 35.238074627405702 ], [ -112.086441301193005, 34.249018117524301 ], [ -111.171637815685003, 34.254252630787903 ], [ -110.892875963752005, 35.136386493702297 ], [ -110.892735914686995, 35.243021163215097 ], [ -112.099469246340007, 35.238074627405702 ] ] ] } }
+{ "type": "Feature", "properties": { "type": "Feature", "stac_version": "1.1.0", "assets.rededge3.href": "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/12/S/VD/2023/12/S2A_T12SVD_20231213T181818_L2A/B07.tif", "assets.rededge3.title": "Red Edge 3 - 20m", "assets.rededge3.type": "image/tiff; application=geotiff; profile=cloud-optimized", "assets.rededge3.roles": [ "data", "reflectance" ], "assets.rededge3.eo:bands": [ { "name": "B07", "common_name": "rededge", "center_wavelength": 0.783, "full_width_half_max": 0.028 } ], "assets.rededge3.gsd": 20.0, "assets.rededge3.proj:shape": [ 5490, 5490 ], "assets.rededge3.proj:transform": [ 20, 0, 399960, 0, -20, 3900000 ], "assets.rededge3.raster:bands": [ { "nodata": 0, "data_type": "uint16", "spatial_resolution": 20, "scale": 0.0001, "offset": -0.1 } ], "assets.rededge3.file:checksum": "1220f2623a57da9b09fe8253c2dc9a6377227041e86d192715c6d753638cb0977d99", "assets.rededge3.file:size": 53583888, "proj:epsg": 32612 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -112.099469246340007, 35.238074627405702 ], [ -112.086441301193005, 34.249018117524301 ], [ -111.171637815685003, 34.254252630787903 ], [ -110.892875963752005, 35.136386493702297 ], [ -110.892735914686995, 35.243021163215097 ], [ -112.099469246340007, 35.238074627405702 ] ] ] } }
 ]
 }
diff --git a/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_projjson.geojson b/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_projjson.geojson
@@ -0,0 +1,194 @@
+{
+"type": "FeatureCollection",
+"name": "sentinel2_stac_geoparquet",
+"features": [
+{ "type": "Feature", "properties": { "type": "Feature", "stac_version": "1.1.0", "assets.rededge3.href": "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/12/S/VD/2023/12/S2A_T12SVD_20231213T181818_L2A/B07.tif", "assets.rededge3.title": "Red Edge 3 - 20m", "assets.rededge3.type": "image/tiff; application=geotiff; profile=cloud-optimized", "assets.rededge3.roles": [ "data", "reflectance" ], "assets.rededge3.eo:bands": [ { "name": "B07", "common_name": "rededge", "center_wavelength": 0.783, "full_width_half_max": 0.028 } ], "assets.rededge3.gsd": 20.0, "assets.rededge3.proj:shape": [ 5490, 5490 ], "assets.rededge3.proj:transform": [ 20, 0, 399960, 0, -20, 3900000 ], "assets.rededge3.raster:bands": [ { "nodata": 0, "data_type": "uint16", "spatial_resolution": 20, "scale": 0.0001, "offset": -0.1 } ], "assets.rededge3.file:checksum": "1220f2623a57da9b09fe8253c2dc9a6377227041e86d192715c6d753638cb0977d99", "assets.rededge3.file:size": 53583888, "proj:projjson": {
+  "$schema": "https://proj.org/schemas/v0.7/projjson.schema.json",
+  "type": "ProjectedCRS",
+  "name": "WGS 84 / UTM zone 12N",
+  "base_crs": {
+    "type": "GeographicCRS",
+    "name": "WGS 84",
+    "datum_ensemble": {
+      "name": "World Geodetic System 1984 ensemble",
+      "members": [
+        {
+          "name": "World Geodetic System 1984 (Transit)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1166
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G730)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1152
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G873)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1153
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G1150)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1154
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G1674)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1155
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G1762)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1156
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G2139)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1309
+          }
+        },
+        {
+          "name": "World Geodetic System 1984 (G2296)",
+          "id": {
+            "authority": "EPSG",
+            "code": 1383
+          }
+        }
+      ],
+      "ellipsoid": {
+        "name": "WGS 84",
+        "semi_major_axis": 6378137,
+        "inverse_flattening": 298.257223563
+      },
+      "accuracy": "2.0",
+      "id": {
+        "authority": "EPSG",
+        "code": 6326
+      }
+    },
+    "coordinate_system": {
+      "subtype": "ellipsoidal",
+      "axis": [
+        {
+          "name": "Geodetic latitude",
+          "abbreviation": "Lat",
+          "direction": "north",
+          "unit": "degree"
+        },
+        {
+          "name": "Geodetic longitude",
+          "abbreviation": "Lon",
+          "direction": "east",
+          "unit": "degree"
+        }
+      ]
+    },
+    "id": {
+      "authority": "EPSG",
+      "code": 4326
+    }
+  },
+  "conversion": {
+    "name": "UTM zone 12N",
+    "method": {
+      "name": "Transverse Mercator",
+      "id": {
+        "authority": "EPSG",
+        "code": 9807
+      }
+    },
+    "parameters": [
+      {
+        "name": "Latitude of natural origin",
+        "value": 0,
+        "unit": "degree",
+        "id": {
+          "authority": "EPSG",
+          "code": 8801
+        }
+      },
+      {
+        "name": "Longitude of natural origin",
+        "value": -111,
+        "unit": "degree",
+        "id": {
+          "authority": "EPSG",
+          "code": 8802
+        }
+      },
+      {
+        "name": "Scale factor at natural origin",
+        "value": 0.9996,
+        "unit": "unity",
+        "id": {
+          "authority": "EPSG",
+          "code": 8805
+        }
+      },
+      {
+        "name": "False easting",
+        "value": 500000,
+        "unit": "metre",
+        "id": {
+          "authority": "EPSG",
+          "code": 8806
+        }
+      },
+      {
+        "name": "False northing",
+        "value": 0,
+        "unit": "metre",
+        "id": {
+          "authority": "EPSG",
+          "code": 8807
+        }
+      }
+    ]
+  },
+  "coordinate_system": {
+    "subtype": "Cartesian",
+    "axis": [
+      {
+        "name": "Easting",
+        "abbreviation": "E",
+        "direction": "east",
+        "unit": "metre"
+      },
+      {
+        "name": "Northing",
+        "abbreviation": "N",
+        "direction": "north",
+        "unit": "metre"
+      }
+    ]
+  },
+  "scope": "Navigation and medium accuracy spatial referencing.",
+  "area": "Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA).",
+  "bbox": {
+    "south_latitude": 0,
+    "west_longitude": -114,
+    "north_latitude": 84,
+    "east_longitude": -108
+  },
+  "id": {
+    "authority": "EPSG",
+    "code": 32612
+  }
+}
+ }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -112.099469246340007, 35.238074627405702 ], [ -112.086441301193005, 34.249018117524301 ], [ -111.171637815685003, 34.254252630787903 ], [ -110.892875963752005, 35.136386493702297 ], [ -110.892735914686995, 35.243021163215097 ], [ -112.099469246340007, 35.238074627405702 ] ] ] } }
+]
+}
diff --git a/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_wkt2.geojson b/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_wkt2.geojson
@@ -0,0 +1,7 @@
+{
+"type": "FeatureCollection",
+"name": "sentinel2_stac_geoparquet",
+"features": [
+{ "type": "Feature", "properties": { "type": "Feature", "stac_version": "1.1.0", "assets.rededge3.href": "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/12/S/VD/2023/12/S2A_T12SVD_20231213T181818_L2A/B07.tif", "assets.rededge3.title": "Red Edge 3 - 20m", "assets.rededge3.type": "image/tiff; application=geotiff; profile=cloud-optimized", "assets.rededge3.roles": [ "data", "reflectance" ], "assets.rededge3.eo:bands": [ { "name": "B07", "common_name": "rededge", "center_wavelength": 0.783, "full_width_half_max": 0.028 } ], "assets.rededge3.gsd": 20.0, "assets.rededge3.proj:shape": [ 5490, 5490 ], "assets.rededge3.proj:transform": [ 20, 0, 399960, 0, -20, 3900000 ], "assets.rededge3.raster:bands": [ { "nodata": 0, "data_type": "uint16", "spatial_resolution": 20, "scale": 0.0001, "offset": -0.1 } ], "assets.rededge3.file:checksum": "1220f2623a57da9b09fe8253c2dc9a6377227041e86d192715c6d753638cb0977d99", "assets.rededge3.file:size": 53583888, "proj:wkt2": "PROJCRS[\"WGS 84 / UTM zone 12N\",BASEGEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 12N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-111,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[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA).\"],BBOX[0,-114,84,-108]],ID[\"EPSG\",32612]]" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -112.099469246340007, 35.238074627405702 ], [ -112.086441301193005, 34.249018117524301 ], [ -111.171637815685003, 34.254252630787903 ], [ -110.892875963752005, 35.136386493702297 ], [ -110.892735914686995, 35.243021163215097 ], [ -112.099469246340007, 35.238074627405702 ] ] ] } }
+]
+}
diff --git a/autotest/gdrivers/gti.py b/autotest/gdrivers/gti.py
@@ -2991,15 +2991,24 @@ def test_gti_stac_geoparquet():
 
 @pytest.mark.require_curl()
 @pytest.mark.require_driver("GeoJSON")
-def test_gti_stac_geoparquet_sentinel2():
+@pytest.mark.parametrize(
+    "filename",
+    [
+        "sentinel2_stac_geoparquet_proj_code.geojson",
+        "sentinel2_stac_geoparquet_proj_epsg.geojson",
+        "sentinel2_stac_geoparquet_proj_wkt2.geojson",
+        "sentinel2_stac_geoparquet_proj_projjson.geojson",
+    ],
+)
+def test_gti_stac_geoparquet_sentinel2(filename):
 
     url = "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/12/S/VD/2023/12/S2A_T12SVD_20231213T181818_L2A/B07.tif"
 
     conn = gdaltest.gdalurlopen(url, timeout=4)
     if conn is None:
         pytest.skip("cannot open URL")
 
-    ds = gdal.Open("GTI:data/gti/sentinel2_stac_geoparquet.geojson")
+    ds = gdal.Open(f"GTI:data/gti/{filename}")
     assert ds.RasterXSize == 5556
     assert ds.RasterYSize == 5540
     assert ds.GetSpatialRef().GetAuthorityCode(None) == "32612"

diff --git a/doc/source/drivers/raster/gti.rst b/doc/source/drivers/raster/gti.rst
@@ -81,7 +81,7 @@ STAC GeoParquet support
 
 The driver can support `STAC GeoParquet catalogs <https://stac-utils.github.io/stac-geoparquet/latest/spec/stac-geoparquet-spec>`_,
 provided GDAL is built with :ref:`vector.parquet` support.
-It can make use of fields ``proj:epsg`` and ``proj:transform`` from the
+It can make use of fields (``proj:code``, ``proj:epsg``, ``proj:wkt2``, or ``proj:projson``) and ``proj:transform`` from the
 `Projection Extension Specification <https://github.com/stac-extensions/projection/>`_,
 to correctly infer the appropriate projection and resolution.
 

diff --git a/frmts/gti/gdaltileindexdataset.cpp b/frmts/gti/gdaltileindexdataset.cpp
@@ -1150,10 +1150,14 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
         }
     }
 
-    // Take into STAC GeoParquet proj:epsg and proj:transform fields
+    // Take into STAC GeoParquet proj:code / proj:epsg / proj:wkt2 / proj:projjson
+    // and proj:transform fields
     std::unique_ptr<OGRFeature> poFeature;
     std::string osResX, osResY, osMinX, osMinY, osMaxX, osMaxY;
+    int iProjCode = -1;
     int iProjEPSG = -1;
+    int iProjWKT2 = -1;
+    int iProjPROJSON = -1;
     int iProjTransform = -1;
 
     const bool bIsStacGeoParquet =
@@ -1168,29 +1172,61 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
             strlen("assets."),
             osLocationFieldName.size() - strlen("assets.") - strlen(".href"));
     }
+
+    const auto GetAssetFieldIndex =
+        [poLayerDefn, &osAssetName](const char *pszFieldName)
+    {
+        const int idx = poLayerDefn->GetFieldIndex(
+            CPLSPrintf("assets.%s.%s", osAssetName.c_str(), pszFieldName));
+        if (idx >= 0)
+            return idx;
+        return poLayerDefn->GetFieldIndex(pszFieldName);
+    };
+
     if (bIsStacGeoParquet && !pszSRS && !pszResX && !pszResY && !pszMinX &&
         !pszMinY && !pszMaxX && !pszMaxY &&
-        ((iProjEPSG = poLayerDefn->GetFieldIndex(
-              CPLSPrintf("assets.%s.proj:epsg", osAssetName.c_str()))) >= 0 ||
-         (iProjEPSG = poLayerDefn->GetFieldIndex("proj:epsg")) >= 0) &&
-        ((iProjTransform = poLayerDefn->GetFieldIndex(CPLSPrintf(
-              "assets.%s.proj:transform", osAssetName.c_str()))) >= 0 ||
-         (iProjTransform = poLayerDefn->GetFieldIndex("proj:transform")) >= 0))
+        ((iProjCode = GetAssetFieldIndex("proj:code")) >= 0 ||
+         (iProjEPSG = GetAssetFieldIndex("proj:epsg")) >= 0 ||
+         (iProjWKT2 = GetAssetFieldIndex("proj:wkt2")) >= 0 ||
+         (iProjPROJSON = GetAssetFieldIndex("proj:projjson")) >= 0) &&
+        ((iProjTransform = GetAssetFieldIndex("proj:transform")) >= 0))
     {
         poFeature.reset(m_poLayer->GetNextFeature());
         const auto poProjTransformField =
             poLayerDefn->GetFieldDefn(iProjTransform);
-        if (poFeature && poFeature->IsFieldSet(iProjEPSG) &&
+        if (poFeature &&
+            ((iProjCode >= 0 && poFeature->IsFieldSet(iProjCode)) ||
+             (iProjEPSG >= 0 && poFeature->IsFieldSet(iProjEPSG)) ||
+             (iProjWKT2 >= 0 && poFeature->IsFieldSet(iProjWKT2)) ||
+             (iProjPROJSON >= 0 && poFeature->IsFieldSet(iProjPROJSON))) &&
             poFeature->IsFieldSet(iProjTransform) &&
             (poProjTransformField->GetType() == OFTRealList ||
              poProjTransformField->GetType() == OFTIntegerList ||
              poProjTransformField->GetType() == OFTInteger64List))
         {
-            const int nEPSGCode = poFeature->GetFieldAsInteger(iProjEPSG);
             OGRSpatialReference oSTACSRS;
             oSTACSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
-            if (nEPSGCode > 0 &&
-                oSTACSRS.importFromEPSG(nEPSGCode) == OGRERR_NONE)
+
+            if (iProjCode >= 0 && poFeature->IsFieldSet(iProjCode))
+                oSTACSRS.SetFromUserInput(
+                    poFeature->GetFieldAsString(iProjCode),
+                    OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
+
+            else if (iProjEPSG >= 0 && poFeature->IsFieldSet(iProjEPSG))
+                oSTACSRS.importFromEPSG(
+                    poFeature->GetFieldAsInteger(iProjEPSG));
+
+            else if (iProjWKT2 >= 0 && poFeature->IsFieldSet(iProjWKT2))
+                oSTACSRS.SetFromUserInput(
+                    poFeature->GetFieldAsString(iProjWKT2),
+                    OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
+
+            else if (iProjPROJSON >= 0 && poFeature->IsFieldSet(iProjPROJSON))
+                oSTACSRS.SetFromUserInput(
+                    poFeature->GetFieldAsString(iProjPROJSON),
+                    OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get());
+
+            if (!oSTACSRS.IsEmpty())
             {
                 int nTransformCount = 0;
                 double adfGeoTransform[6] = {0, 0, 0, 0, 0, 0};