VeniVidiVici
/
gdal
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #11514 from rouault/fix_11512 

GTI STACGeoParquet: add support for proj:code, proj:wkt2, proj:projjson
This commit is contained in:
Even Rouault 2024-12-19 19:02:22 +01:00 committed by GitHub
parent fbc0ed6dc8 27ff763aea
commit 989e655e02
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
7 changed files with 268 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_code.geojson Normal file
View File

@ -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 ] ] ] } }
]
}

2
autotest/gdrivers/data/gti/sentinel2_stac_geoparquet.geojson → autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_epsg.geojson
View File

@ -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 ] ] ] } }
]
}

194
autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_projjson.geojson Normal file
View File

@ -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 ] ] ] } }
]
}

7
autotest/gdrivers/data/gti/sentinel2_stac_geoparquet_proj_wkt2.geojson Normal file
View File

@ -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 ] ] ] } }
]
}

13
autotest/gdrivers/gti.py
View File

@ -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"

2
doc/source/drivers/raster/gti.rst
View File

@ -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.

58
frmts/gti/gdaltileindexdataset.cpp
View File

@ -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};