Skip to content

Commit

Permalink
GTI driver: recognize STAC GeoParquet catalogs
Browse files Browse the repository at this point in the history
  • Loading branch information
@rouault
rouault committed Aug 31, 2024
1 parent 9b1ab02 commit 4db9129
Show file tree
Hide file tree
Showing 4 changed files with 165 additions and 11 deletions.
23 changes: 23 additions & 0 deletions autotest/gdrivers/gti.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2947,3 +2947,26 @@ def test_gti_read_multi_threaded_disabled_because_truncated_source(tmp_vsimem):
assert vrt_ds.GetMetadataItem("MULTI_THREADED_RASTERIO_LAST_USED", "__DEBUG__") == (
"1" if gdal.GetNumCPUs() >= 2 else "0"
)


###############################################################################


@pytest.mark.require_curl()
@pytest.mark.require_driver("Parquet")
def test_gti_stac_geoparquet():

url = (
"https://github.com/stac-utils/stac-geoparquet/raw/main/tests/data/naip.parquet"
)

conn = gdaltest.gdalurlopen(url, timeout=4)
if conn is None:
pytest.skip("cannot open URL")

ds = gdal.Open("GTI:/vsicurl/" + url)
assert ds.GetSpatialRef().GetAuthorityCode(None) == "26914"
assert ds.GetGeoTransform() == pytest.approx(
(408231.0, 1.0, 0.0, 3873862.0, 0.0, -1.0), rel=1e-5
)
assert ds.RasterCount == 4
13 changes: 13 additions & 0 deletions doc/source/drivers/raster/gti.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,19 @@ The GTI driver accepts different types of connection strings:

For example: ``tileindex.gti``

STAC GeoParquet support
-----------------------

.. versionadded:: 3.10

The driver can support `STAC GeoParquet catalogs <https://stac-utils.github.io/stac-geoparquet/latest/spec/stac-geoparquet-spec>`_,
provided GDAL is build with :ref:`vector.parquet` support.
It can make use of fields ``proj:epsg`` and ``proj:transform`` from the
`Projection Extension Specification <https://github.com/stac-extensions/projection/>`_,
to correctly infer the appropriate projection and resolution.

Example of a valid connection string: ``GTI:/vsicurl/https://github.com/stac-utils/stac-geoparquet/raw/main/tests/data/naip.parquet``

Tile index requirements
-----------------------

Expand Down
3 changes: 2 additions & 1 deletion frmts/vrt/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,8 @@ add_gdal_driver(
gdaltileindexdataset.cpp
STRONG_CXX_WFLAGS)
gdal_standard_includes(gdal_vrt)
target_include_directories(gdal_vrt PRIVATE ${GDAL_RASTER_FORMAT_SOURCE_DIR}/raw)
target_include_directories(gdal_vrt PRIVATE ${GDAL_RASTER_FORMAT_SOURCE_DIR}/raw
$<TARGET_PROPERTY:ogrsf_generic,SOURCE_DIR>)

