diff --git a/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_code.geojson b/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_code.geojson new file mode 100644 index 000000000000..b8c8a7b5da4e --- /dev/null +++ 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/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet.geojson b/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_epsg.geojson similarity index 73% rename from autotest/gdrivers/data/gti/sentinel2_stac_geoparquet.geojson rename to autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_epsg.geojson index 40048d297f81..9edeb19e291e 100644 --- a/autotest/gdrivers/data/gti/sentinel2_stac_geoparquet.geojson +++ b/autotest/gdrivers/data/gti/sentinel2_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 new file mode 100644 index 000000000000..c3b6f4094814 --- /dev/null +++ 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 new file mode 100644 index 000000000000..4686c8de852d --- /dev/null +++ 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 index e046df171b95..3c57323bc8ec 100755 --- a/autotest/gdrivers/gti.py +++ b/autotest/gdrivers/gti.py @@ -2991,7 +2991,16 @@ 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" @@ -2999,7 +3008,7 @@ def test_gti_stac_geoparquet_sentinel2(): 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 index f169e2f82298..ec4bd45929f2 100644 --- 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 `_, 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 `_, to correctly infer the appropriate projection and resolution. diff --git a/frmts/gti/gdaltileindexdataset.cpp b/frmts/gti/gdaltileindexdataset.cpp index f8a5217f456f..3877324665f8 100644 --- 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 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};