[geos-commits] [SCM] GEOS branch main updated. a13c2aab335ac237b89e0748cf74e0c15c8498cc
git at osgeo.org
git at osgeo.org
Wed May 22 17:12:51 PDT 2024
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, main has been updated
via a13c2aab335ac237b89e0748cf74e0c15c8498cc (commit)
from 19455699c7fe48ae54d0bdb4fa914bf40b46e247 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit a13c2aab335ac237b89e0748cf74e0c15c8498cc
Author: vadimcn <vadimcn at gmail.com>
Date: Wed May 22 17:12:28 2024 -0700
GEOSClipByRect: preserve Z coordinate (#1095)
For [multi]linestrings/linerings, Z coordinate is interpolated
from that of the nearest existing vertices.
For [multi]polygons, Z coordinate of points created on the boundary
is interpolated as above. For interior points, however, it is populated
using the ElevationModel.
diff --git a/src/operation/intersection/RectangleIntersection.cpp b/src/operation/intersection/RectangleIntersection.cpp
index e9ddc3142..e3f7750af 100644
--- a/src/operation/intersection/RectangleIntersection.cpp
+++ b/src/operation/intersection/RectangleIntersection.cpp
@@ -17,6 +17,7 @@
#include <geos/operation/intersection/RectangleIntersection.h>
#include <geos/operation/intersection/Rectangle.h>
#include <geos/operation/intersection/RectangleIntersectionBuilder.h>
+#include <geos/operation/overlayng/ElevationModel.h>
#include <geos/operation/predicate/RectangleIntersects.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/CoordinateSequence.h>
@@ -34,6 +35,7 @@
using geos::operation::intersection::Rectangle;
using geos::operation::intersection::RectangleIntersectionBuilder;
+using geos::operation::overlayng::ElevationModel;
using namespace geos::geom;
using namespace geos::algorithm;
namespace geos {
@@ -62,15 +64,18 @@ different(double x1, double y1, double x2, double y2)
inline
void
-clip_one_edge(double& x1, double& y1, double x2, double y2, double limit)
+clip_one_edge(double& x1, double& y1, double& z1, double x2, double y2, double z2, double limit)
{
if(x2 == limit) {
y1 = y2;
x1 = x2;
+ z1 = z2;
}
if(x1 != x2) {
- y1 += (y2 - y1) * (limit - x1) / (x2 - x1);
+ double fraction = (limit - x1) / (x2 - x1);
+ y1 += (y2 - y1) * fraction;
+ z1 += (z2 - z1) * fraction;
x1 = limit;
}
}
@@ -86,22 +91,22 @@ clip_one_edge(double& x1, double& y1, double x2, double y2, double limit)
*/
void
-clip_to_edges(double& x1, double& y1,
- double x2, double y2,
+clip_to_edges(double& x1, double& y1, double& z1,
+ double x2, double y2, double z2,
const Rectangle& rect)
{
if(x1 < rect.xmin()) {
- clip_one_edge(x1, y1, x2, y2, rect.xmin());
+ clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmin());
}
else if(x1 > rect.xmax()) {
- clip_one_edge(x1, y1, x2, y2, rect.xmax());
+ clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmax());
}
if(y1 < rect.ymin()) {
- clip_one_edge(y1, x1, y2, x2, rect.ymin());
+ clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymin());
}
else if(y1 > rect.ymax()) {
- clip_one_edge(y1, x1, y2, x2, rect.ymax());
+ clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymax());
}
}
@@ -153,6 +158,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
double x0 = 0;
double y0 = 0;
+ double z0 = 0;
bool add_start = false;
// Start iterating
@@ -164,6 +170,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
double x = cs[i].x;
double y = cs[i].y;
+ double z = cs[i].z;
Rectangle::Position pos = rect.position(x, y);
if(pos == Rectangle::Outside) {
@@ -199,12 +206,14 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
// Establish new position
x = cs[i].x;
y = cs[i].y;
+ z = cs[i].z;
pos = rect.position(x, y);
// Handle all possible cases
x0 = cs[i - 1].x;
y0 = cs[i - 1].y;
- clip_to_edges(x0, y0, x, y, rect);
+ z0 = cs[i - 1].z;
+ clip_to_edges(x0, y0, z0, x, y, z, rect);
if(pos == Rectangle::Inside) {
add_start = true; // x0,y0 must have clipped the rectangle
@@ -218,7 +227,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
// which will then enter the Outside section.
// Clip the other end too
- clip_to_edges(x, y, x0, y0, rect);
+ clip_to_edges(x, y, z, x0, y0, z0, rect);
Rectangle::Position prev_pos = rect.position(x0, y0);
pos = rect.position(x, y);
@@ -229,8 +238,8 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
!Rectangle::onSameEdge(prev_pos, pos) // discard if travels along edge
) {
auto coords = detail::make_unique<geom::CoordinateSequence>(2u);
- coords->setAt(Coordinate(x0, y0), 0);
- coords->setAt(Coordinate(x, y), 1);
+ coords->setAt(Coordinate(x0, y0, z0), 0);
+ coords->setAt(Coordinate(x, y, z), 1);
auto line = _gf->createLineString(std::move(coords));
parts.add(line.release());
}
@@ -267,6 +276,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
while(!go_outside && ++i < n) {
x = cs[i].x;
y = cs[i].y;
+ z = cs[i].z;
Rectangle::Position prev_pos = pos;
pos = rect.position(x, y);
@@ -278,7 +288,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
go_outside = true;
// Clip the outside point to edges
- clip_to_edges(x, y, cs[i - 1].x, cs[i - 1].y, rect);
+ clip_to_edges(x, y, z, cs[i - 1].x, cs[i - 1].y, cs[i - 1].z, rect);
pos = rect.position(x, y);
// Does the line exit through the inside of the box?
@@ -291,7 +301,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
if(start_index < i - 1 || add_start || through_box) {
auto coords = detail::make_unique<CoordinateSequence>();
if(add_start) {
- coords->add(Coordinate(x0, y0));
+ coords->add(Coordinate(x0, y0, z0));
add_start = false;
}
//line->addSubLineString(&g, start_index, i-1);
@@ -299,7 +309,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
cs.begin() + static_cast<long>(i));
if(through_box) {
- coords->add(Coordinate(x, y));
+ coords->add(Coordinate(x, y, z));
}
auto line = _gf->createLineString(std::move(coords));
@@ -316,7 +326,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
//geom::LineString * line = new geom::LineString();
if(add_start) {
//line->addPoint(x0,y0);
- coords->add(Coordinate(x0, y0));
+ coords->add(Coordinate(x0, y0, z0));
add_start = false;
}
//line->addSubLineString(&g, start_index, i-1);
@@ -349,7 +359,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
//geom::LineString * line = new geom::LineString();
if(add_start) {
//line->addPoint(x0,y0);
- coords->add(Coordinate(x0, y0));
+ coords->add(Coordinate(x0, y0, z0));
add_start = false;
}
//line->addSubLineString(&g, start_index, i-1);
@@ -521,7 +531,6 @@ RectangleIntersection::clip_polygon_to_polygons(const geom::Polygon* g,
parts.reconnectPolygons(rect);
parts.release(toParts);
-
}
/**
@@ -682,7 +691,14 @@ std::unique_ptr<geom::Geometry>
RectangleIntersection::clip(const geom::Geometry& g, const Rectangle& rect)
{
RectangleIntersection ri(g, rect);
- return ri.clip();
+ std::unique_ptr<geom::Geometry> result = ri.clip();
+
+ if (g.hasZ()) {
+ std::unique_ptr<ElevationModel> elevModel = ElevationModel::create(g);
+ elevModel->populateZ(*result);
+ }
+
+ return result;
}
std::unique_ptr<geom::Geometry>
diff --git a/tests/unit/capi/GEOSClipByRectTest.cpp b/tests/unit/capi/GEOSClipByRectTest.cpp
index 13fde6f73..be47a4c0f 100644
--- a/tests/unit/capi/GEOSClipByRectTest.cpp
+++ b/tests/unit/capi/GEOSClipByRectTest.cpp
@@ -42,6 +42,27 @@ struct test_capigeosclipbyrect_data : public capitest::utility {
ensure(equal);
}
+ void
+ checkClipByRectIdentical(
+ const char* wktInput,
+ double xmin, double ymin,
+ double xmax, double ymax,
+ const char* wktExpected)
+ {
+ input_ = fromWKT(wktInput);
+ result_ = GEOSClipByRect(input_, xmin, ymin, xmax, ymax);
+ expected_ = fromWKT(wktExpected);
+
+ bool equal = GEOSEqualsIdentical(result_, expected_) == 1;
+ if (!equal) {
+ wkt_ = GEOSWKTWriter_write(wktw_, expected_);
+ std::printf("EXP: %s\n", wkt_);
+ GEOSFree(wkt_);
+ wkt_ = GEOSWKTWriter_write(wktw_, result_);
+ std::printf("OBT: %s\n", wkt_);
+ }
+ ensure(equal);
+ }
};
typedef test_group<test_capigeosclipbyrect_data> group;
@@ -70,7 +91,7 @@ template<>
template<>
void object::test<2>()
{
- checkClipByRect(
+ checkClipByRectIdentical(
"POINT(15 15)",
10, 10, 20, 20,
"POINT(15 15)"
@@ -106,7 +127,7 @@ template<>
template<>
void object::test<5>()
{
- checkClipByRect(
+ checkClipByRectIdentical(
"LINESTRING(15 15, 16 15)",
10, 10, 20, 20,
"LINESTRING(15 15, 16 15)"
@@ -130,7 +151,7 @@ template<>
template<>
void object::test<7>()
{
- checkClipByRect(
+ checkClipByRectIdentical(
"LINESTRING(10 5, 25 20)",
10, 10, 20, 20,
"LINESTRING (15 10, 20 15)"
@@ -191,7 +212,7 @@ template<>
void object::test<12>()
{
const char* wkt = "POLYGON((1 1, 1 30, 30 30, 30 1, 1 1),(10 10, 20 10, 20 20, 10 20, 10 10))";
- checkClipByRect(
+ checkClipByRectIdentical(
wkt,
0, 0, 40, 40,
wkt);
@@ -202,7 +223,7 @@ template<>
template<>
void object::test<13>()
{
- checkClipByRect(
+ checkClipByRectIdentical(
"POLYGON((0 0, 0 30, 30 30, 30 0, 0 0),(10 10, 20 10, 20 20, 10 20, 10 10))",
5, 5, 15, 15,
"POLYGON ((5 5, 5 15, 10 15, 10 10, 15 10, 15 5, 5 5))"
@@ -282,6 +303,66 @@ void object::test<17>()
ensure(result_ == nullptr);
}
+/// Point Z inside
+template<>
+template<>
+void object::test<18>()
+{
+ checkClipByRectIdentical(
+ "POINT Z (15 15 100)",
+ 10, 10, 20, 20,
+ "POINT Z (15 15 100)"
+ );
+}
+
+/// Line Z inside
+template<>
+template<>
+void object::test<19>()
+{
+ checkClipByRectIdentical(
+ "LINESTRING Z (15 15 0, 16 15 100)",
+ 10, 10, 20, 20,
+ "LINESTRING Z (15 15 0, 16 15 100)"
+ );
+}
+
+/// Line Z splitting rectangle
+template<>
+template<>
+void object::test<20>()
+{
+ checkClipByRectIdentical(
+ "LINESTRING Z (0 15 0, 100 15 100)",
+ 10, 10, 20, 20,
+ "LINESTRING Z (10 15 10, 20 15 20)"
+ );
+}
+
+/// Polygon Z overlapping rectangle
+template<>
+template<>
+void object::test<21>()
+{
+ checkClipByRectIdentical(
+ "POLYGON Z ((0 0 100, 0 30 100, 30 30 100, 30 0 100, 0 0 100),(10 10 100, 20 10 100, 20 20 100, 10 20 100, 10 10 100))",
+ 5, 5, 15, 15,
+ "POLYGON Z ((5 5 100, 5 15 100, 10 15 100, 10 10 100, 15 10 100, 15 5 100, 5 5 100))"
+ );
+}
+
+/// Polygon Z enclosing rectangle
+template<>
+template<>
+void object::test<22>()
+{
+ checkClipByRectIdentical(
+ "POLYGON Z ((0 0 100, 0 30 100, 30 30 100, 30 0 100, 0 0 100))",
+ 5, 5, 15, 15,
+ "POLYGON Z ((5 5 100, 5 15 100, 15 15 100, 15 5 100, 5 5 100))"
+ );
+}
+
} // namespace tut
-----------------------------------------------------------------------
Summary of changes:
.../intersection/RectangleIntersection.cpp | 54 ++++++++-----
tests/unit/capi/GEOSClipByRectTest.cpp | 91 ++++++++++++++++++++--
2 files changed, 121 insertions(+), 24 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list