Skip to content

Commit

Permalink
Merge pull request OSGeo#10459 from rouault/gti_GetDataCoverageStatus
Browse files Browse the repository at this point in the history 
GTI: implement GetDataCoverageStatus()
  • Loading branch information
@rouault
rouault authored Jul 22, 2024
2 parents 6e095f3 + 172f0b0 commit 15866aa
Show file tree
Hide file tree
Showing 2 changed files with 210 additions and 0 deletions.
62 changes: 62 additions & 0 deletions autotest/gdrivers/gti.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
import struct

import gdaltest
import ogrtest
import pytest

from osgeo import gdal, ogr
Expand Down Expand Up @@ -1057,6 +1058,20 @@ def test_gti_rgb_left_right(tmp_vsimem):
== "<LocationInfo><File>/vsimem/test_gti_rgb_left_right/left.tif</File></LocationInfo>"
)

if ogrtest.have_geos():
(flags, pct) = vrt_ds.GetRasterBand(1).GetDataCoverageStatus(
0, 0, vrt_ds.RasterXSize, vrt_ds.RasterYSize
)
assert flags == gdal.GDAL_DATA_COVERAGE_STATUS_DATA and pct == 100.0

(flags, pct) = vrt_ds.GetRasterBand(1).GetDataCoverageStatus(1, 2, 3, 4)
assert flags == gdal.GDAL_DATA_COVERAGE_STATUS_DATA and pct == 100.0

(flags, pct) = vrt_ds.GetRasterBand(1).GetDataCoverageStatus(
vrt_ds.RasterXSize // 2 - 1, 2, 2, 4
)
assert flags == gdal.GDAL_DATA_COVERAGE_STATUS_DATA and pct == 100.0


def test_gti_overlapping_sources(tmp_vsimem):

Expand All @@ -1082,6 +1097,12 @@ def test_gti_overlapping_sources(tmp_vsimem):
vrt_ds = gdal.Open(index_filename)
assert vrt_ds.GetRasterBand(1).Checksum() == 2

if ogrtest.have_geos():
(flags, pct) = vrt_ds.GetRasterBand(1).GetDataCoverageStatus(
0, 0, vrt_ds.RasterXSize, vrt_ds.RasterYSize
)
assert flags == gdal.GDAL_DATA_COVERAGE_STATUS_DATA and pct == 100.0

# Test unsupported sort_field_type = OFTBinary
index_filename = str(tmp_vsimem / "index.gti.gpkg")
sort_values = [None, None]
Expand Down Expand Up @@ -1319,6 +1340,41 @@ def test_gti_overlapping_sources(tmp_vsimem):
assert vrt_ds.GetRasterBand(1).Checksum() == 2, sort_values


def test_gti_gap_between_sources(tmp_vsimem):

filename1 = str(tmp_vsimem / "one.tif")
ds = gdal.GetDriverByName("GTiff").Create(filename1, 1, 1)
ds.SetGeoTransform([2, 1, 0, 49, 0, -1])
ds.GetRasterBand(1).Fill(1)
del ds

filename2 = str(tmp_vsimem / "two.tif")
ds = gdal.GetDriverByName("GTiff").Create(filename2, 1, 1)
ds.SetGeoTransform([4, 1, 0, 49, 0, -1])
ds.GetRasterBand(1).Fill(2)
del ds

index_filename = str(tmp_vsimem / "index.gti.gpkg")
index_ds, _ = create_basic_tileindex(
index_filename, [gdal.Open(filename1), gdal.Open(filename2)]
)
del index_ds

vrt_ds = gdal.Open(index_filename)
assert vrt_ds.GetRasterBand(1).Checksum() == 3

if ogrtest.have_geos():
(flags, pct) = vrt_ds.GetRasterBand(1).GetDataCoverageStatus(
0, 0, vrt_ds.RasterXSize, vrt_ds.RasterYSize
)
assert (
flags
== gdal.GDAL_DATA_COVERAGE_STATUS_DATA
| gdal.GDAL_DATA_COVERAGE_STATUS_EMPTY
and pct == pytest.approx(100.0 * 2 / 3)
)