set(GDAL_DATA_FILES
${CMAKE_CURRENT_SOURCE_DIR}/data/gdalvrt.xsd
Expand Down
137 changes: 127 additions & 10 deletions frmts/vrt/gdaltileindexdataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
#include "vrtdataset.h"
#include "vrt_priv.h"
#include "ogrsf_frmts.h"
#include "ogrwarpedlayer.h"
#include "gdal_proxy.h"
#include "gdal_thread_pool.h"
#include "gdal_utils.h"
Expand Down Expand Up @@ -215,6 +216,9 @@ class GDALTileIndexDataset final : public GDALPamDataset
//! Vector layer with the sources
OGRLayer *m_poLayer = nullptr;

//! When the SRS of m_poLayer is not the one we expose
std::unique_ptr<OGRLayer> m_poWarpedLayerKeeper{};

//! Geotransform matrix of the tile index
std::array<double, 6> m_adfGeoTransform{0, 0, 0, 0, 0, 0};

Expand Down Expand Up @@ -583,6 +587,9 @@ GDALTileIndexDataset::GDALTileIndexDataset()
static std::string GetAbsoluteFileName(const char *pszTileName,
const char *pszVRTName)
{
if (STARTS_WITH(pszTileName, "https://"))
return std::string("/vsicurl/").append(pszTileName);

if (CPLIsFilenameRelative(pszTileName) &&
!STARTS_WITH(pszTileName, "<VRTDataset") &&
!STARTS_WITH(pszVRTName, "<GDALTileIndexDataset"))
Expand Down Expand Up @@ -897,13 +904,26 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
return false;
}

const auto poLayerDefn = m_poLayer->GetLayerDefn();

// Is this a https://stac-utils.github.io/stac-geoparquet/latest/spec/stac-geoparquet-spec ?
const bool bIsStacGeoParquet =
poLayerDefn->GetFieldIndex("assets.image.href") >= 0;

const char *pszLocationFieldName = GetOption(MD_LOCATION_FIELD);
if (!pszLocationFieldName)
{
constexpr const char *DEFAULT_LOCATION_FIELD_NAME = "location";
pszLocationFieldName = DEFAULT_LOCATION_FIELD_NAME;
if (bIsStacGeoParquet)
{
pszLocationFieldName = "assets.image.href";
}
else
{
constexpr const char *DEFAULT_LOCATION_FIELD_NAME = "location";
pszLocationFieldName = DEFAULT_LOCATION_FIELD_NAME;
}
}
auto poLayerDefn = m_poLayer->GetLayerDefn();

m_nLocationFieldIndex = poLayerDefn->GetFieldIndex(pszLocationFieldName);
if (m_nLocationFieldIndex < 0)
{
Expand Down Expand Up @@ -979,8 +999,8 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
const char *pszMinY = GetOption(MD_MINY);
const char *pszMaxX = GetOption(MD_MAXX);
const char *pszMaxY = GetOption(MD_MAXY);
const int nCountMinMaxXY = (pszMinX ? 1 : 0) + (pszMinY ? 1 : 0) +
(pszMaxX ? 1 : 0) + (pszMaxY ? 1 : 0);
int nCountMinMaxXY = (pszMinX ? 1 : 0) + (pszMinY ? 1 : 0) +
(pszMaxX ? 1 : 0) + (pszMaxY ? 1 : 0);
if (nCountMinMaxXY != 0 && nCountMinMaxXY != 4)
{
CPLError(CE_Failure, CPLE_AppDefined,
Expand Down Expand Up @@ -1022,11 +1042,10 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
return false;
}
}
else
else if (const auto poSRS = m_poLayer->GetSpatialRef())
{
const auto poSRS = m_poLayer->GetSpatialRef();
// Ignore GPKG "Undefined geographic SRS" and "Undefined Cartesian SRS"
if (poSRS && !STARTS_WITH(poSRS->GetName(), "Undefined "))
if (!STARTS_WITH(poSRS->GetName(), "Undefined "))
m_oSRS = *poSRS;
}

Expand Down Expand Up @@ -1086,15 +1105,113 @@ bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
}
}

