Skip to content

Commit

Permalink
ogr2ogr: fix -clipsrc/-clipdst when a input geometry is within the en…
Browse files Browse the repository at this point in the history 
…velope of a non-rectangle clip geometry, but doesn't intersect it (3.10.0 regression)

Fixes #11652
Fixes #10341

This fixes commit
388d3ae (PR #10341) "ogr2ogr: speed-up
-clipsrc/-clipdst by avoiding GEOS when possible", which only worked if
the clipping geometry was a rectangle.
  • Loading branch information
@rouault @github-actions
rouault authored and github-actions[bot] committed Jan 15, 2025
1 parent bfd01f5 commit 9d52f68
Show file tree
Hide file tree
Showing 2 changed files with 96 additions and 27 deletions.
75 changes: 48 additions & 27 deletions apps/ogr2ogr_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -581,12 +581,14 @@ class LayerTranslator
std::unique_ptr<OGRGeometry> m_poClipSrcReprojectedToSrcSRS{};
const OGRSpatialReference *m_poClipSrcReprojectedToSrcSRS_SRS = nullptr;
OGREnvelope m_oClipSrcEnv{};
bool m_bClipSrcIsRectangle = false;

OGRGeometry *m_poClipDstOri = nullptr;
bool m_bWarnedClipDstSRS = false;
std::unique_ptr<OGRGeometry> m_poClipDstReprojectedToDstSRS{};
const OGRSpatialReference *m_poClipDstReprojectedToDstSRS_SRS = nullptr;
OGREnvelope m_oClipDstEnv{};
bool m_bClipDstIsRectangle = false;

bool m_bExplodeCollections = false;
bool m_bNativeData = false;
Expand All @@ -600,10 +602,15 @@ class LayerTranslator
const GDALVectorTranslateOptions *psOptions);

private:
std::pair<const OGRGeometry *, const OGREnvelope *>
GetDstClipGeom(const OGRSpatialReference *poGeomSRS);
std::pair<const OGRGeometry *, const OGREnvelope *>
GetSrcClipGeom(const OGRSpatialReference *poGeomSRS);
struct ClipGeomDesc
{
const OGRGeometry *poGeom = nullptr;
const OGREnvelope *poEnv = nullptr;
bool bGeomIsRectangle = false;
};

ClipGeomDesc GetDstClipGeom(const OGRSpatialReference *poGeomSRS);
ClipGeomDesc GetSrcClipGeom(const OGRSpatialReference *poGeomSRS);
};