def test_gti_no_source(tmp_vsimem):

index_filename = str(tmp_vsimem / "index.gti.gpkg")
Expand Down Expand Up @@ -1359,6 +1415,12 @@ def test_gti_no_source(tmp_vsimem):
is None
)

if ogrtest.have_geos():
(flags, pct) = vrt_ds.GetRasterBand(1).GetDataCoverageStatus(
0, 0, vrt_ds.RasterXSize, vrt_ds.RasterYSize
)
assert flags == gdal.GDAL_DATA_COVERAGE_STATUS_EMPTY and pct == 0.0


def test_gti_invalid_source(tmp_vsimem):

Expand Down
148 changes: 148 additions & 0 deletions frmts/vrt/gdaltileindexdataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,6 +370,9 @@ class GDALTileIndexBand final : public GDALPamRasterBand
GSpacing nLineSpace,
GDALRasterIOExtraArg *psExtraArg) override;

int IGetDataCoverageStatus(int nXOff, int nYOff, int nXSize, int nYSize,
int nMaskFlagStop, double *pdfDataPct) override;

int GetMaskFlags() override
{
if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)
Expand Down Expand Up @@ -2328,6 +2331,151 @@ CPLErr GDALTileIndexBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
nPixelSpace, nLineSpace, 0, psExtraArg);
}

/************************************************************************/
/* IGetDataCoverageStatus() */
/************************************************************************/