// Take into STAC GeoParquet proj:epsg and proj:transform fields
std::unique_ptr<OGRFeature> poFeature;
std::string osResX, osResY, osMinX, osMinY, osMaxX, osMaxY;
int iProjEPSG = -1;
int iProjTransform = -1;
if (bIsStacGeoParquet && !pszSRS && !pszResX && !pszResY && !pszMinX &&
!pszMinY && !pszMaxX && !pszMaxY &&
(iProjEPSG = poLayerDefn->GetFieldIndex("proj:epsg")) >= 0 &&
(iProjTransform = poLayerDefn->GetFieldIndex("proj:transform")) >= 0)
{
poFeature.reset(m_poLayer->GetNextFeature());
if (poFeature && poFeature->IsFieldSet(iProjEPSG) &&
poFeature->IsFieldSet(iProjTransform))
{
const int nEPSGCode = poFeature->GetFieldAsInteger(iProjEPSG);
OGRSpatialReference oSTACSRS;
oSTACSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
if (nEPSGCode > 0 &&
oSTACSRS.importFromEPSG(nEPSGCode) == OGRERR_NONE)
{
int nTransformCount = 0;
const auto padfFeatureTransform =
poFeature->GetFieldAsDoubleList(iProjTransform,
&nTransformCount);
OGREnvelope sEnvelope;
if (nTransformCount >= 6 && m_poLayer->GetSpatialRef() &&
m_poLayer->GetExtent(&sEnvelope, /* bForce = */ true) ==
OGRERR_NONE)
{
const double dfResX = padfFeatureTransform[0];
osResX = CPLSPrintf("%.17g", dfResX);
const double dfResY = std::fabs(padfFeatureTransform[4]);
osResY = CPLSPrintf("%.17g", dfResY);

auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
OGRCreateCoordinateTransformation(
m_poLayer->GetSpatialRef(), &oSTACSRS));
auto poInvCT = std::unique_ptr<OGRCoordinateTransformation>(
poCT ? poCT->GetInverse() : nullptr);
double dfOutMinX = 0;
double dfOutMinY = 0;
double dfOutMaxX = 0;
double dfOutMaxY = 0;
if (dfResX > 0 && dfResY > 0 && poCT && poInvCT &&
poCT->TransformBounds(sEnvelope.MinX, sEnvelope.MinY,
sEnvelope.MaxX, sEnvelope.MaxY,
&dfOutMinX, &dfOutMinY,
&dfOutMaxX, &dfOutMaxY, 21))
{
constexpr double EPSILON = 1e-3;
const bool bTileAlignedOnRes =
(fmod(std::fabs(padfFeatureTransform[3]), dfResX) <=
EPSILON * dfResX &&
fmod(std::fabs(padfFeatureTransform[5]), dfResY) <=
EPSILON * dfResY);

osMinX = CPLSPrintf(
"%.17g",
!bTileAlignedOnRes
? dfOutMinX
: std::floor(dfOutMinX / dfResX) * dfResX);
osMinY = CPLSPrintf(
"%.17g",
!bTileAlignedOnRes
? dfOutMinY
: std::floor(dfOutMinY / dfResY) * dfResY);
osMaxX = CPLSPrintf(
"%.17g",
!bTileAlignedOnRes
? dfOutMaxX
: std::ceil(dfOutMaxX / dfResX) * dfResX);
osMaxY = CPLSPrintf(
"%.17g",
!bTileAlignedOnRes
? dfOutMaxY
: std::ceil(dfOutMaxY / dfResY) * dfResY);

m_oSRS = oSTACSRS;
pszResX = osResX.c_str();
pszResY = osResY.c_str();
pszMinX = osMinX.c_str();
pszMinY = osMinY.c_str();
pszMaxX = osMaxX.c_str();
pszMaxY = osMaxY.c_str();
nCountMinMaxXY = 4;

m_poWarpedLayerKeeper =
std::make_unique<OGRWarpedLayer>(
m_poLayer, /* iGeomField = */ 0,
/* bTakeOwnership = */ false, poCT.release(),
poInvCT.release());
m_poLayer = m_poWarpedLayerKeeper.get();
}
}
}
}
}

bool bHasMaskBand = false;
if ((!pszBandCount && apoXMLNodeBands.empty()) ||
(!(pszResX && pszResY) && nCountXSizeYSizeGT == 0))
{
CPLDebug("VRT", "Inspecting one feature due to missing metadata items");
m_bScannedOneFeatureAtOpening = true;

auto poFeature =
std::unique_ptr<OGRFeature>(m_poLayer->GetNextFeature());
if (!poFeature)
poFeature.reset(m_poLayer->GetNextFeature());
if (!poFeature ||
!poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
{
Expand Down

0 comments on commit 4db9129

Please sign in to comment.