static OGRLayer *GetLayerAndOverwriteIfNecessary(GDALDataset *poDstDS,
Expand Down Expand Up @@ -6317,16 +6324,17 @@ bool LayerTranslator::Translate(
if (poStolenGeometry->IsEmpty())
goto end_loop;

const auto [poClipGeom, poClipGeomEnvelope] =
const auto clipGeomDesc =
GetSrcClipGeom(poStolenGeometry->getSpatialReference());

if (poClipGeom && poClipGeomEnvelope)
if (clipGeomDesc.poGeom && clipGeomDesc.poEnv)
{
OGREnvelope oEnv;
poStolenGeometry->getEnvelope(&oEnv);
if (!poClipGeomEnvelope->Contains(oEnv) &&
!(poClipGeomEnvelope->Intersects(oEnv) &&
poClipGeom->Intersects(poStolenGeometry.get())))
if (!clipGeomDesc.poEnv->Contains(oEnv) &&
!(clipGeomDesc.poEnv->Intersects(oEnv) &&
clipGeomDesc.poGeom->Intersects(
poStolenGeometry.get())))
{
goto end_loop;
}
Expand Down Expand Up @@ -6568,22 +6576,23 @@ bool LayerTranslator::Translate(
if (poDstGeometry->IsEmpty())
goto end_loop;

const auto [poClipGeom, poClipGeomEnvelope] =
const auto clipGeomDesc =
GetSrcClipGeom(poDstGeometry->getSpatialReference());

if (!(poClipGeom && poClipGeomEnvelope))
if (!(clipGeomDesc.poGeom && clipGeomDesc.poEnv))
goto end_loop;

OGREnvelope oDstEnv;
poDstGeometry->getEnvelope(&oDstEnv);

if (!poClipGeomEnvelope->Contains(oDstEnv))
if (!(clipGeomDesc.bGeomIsRectangle &&
clipGeomDesc.poEnv->Contains(oDstEnv)))
{
std::unique_ptr<OGRGeometry> poClipped;
if (poClipGeomEnvelope->Intersects(oDstEnv))
if (clipGeomDesc.poEnv->Intersects(oDstEnv))
{
poClipped.reset(
poClipGeom->Intersection(poDstGeometry.get()));
poClipped.reset(clipGeomDesc.poGeom->Intersection(
poDstGeometry.get()));
}
if (poClipped == nullptr || poClipped->IsEmpty())
{
Expand Down Expand Up @@ -6750,23 +6759,25 @@ bool LayerTranslator::Translate(
if (poDstGeometry->IsEmpty())
goto end_loop;

auto [poClipGeom, poClipGeomEnvelope] = GetDstClipGeom(
const auto clipGeomDesc = GetDstClipGeom(
poDstGeometry->getSpatialReference());
if (!poClipGeom || !poClipGeomEnvelope)
if (!clipGeomDesc.poGeom || !clipGeomDesc.poEnv)
{
goto end_loop;
}

OGREnvelope oDstEnv;
poDstGeometry->getEnvelope(&oDstEnv);

if (!poClipGeomEnvelope->Contains(oDstEnv))
if (!(clipGeomDesc.bGeomIsRectangle &&
clipGeomDesc.poEnv->Contains(oDstEnv)))
{
std::unique_ptr<OGRGeometry> poClipped;
if (poClipGeomEnvelope->Intersects(oDstEnv))
if (clipGeomDesc.poEnv->Intersects(oDstEnv))
{
poClipped.reset(poClipGeom->Intersection(
poDstGeometry.get()));
poClipped.reset(
clipGeomDesc.poGeom->Intersection(
poDstGeometry.get()));
}

if (poClipped == nullptr || poClipped->IsEmpty())
Expand Down Expand Up @@ -6975,7 +6986,7 @@ bool LayerTranslator::Translate(
* expressed.
* @return the destination clip geometry and its envelope, or (nullptr, nullptr)
*/
std::pair<const OGRGeometry *, const OGREnvelope *>
LayerTranslator::ClipGeomDesc
LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
{
if (m_poClipDstReprojectedToDstSRS_SRS != poGeomSRS)
Expand All @@ -6988,7 +6999,7 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
if (m_poClipDstReprojectedToDstSRS->transformTo(poGeomSRS) !=
OGRERR_NONE)
{
return std::make_pair(nullptr, nullptr);
return ClipGeomDesc();
}
m_poClipDstReprojectedToDstSRS_SRS = poGeomSRS;
}
Expand All @@ -7014,8 +7025,13 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
if (poGeom && !m_oClipDstEnv.IsInit())
{
poGeom->getEnvelope(&m_oClipDstEnv);
m_bClipDstIsRectangle = poGeom->IsRectangle();
}
return std::make_pair(poGeom, poGeom ? &m_oClipDstEnv : nullptr);
ClipGeomDesc ret;
ret.poGeom = poGeom;
ret.poEnv = poGeom ? &m_oClipDstEnv : nullptr;
ret.bGeomIsRectangle = m_bClipDstIsRectangle;
return ret;
}

/************************************************************************/
Expand All @@ -7028,7 +7044,7 @@ LayerTranslator::GetDstClipGeom(const OGRSpatialReference *poGeomSRS)
* expressed.
* @return the source clip geometry and its envelope, or (nullptr, nullptr)
*/
std::pair<const OGRGeometry *, const OGREnvelope *>
LayerTranslator::ClipGeomDesc
LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS)
{
if (m_poClipSrcReprojectedToSrcSRS_SRS != poGeomSRS)
Expand All @@ -7041,7 +7057,7 @@ LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS)
if (m_poClipSrcReprojectedToSrcSRS->transformTo(poGeomSRS) !=
OGRERR_NONE)
{
return std::make_pair(nullptr, nullptr);
return ClipGeomDesc();
}
m_poClipSrcReprojectedToSrcSRS_SRS = poGeomSRS;
}
Expand All @@ -7066,8 +7082,13 @@ LayerTranslator::GetSrcClipGeom(const OGRSpatialReference *poGeomSRS)
if (poGeom && !m_oClipSrcEnv.IsInit())
{
poGeom->getEnvelope(&m_oClipSrcEnv);
m_bClipSrcIsRectangle = poGeom->IsRectangle();
}
return std::make_pair(poGeom, poGeom ? &m_oClipSrcEnv : nullptr);
ClipGeomDesc ret;
ret.poGeom = poGeom;
ret.poEnv = poGeom ? &m_oClipSrcEnv : nullptr;
ret.bGeomIsRectangle = m_bClipDstIsRectangle;
return ret;
}

/************************************************************************/
Expand Down
48 changes: 48 additions & 0 deletions autotest/utilities/test_ogr2ogr_lib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1116,6 +1116,17 @@ def test_ogr2ogr_lib_clipsrc_datasource(tmp_vsimem):
f.SetField("filter_field", "exact_overlap_full_result")
f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))"))
clip_layer.CreateFeature(f)
# Clip geometry envelope contains envelope of input geometry, but does not intersect it
f = ogr.Feature(clip_layer.GetLayerDefn())
f.SetField(
"filter_field", "clip_geometry_envelope_contains_envelope_but_no_intersect"
)
f.SetGeometry(
ogr.CreateGeometryFromWkt(
"POLYGON ((-2 -1,-2 4,4 4,4 -1,3 -1,3 3,-1 3,-1 -1,-2 -1))"
)
)
clip_layer.CreateFeature(f)
clip_ds = None

# Test clip with 'half_overlap_line_result' using sql statement
Expand Down Expand Up @@ -1162,6 +1173,19 @@ def test_ogr2ogr_lib_clipsrc_datasource(tmp_vsimem):
dst_ds = None
gdal.Unlink(dst_filename)

# Test clip with the "clip_geometry_envelope_contains_envelope_but_no_intersect" using only clipSrcWhere
dst_ds = gdal.VectorTranslate(
dst_filename,
src_filename,
format="GPKG",
clipSrc=clip_path,
clipSrcWhere="filter_field = 'clip_geometry_envelope_contains_envelope_but_no_intersect'",
)
dst_lyr = dst_ds.GetLayer(0)
assert dst_lyr.GetFeatureCount() == 0
dst_ds = None
gdal.Unlink(dst_filename)

# Cleanup
gdal.Unlink(clip_path)

Expand Down Expand Up @@ -1390,6 +1414,17 @@ def test_ogr2ogr_lib_clipdst_datasource(tmp_vsimem):
f.SetField("filter_field", "exact_overlap_full_result")
f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0, 0 2, 2 2, 2 0, 0 0))"))
clip_layer.CreateFeature(f)
# Clip geometry envelope contains envelope of input geometry, but does not intersect it
f = ogr.Feature(clip_layer.GetLayerDefn())
f.SetField(
"filter_field", "clip_geometry_envelope_contains_envelope_but_no_intersect"
)
f.SetGeometry(
ogr.CreateGeometryFromWkt(
"POLYGON ((-2 -1,-2 4,4 4,4 -1,3 -1,3 3,-1 3,-1 -1,-2 -1))"
)
)
clip_layer.CreateFeature(f)
clip_ds = None

# Test clip with 'half_overlap_line_result' using sql statement
Expand Down Expand Up @@ -1436,6 +1471,19 @@ def test_ogr2ogr_lib_clipdst_datasource(tmp_vsimem):
dst_ds = None
gdal.Unlink(dst_filename)

# Test clip with the "clip_geometry_envelope_contains_envelope_but_no_intersect" using only clipSrcWhere
dst_ds = gdal.VectorTranslate(
dst_filename,
src_filename,
format="GPKG",
clipDst=clip_path,
clipDstWhere="filter_field = 'clip_geometry_envelope_contains_envelope_but_no_intersect'",
)
dst_lyr = dst_ds.GetLayer(0)
assert dst_lyr.GetFeatureCount() == 0
dst_ds = None
gdal.Unlink(dst_filename)

# Cleanup
gdal.Unlink(clip_path)

Expand Down

0 comments on commit 9d52f68

Please sign in to comment.