#ifndef HAVE_GEOS
int GDALTileIndexBand::IGetDataCoverageStatus(int /* nXOff */, int /* nYOff */,
int /* nXSize */,
int /* nYSize */,
int /* nMaskFlagStop */,
double *pdfDataPct)
{
if (pdfDataPct != nullptr)
*pdfDataPct = -1.0;
return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
GDAL_DATA_COVERAGE_STATUS_DATA;
}
#else
int GDALTileIndexBand::IGetDataCoverageStatus(int nXOff, int nYOff, int nXSize,
int nYSize, int nMaskFlagStop,
double *pdfDataPct)
{
if (pdfDataPct != nullptr)
*pdfDataPct = -1.0;

const double dfMinX = m_poDS->m_adfGeoTransform[GT_TOPLEFT_X] +
nXOff * m_poDS->m_adfGeoTransform[GT_WE_RES];
const double dfMaxX =
dfMinX + nXSize * m_poDS->m_adfGeoTransform[GT_WE_RES];
const double dfMaxY = m_poDS->m_adfGeoTransform[GT_TOPLEFT_Y] +
nYOff * m_poDS->m_adfGeoTransform[GT_NS_RES];
const double dfMinY =
dfMaxY + nYSize * m_poDS->m_adfGeoTransform[GT_NS_RES];
m_poDS->m_poLayer->SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
m_poDS->m_poLayer->ResetReading();

int nStatus = 0;

auto poPolyNonCoveredBySources = std::make_unique<OGRPolygon>();
{
auto poLR = std::make_unique<OGRLinearRing>();
poLR->addPoint(nXOff, nYOff);
poLR->addPoint(nXOff, nYOff + nYSize);
poLR->addPoint(nXOff + nXSize, nYOff + nYSize);
poLR->addPoint(nXOff + nXSize, nYOff);
poLR->addPoint(nXOff, nYOff);
poPolyNonCoveredBySources->addRingDirectly(poLR.release());
}
while (true)
{
auto poFeature =
std::unique_ptr<OGRFeature>(m_poDS->m_poLayer->GetNextFeature());
if (!poFeature)
break;
if (!poFeature->IsFieldSetAndNotNull(m_poDS->m_nLocationFieldIndex))
{
continue;
}

const auto poGeom = poFeature->GetGeometryRef();
if (!poGeom || poGeom->IsEmpty())
continue;

OGREnvelope sSourceEnvelope;
poGeom->getEnvelope(&sSourceEnvelope);

const double dfDstXOff =
std::max<double>(nXOff, (sSourceEnvelope.MinX -
m_poDS->m_adfGeoTransform[GT_TOPLEFT_X]) /
m_poDS->m_adfGeoTransform[GT_WE_RES]);
const double dfDstXOff2 = std::min<double>(
nXOff + nXSize,
(sSourceEnvelope.MaxX - m_poDS->m_adfGeoTransform[GT_TOPLEFT_X]) /
m_poDS->m_adfGeoTransform[GT_WE_RES]);
const double dfDstYOff =
std::max<double>(nYOff, (sSourceEnvelope.MaxY -
m_poDS->m_adfGeoTransform[GT_TOPLEFT_Y]) /
m_poDS->m_adfGeoTransform[GT_NS_RES]);
const double dfDstYOff2 = std::min<double>(
nYOff + nYSize,
(sSourceEnvelope.MinY - m_poDS->m_adfGeoTransform[GT_TOPLEFT_Y]) /
m_poDS->m_adfGeoTransform[GT_NS_RES]);

// CPLDebug("GTI", "dfDstXOff=%f, dfDstXOff2=%f, dfDstYOff=%f, dfDstYOff2=%f",
// dfDstXOff, dfDstXOff2, dfDstYOff, dfDstXOff2);

// Check if the AOI is fully inside the source
if (nXOff >= dfDstXOff && nYOff >= dfDstYOff &&
nXOff + nXSize <= dfDstXOff2 && nYOff + nYSize <= dfDstYOff2)
{
if (pdfDataPct)
*pdfDataPct = 100.0;
return GDAL_DATA_COVERAGE_STATUS_DATA;
}

// Check intersection of bounding boxes.
if (dfDstXOff2 > nXOff && dfDstYOff2 > nYOff &&
dfDstXOff < nXOff + nXSize && dfDstYOff < nYOff + nYSize)
{
nStatus |= GDAL_DATA_COVERAGE_STATUS_DATA;
if (poPolyNonCoveredBySources)
{
OGRPolygon oPolySource;
auto poLR = std::make_unique<OGRLinearRing>();
poLR->addPoint(dfDstXOff, dfDstYOff);
poLR->addPoint(dfDstXOff, dfDstYOff2);
poLR->addPoint(dfDstXOff2, dfDstYOff2);
poLR->addPoint(dfDstXOff2, dfDstYOff);
poLR->addPoint(dfDstXOff, dfDstYOff);
oPolySource.addRingDirectly(poLR.release());
auto poRes = std::unique_ptr<OGRGeometry>(
poPolyNonCoveredBySources->Difference(&oPolySource));
if (poRes && poRes->IsEmpty())
{
if (pdfDataPct)
*pdfDataPct = 100.0;
return GDAL_DATA_COVERAGE_STATUS_DATA;
}
else if (poRes && poRes->getGeometryType() == wkbPolygon)
{
poPolyNonCoveredBySources.reset(
poRes.release()->toPolygon());
}
else
{
poPolyNonCoveredBySources.reset();
}
}
}
if (nMaskFlagStop != 0 && (nStatus & nMaskFlagStop) != 0)
{
return nStatus;
}
}
if (poPolyNonCoveredBySources)
{
if (!poPolyNonCoveredBySources->IsEmpty())
nStatus |= GDAL_DATA_COVERAGE_STATUS_EMPTY;
if (pdfDataPct)
*pdfDataPct = 100.0 * (1.0 - poPolyNonCoveredBySources->get_Area() /
nXSize / nYSize);
}
return nStatus;
}
#endif // HAVE_GEOS

/************************************************************************/
/* GetMetadataDomainList() */
/************************************************************************/
Expand Down

0 comments on commit 15866aa

Please sign in to comment.