[geos-commits] [SCM] GEOS branch master updated. 2ad9efe1479beccf1fc204094f5fb8d39a82053e
git at osgeo.org
git at osgeo.org
Thu Nov 26 14:05:33 PST 2020
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, master has been updated
via 2ad9efe1479beccf1fc204094f5fb8d39a82053e (commit)
from 7d4e17336f45d39aa2debd91237eed461aa5da09 (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 2ad9efe1479beccf1fc204094f5fb8d39a82053e
Author: Sandro Santilli <strk at kbt.io>
Date: Mon Nov 23 13:19:56 2020 +0100
Port ElevationModel class from JTS
Closes #1070
Differences from JTS:
- Cells are all initialized upfront, no heap allocations
- CoordinateSequence::hasZ implemented as non-virtual
- ElevationModel::populateZ uses CoordinateFilter instead
of CoordinateSequenceFilter (see note 1)
Note1:
Using CoordinateFilter instead of CoordinateSequenceFilter
allows CoordinateSequence to clear the dimension cache which
would otherwise not properly set the sequence as having Z
(see #435)
diff --git a/include/geos/geom/CoordinateSequence.h b/include/geos/geom/CoordinateSequence.h
index 976b2cc..e8eeca5 100644
--- a/include/geos/geom/CoordinateSequence.h
+++ b/include/geos/geom/CoordinateSequence.h
@@ -222,6 +222,10 @@ public:
*/
virtual std::size_t getDimension() const = 0;
+ bool hasZ() const {
+ return getDimension() > 2;
+ }
+
/**
* Returns the ordinate of a coordinate in this sequence.
* Ordinate indices 0 and 1 are assumed to be X and Y.
diff --git a/include/geos/operation/overlayng/ElevationModel.h b/include/geos/operation/overlayng/ElevationModel.h
new file mode 100644
index 0000000..199efbc
--- /dev/null
+++ b/include/geos/operation/overlayng/ElevationModel.h
@@ -0,0 +1,170 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Sandro Santilli <strk at kbt.io>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last Port: operation/overlayng/ElevationModel.java a0e53557
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+
+#include <geos/geom/Envelope.h> // for composition
+
+// Forward declarations
+namespace geos {
+ namespace geom {
+ class Geometry;
+ }
+}
+
+namespace geos { // geos.
+namespace operation { // geos.operation
+namespace overlayng { // geos.operation.overlayng
+
+
+/**
+ * \brief
+ * A simple elevation model used to populate missing Z values
+ * in overlay results.
+ *
+ * The model divides the extent of the input geometry(s)
+ * into an NxM grid.
+ * The default grid size is 3x3.
+ * If the input has no extent in the X or Y dimension,
+ * that dimension is given grid size 1.
+ * The elevation of each grid cell is computed as the average of the Z
+ * values
+ * of the input vertices in that cell (if any).
+ * If a cell has no input vertices within it, it is assigned
+ * the average elevation over all cells.
+ *
+ * If no input vertices have Z values, the model does not assign a Z
+ * value.
+ *
+ * The elevation of an arbitrary location is determined as the
+ * Z value of the nearest grid cell.
+ *
+ * An elevation model can be used to populate missing Z values
+ * in an overlay result geometry.
+ *
+ * @author Martin Davis
+ *
+ */
+class GEOS_DLL ElevationModel {
+
+private:
+
+ class ElevationCell {
+ private:
+
+ int numZ = 0;
+ double sumZ = 0.0;
+ double avgZ;
+
+ public:
+
+ bool isNull() const
+ {
+ return numZ == 0;
+ }
+
+ void add(double z)
+ {
+ ++numZ;
+ sumZ += z;
+ }
+
+ void compute()
+ {
+ avgZ = DoubleNotANumber;
+ if (numZ > 0)
+ avgZ = sumZ / numZ;
+ }
+
+ double getZ() const
+ {
+ return avgZ;
+ }
+ };
+
+
+ static const int DEFAULT_CELL_NUM;
+ geom::Envelope extent;
+ int numCellX;
+ int numCellY;
+ double cellSizeX;
+ double cellSizeY;
+ std::vector<ElevationCell> cells;
+ bool isInitialized = false;
+ bool hasZValue = false;
+ double averageZ = DoubleNotANumber;
+
+ void init();
+
+ ElevationCell& getCell(double x, double y); //, bool isCreateIfMissing);
+
+ int getCellOffset(int ix, int iy) {
+ return (numCellX * iy + ix);
+ }
+
+protected:
+
+ void add(double x, double y, double z);
+
+
+public:
+
+ static std::unique_ptr<ElevationModel> create(const geom::Geometry& geom1,
+ const geom::Geometry& geom2);
+
+ static std::unique_ptr<ElevationModel> create(const geom::Geometry& geom1);
+
+ ElevationModel(const geom::Envelope& extent, int numCellX, int numCellY);
+
+ void add(const geom::Geometry& geom);
+
+
+ /**
+ * Gets the model Z value at a given location.
+ * If the location lies outside the model grid extent,
+ * this returns the Z value of the nearest grid cell.
+ * If the model has no elevation computed (i.e. due
+ * to empty input), the value is returned as {@link Double#NaN}.
+ *
+ * @param x the x ordinate of the location
+ * @param y the y ordinate of the location
+ * @return the computed model Z value
+ */
+ double getZ(double x, double y);
+
+
+ /**
+ * \brief
+ * Computes Z values for any missing Z values in a geometry,
+ * using the computed model.
+ *
+ * If the model has no Z value, or the geometry coordinate dimension
+ * does not include Z, no action is taken.
+ *
+ * @param geom the geometry to elevate
+ */
+ void populateZ(geom::Geometry& geom);
+
+
+};
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
diff --git a/include/geos/operation/overlayng/Makefile.am b/include/geos/operation/overlayng/Makefile.am
index 0dd90d8..8efcceb 100644
--- a/include/geos/operation/overlayng/Makefile.am
+++ b/include/geos/operation/overlayng/Makefile.am
@@ -13,6 +13,7 @@ geos_HEADERS = \
EdgeMerger.h \
EdgeNodingBuilder.h \
EdgeSourceInfo.h \
+ ElevationModel.h \
IndexedPointOnLineLocator.h \
InputGeometry.h \
IntersectionPointBuilder.h \
diff --git a/include/geos/operation/overlayng/OverlayNG.h b/include/geos/operation/overlayng/OverlayNG.h
index 505e918..30219e4 100644
--- a/include/geos/operation/overlayng/OverlayNG.h
+++ b/include/geos/operation/overlayng/OverlayNG.h
@@ -3,6 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
+ * Copyright (C) 2020 Sandro Santilli <strk at kbt.io>
* Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
*
* This is free software; you can redistribute and/or modify it under
@@ -10,6 +11,10 @@
* by the Free Software Foundation.
* See the COPYING file for more information.
*
+ **********************************************************************
+ *
+ * Last Port: operation/overlayng/OverlayNG.java 459bc7f8e
+ *
**********************************************************************/
#pragma once
diff --git a/src/operation/overlayng/ElevationModel.cpp b/src/operation/overlayng/ElevationModel.cpp
new file mode 100644
index 0000000..d672dee
--- /dev/null
+++ b/src/operation/overlayng/ElevationModel.cpp
@@ -0,0 +1,277 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Sandro Santilli <strk at kbt.io>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last Port: operation/overlayng/ElevationModel.java 459bc7f8
+ *
+ **********************************************************************/
+
+#include <geos/operation/overlayng/ElevationModel.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateFilter.h>
+#include <geos/geom/CoordinateSequenceFilter.h>
+#include <geos/geom/Envelope.h>
+
+#include <memory> // for unique_ptr
+
+#ifndef GEOS_DEBUG
+# define GEOS_DEBUG 0
+#endif
+
+#if GEOS_DEBUG
+# include <iostream>
+#endif
+
+namespace geos { // geos.
+namespace operation { // geos.operation
+namespace overlayng { // geos.operation.overlayng
+
+namespace {
+ template<typename T>
+ T clamp(const T& v, const T& low, const T& high) {
+ return v < low ? low : high < v ? high : v;
+ }
+}
+
+using geos::geom::Coordinate;
+using geos::geom::CoordinateSequence;
+using geos::geom::Geometry;
+using geos::geom::Envelope;
+
+/* public static */
+std::unique_ptr<ElevationModel>
+ElevationModel::create(const geom::Geometry& geom1, const geom::Geometry& geom2)
+{
+ Envelope extent;
+ if (! geom1.isEmpty() ) {
+ extent.expandToInclude(geom1.getEnvelopeInternal());
+ }
+ if (! geom2.isEmpty() ) {
+ extent.expandToInclude(geom2.getEnvelopeInternal());
+ }
+ std::unique_ptr<ElevationModel> model ( new ElevationModel(extent, DEFAULT_CELL_NUM, DEFAULT_CELL_NUM) );
+ if (! geom1.isEmpty() ) model->add(geom1);
+ if (! geom2.isEmpty() ) model->add(geom2);
+ return model;
+}
+
+/* public static */
+std::unique_ptr<ElevationModel>
+ElevationModel::create(const geom::Geometry& geom1)
+{
+ Envelope extent;
+ if (! geom1.isEmpty() ) {
+ extent.expandToInclude(geom1.getEnvelopeInternal());
+ }
+ std::unique_ptr<ElevationModel> model ( new ElevationModel(extent, DEFAULT_CELL_NUM, DEFAULT_CELL_NUM) );
+ if (! geom1.isEmpty() ) model->add(geom1);
+ return model;
+}
+
+const int ElevationModel::DEFAULT_CELL_NUM = 3;
+
+ElevationModel::ElevationModel(const Envelope& nExtent, int nNumCellX, int nNumCellY)
+ :
+ extent(nExtent),
+ numCellX(nNumCellX),
+ numCellY(nNumCellY)
+{
+
+ cellSizeX = extent.getWidth() / numCellX;
+ cellSizeY = extent.getHeight() / numCellY;
+ if(cellSizeX <= 0.0) {
+ numCellX = 1;
+ }
+ if(cellSizeY <= 0.0) {
+ numCellY = 1;
+ }
+ cells.resize(numCellX * numCellY);
+}
+
+/* public */
+void
+ElevationModel::add(const Geometry& geom)
+{
+ class Filter: public geom::CoordinateSequenceFilter {
+ ElevationModel& model;
+ bool hasZ;
+ public:
+ Filter(ElevationModel& nModel) : model(nModel), hasZ(true) {}
+
+ void filter_ro(const geom::CoordinateSequence& seq, std::size_t i) override
+ {
+#if 0
+ if (! seq.hasZ()) {
+ hasZ = false;;
+ return;
+ }
+#endif
+ const Coordinate& c = seq.getAt(i);
+#if GEOS_DEBUG
+ std::cout << "Coordinate " << i << " of added geom is: "
+<< c << std::endl;
+#endif
+ model.add(c.x, c.y, c.z);
+ }
+
+ bool isDone() const override {
+ // no need to scan if no Z present
+ return ! hasZ;
+ }
+
+ bool isGeometryChanged() const override {
+ return false;
+ }
+
+ };
+
+ Filter filter(*this);
+ geom.apply_ro(filter);
+
+}
+
+/* protected */
+void
+ElevationModel::add(double x, double y, double z)
+{
+ if (std::isnan(z))
+ return;
+ hasZValue = true;
+ ElevationCell& cell = getCell(x, y); //, true);
+ cell.add(z);
+}
+
+/* private */
+void
+ElevationModel::init()
+{
+ isInitialized = true;
+ int numCells = 0;
+ double sumZ = 0.0;
+
+#if GEOS_DEBUG
+ int offset = 0;
+#endif
+ for (ElevationCell& cell : cells )
+ {
+ if (!cell.isNull()) {
+ cell.compute();
+#if GEOS_DEBUG
+ std::cout << "init: cell" << offset
+ << " getZ: " << cell.getZ()
+ << std::endl;
+#endif
+ numCells++;
+ sumZ += cell.getZ();
+ }
+#if GEOS_DEBUG
+ ++offset;
+#endif
+ }
+ averageZ = DoubleNotANumber;
+ if (numCells > 0) {
+ averageZ = sumZ / numCells;
+ }
+#if GEOS_DEBUG
+ std::cout << "init: numCells: " << numCells
+ << " averageZ: " << averageZ << std::endl;
+#endif
+}
+
+/* public */
+double
+ElevationModel::getZ(double x, double y)
+{
+ if (! isInitialized)
+ {
+ init();
+ }
+ const ElevationModel::ElevationCell& cell = getCell(x, y); //, false);
+ if (cell.isNull())
+ {
+ return averageZ;
+ }
+ return cell.getZ();
+}
+
+/* public */
+void
+ElevationModel::populateZ(Geometry& geom)
+{
+ // short-circuit if no Zs are present in model
+ if (! hasZValue)
+ return;
+
+ if (! isInitialized)
+ init();
+
+ class Filter: public geom::CoordinateFilter {
+ ElevationModel& model;
+ public:
+ Filter(ElevationModel& nModel) : model(nModel) {}
+
+ void filter_rw(geom::Coordinate* c) const override
+ {
+#if GEOS_DEBUG
+ std::cout << "Input coord: " << *c << std::endl;
+#endif
+ if (std::isnan( c->z )) {
+ c->z = model.getZ(c->x, c->y);
+#if GEOS_DEBUG
+ std::cout << "New coord: " << *c << std::endl;
+#endif
+ }
+ }
+
+ };
+
+ Filter filter(*this);
+#if GEOS_DEBUG
+ std::cout << "ElevationModel::poplateZ calling apply_rw(CoordinateSequenceFilter&) against" <<
+ //std::type_name<decltype(ci)>()
+ typeid(geom).name()
+ << std::endl;
+#endif
+ geom.apply_rw(&filter);
+}
+
+/* private */
+ElevationModel::ElevationCell&
+ElevationModel::getCell(double x, double y) //, bool isCreateIfMissing
+{
+ int ix = 0;
+ if (numCellX > 1) {
+ ix = (int) ((x - extent.getMinX()) / cellSizeX);
+ ix = clamp(ix, 0, numCellX - 1);
+ }
+ int iy = 0;
+ if (numCellY > 1) {
+ iy = (int) ((y - extent.getMinY()) / cellSizeY);
+ iy = clamp(iy, 0, numCellY - 1);
+ }
+ int cellOffset = getCellOffset(ix, iy);
+#if GEOS_DEBUG
+ std::cout << "Cell of coordinate " << x << "," << y
+ << " is " << ix << "," << iy
+ << " offset " << cellOffset << std::endl;
+#endif
+ assert ( cellOffset < numCellX * numCellY );
+ ElevationModel::ElevationCell& cell = cells[cellOffset];
+ return cell;
+}
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
diff --git a/src/operation/overlayng/Makefile.am b/src/operation/overlayng/Makefile.am
index 205273d..ee52fbc 100644
--- a/src/operation/overlayng/Makefile.am
+++ b/src/operation/overlayng/Makefile.am
@@ -12,6 +12,7 @@ libopoverlayng_la_SOURCES = \
EdgeMerger.cpp \
EdgeNodingBuilder.cpp \
EdgeSourceInfo.cpp \
+ ElevationModel.cpp \
IndexedPointOnLineLocator.cpp \
InputGeometry.cpp \
IntersectionPointBuilder.cpp \
diff --git a/src/operation/overlayng/OverlayNG.cpp b/src/operation/overlayng/OverlayNG.cpp
index baba411..6b6ad1d 100644
--- a/src/operation/overlayng/OverlayNG.cpp
+++ b/src/operation/overlayng/OverlayNG.cpp
@@ -3,6 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
+ * Copyright (C) 2020 Sandro Santilli <strk at kbt.io>
* Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
*
* This is free software; you can redistribute and/or modify it under
@@ -10,12 +11,17 @@
* by the Free Software Foundation.
* See the COPYING file for more information.
*
+ **********************************************************************
+ *
+ * Last Port: operation/overlayng/OverlayNG.java 459bc7f8e
+ *
**********************************************************************/
#include <geos/operation/overlayng/OverlayNG.h>
#include <geos/operation/overlayng/Edge.h>
#include <geos/operation/overlayng/EdgeNodingBuilder.h>
+#include <geos/operation/overlayng/ElevationModel.h>
#include <geos/operation/overlayng/InputGeometry.h>
#include <geos/operation/overlayng/IntersectionPointBuilder.h>
#include <geos/operation/overlayng/LineBuilder.h>
@@ -32,11 +38,17 @@
#include <algorithm>
-
#ifndef GEOS_DEBUG
#define GEOS_DEBUG 0
#endif
+#if GEOS_DEBUG
+# include <iostream>
+# include <geos/io/WKTWriter.h>
+#endif
+
+
+
namespace geos { // geos
namespace operation { // geos.operation
namespace overlayng { // geos.operation.overlayng
@@ -135,21 +147,57 @@ OverlayNG::geomunion(const Geometry* geom, const PrecisionModel* pm, noding::Nod
std::unique_ptr<Geometry>
OverlayNG::getResult()
{
- if (OverlayUtil::isEmptyResult(opCode, inputGeom.getGeometry(0), inputGeom.getGeometry(1), pm)) {
+ const Geometry *ig0 = inputGeom.getGeometry(0);
+ const Geometry *ig1 = inputGeom.getGeometry(1);
+
+ if ( OverlayUtil::isEmptyResult(opCode, ig0, ig1, pm) )
+ {
return createEmptyResult();
}
- // special logic for Point-Point inputs
- if (inputGeom.isAllPoints()) {
- return OverlayPoints::overlay(opCode, inputGeom.getGeometry(0), inputGeom.getGeometry(1), pm);
+ /**
+ * The elevation model is only computed if the input geometries have Z values.
+ */
+ std::unique_ptr<ElevationModel> elevModel;
+ if ( ig1 ) {
+ elevModel = ElevationModel::create(*ig0, *ig1);
+ } else {
+ elevModel = ElevationModel::create(*ig0);
}
+ std::unique_ptr<Geometry> result;
- // special logic for Point-nonPoint inputs
- if (! inputGeom.isSingle() && inputGeom.hasPoints()) {
- return OverlayMixedPoints::overlay(opCode, inputGeom.getGeometry(0), inputGeom.getGeometry(1), pm);
+ if (inputGeom.isAllPoints()) {
+ // handle Point-Point inputs
+ result = OverlayPoints::overlay(opCode, ig0, ig1, pm);
+ }
+ else if (! inputGeom.isSingle() && inputGeom.hasPoints()) {
+ // handle Point-nonPoint inputs
+ result = OverlayMixedPoints::overlay(opCode, ig0, ig1, pm);
}
+ else {
+ // handle case where both inputs are formed of edges (Lines and Polygons)
+ result = computeEdgeOverlay();
+ }
+
+#if GEOS_DEBUG
+ io::WKTWriter w;
+ w.setOutputDimension(3);
+ w.setTrim(true);
+
+ std::cout << "Before populatingZ: " << w.write(result.get()) << std::endl;
+#endif
+
+ /**
+ * This is a no-op if the elevation model was not computed due to
+ * Z not present
+ */
+ elevModel->populateZ(*result);
+
+#if GEOS_DEBUG
+ std::cout << " After populatingZ: " << w.write(result.get()) << std::endl;
+#endif
- return computeEdgeOverlay();
+ return result;
}
diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
index 51129c8..adf5d1a 100644
--- a/tests/unit/Makefile.am
+++ b/tests/unit/Makefile.am
@@ -195,6 +195,7 @@ geos_unit_SOURCES = \
operation/overlay/validate/FuzzyPointLocatorTest.cpp \
operation/overlay/validate/OffsetPointGeneratorTest.cpp \
operation/overlay/validate/OverlayResultValidatorTest.cpp \
+ operation/overlayng/ElevationModelTest.cpp \
operation/overlayng/RingClipperTest.cpp \
operation/overlayng/LineLimiterTest.cpp \
operation/overlayng/OverlayGraphTest.cpp \
diff --git a/tests/unit/operation/overlayng/ElevationModelTest.cpp b/tests/unit/operation/overlayng/ElevationModelTest.cpp
new file mode 100644
index 0000000..347f079
--- /dev/null
+++ b/tests/unit/operation/overlayng/ElevationModelTest.cpp
@@ -0,0 +1,254 @@
+//
+// Test Suite for geos::operation::overlayng::ElevationModel class.
+//
+// Last port:
+// modules/core/src/test/java/org/locationtech/jts/operation/overlayng/ElevationModelTest.java
+// 459bc7f8e4057dd05c0576b88d1550d552bf132d
+
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/operation/overlayng/ElevationModel.h>
+#include <geos/constants.h> // for DoubleNotANumber
+
+// std
+#include <memory>
+#include <initializer_list>
+
+using namespace geos::geom;
+using namespace geos::operation::overlayng;
+using geos::io::WKTReader;
+using geos::DoubleNotANumber;
+using geos::io::WKTWriter;
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used by all tests
+struct test_overlayng_elevationmodel_data {
+
+ WKTReader r;
+ WKTWriter w;
+
+ const double TOLERANCE = 0.00001;
+
+ test_overlayng_elevationmodel_data()
+ {
+ w.setTrim(true);
+ w.setOutputDimension(3);
+ }
+
+ void checkElevation(const std::string& wkt1, const std::string& wkt2, std::initializer_list<double> ords)
+ {
+ std::unique_ptr<Geometry> g1 = r.read(wkt1);
+ std::unique_ptr<Geometry> g2 = r.read(wkt2);
+ std::unique_ptr<ElevationModel> model = ElevationModel::create(*g1, *g2);
+ unsigned int numPts = ords.size() / 3;
+ assert ( 3 * numPts == ords.size());
+ for ( std::initializer_list<double>::iterator i=ords.begin(), e=ords.end();
+ i != e; ++i )
+ {
+ double x = *i++;
+ double y = *i++;
+ double expectedZ = *i;
+ double actualZ = model->getZ(x, y);
+ ensure_distance(expectedZ, actualZ, TOLERANCE);
+ }
+ }
+
+ void checkElevation(const std::string& wkt1, std::initializer_list<double> ords)
+ {
+ checkElevation(wkt1, "POINT EMPTY", ords);
+ }
+
+ void checkElevationPopulateZ(const std::string& wkt,
+ const std::string& wktNoZ,
+ const std::string& wktZExpected)
+ {
+ std::unique_ptr<Geometry> geom = r.read(wkt);
+ std::unique_ptr<ElevationModel> model = ElevationModel::create(*geom);
+
+ std::unique_ptr<Geometry> geomNoZ = r.read(wktNoZ);
+ model->populateZ(*geomNoZ);
+
+ std::unique_ptr<Geometry> geomZExpected = r.read(wktZExpected);
+
+ geomNoZ->normalize();
+ geomZExpected->normalize();
+ std::string obtainedWKT = w.write(geomNoZ.get());
+ std::string expectedWKT = w.write(geomZExpected.get());
+ ensure_equals(obtainedWKT, expectedWKT);
+ }
+
+
+};
+
+typedef test_group<test_overlayng_elevationmodel_data> group;
+typedef group::object object;
+
+group test_overlayng_elevationmodel_group("geos::operation::overlayng::ElevationModel");
+
+//
+// Test Cases
+//
+
+// testBox
+template<>
+template<>
+void object::test<1> ()
+{
+ checkElevation(
+ "POLYGON Z ((1 6 50, 9 6 60, 9 4 50, 1 4 40, 1 6 50))",
+ {
+ 0,10, 50, 5,10, 50, 10,10, 60,
+ 0,5, 50, 5,5, 50, 10,5, 50,
+ 0,4, 40, 5,4, 50, 10,4, 50,
+ 0,0, 40, 5,0, 50, 10,0, 50
+ });
+}
+
+
+// testLine
+template<>
+template<>
+void object::test<2> ()
+{
+ checkElevation(
+ "LINESTRING Z (0 0 0, 10 10 10)",
+ {
+ -1,11, 5, 11,11, 10,
+ 0,10, 5, 5,10, 5, 10,10, 10,
+ 0,5, 5, 5,5, 5, 10,5, 5,
+ 0,0, 0, 5,0, 5, 10,0, 5,
+ -1,-1, 0, 5,-1, 5, 11,-1, 5
+ });
+}
+
+
+// testMultiLine
+template<>
+template<>
+void object::test<3> ()
+{
+ checkElevation(
+ "MULTILINESTRING Z ((0 0 0, 10 10 8), (1 2 2, 9 8 6))",
+ {
+ -1,11, 4, 11,11, 7,
+ 0,10, 4, 5,10, 4, 10,10, 7,
+ 0,5, 4, 5,5, 4, 10,5, 4,
+ 0,0, 1, 5,0, 4, 10,0, 4,
+ -1,-1, 1, 5,-1, 4, 11,-1, 4
+ });
+ }
+
+// testTwoLines
+template<>
+template<>
+void object::test<4> ()
+{
+ checkElevation( "LINESTRING Z (0 0 0, 10 10 8)",
+ "LINESTRING Z (1 2 2, 9 8 6))",
+ {
+ -1,11, 4, 11,11, 7,
+ 0,10, 4, 5,10, 4, 10,10, 7,
+ 0,5, 4, 5,5, 4, 10,5, 4,
+ 0,0, 1, 5,0, 4, 10,0, 4,
+ -1,-1, 1, 5,-1, 4, 11,-1, 4
+ });
+}
+
+/**
+ * Tests that XY geometries are scanned correctly (avoiding reading Z)
+ * and that they produce a model Z value of NaN
+ */
+// testLine2D()
+template<>
+template<>
+void object::test<5> ()
+{
+ checkElevation( "LINESTRING(0 0, 10 0)",
+ { 5, 5, DoubleNotANumber }
+ );
+}
+
+// testLineHorizontal
+template<>
+template<>
+void object::test<6> ()
+{
+ checkElevation("LINESTRING Z (0 5 0, 10 5 10)",
+ {
+ 0,10, 0, 5,10, 5, 10,10, 10,
+ 0,5, 0, 5,5, 5, 10,5, 10,
+ 0,0, 0, 5,0, 5, 10,0, 10
+ });
+}
+
+// testLineVertical
+template<>
+template<>
+void object::test<7> ()
+{
+ checkElevation("LINESTRING Z (5 0 0, 5 10 10)",
+ {
+ 0,10, 10, 5,10, 10, 10,10, 10,
+ 0,5, 5, 5,5, 5, 10,5, 5,
+ 0,0, 0, 5,0, 0, 10,0, 0
+ });
+}
+
+// tests that single point Z is used for entire grid and beyond
+// testPoint()
+template<>
+template<>
+void object::test<8> ()
+{
+ checkElevation("POINT Z (5 5 5)",
+ {
+ 0,9, 5, 5,9, 5, 9,9, 5,
+ 0,5, 5, 5,5, 5, 9,5, 5,
+ 0,0, 5, 5,0, 5, 9,0, 5
+ });
+}
+
+// tests that Z is average of input points with same location
+// testMultiPointSame
+template<>
+template<>
+void object::test<9> ()
+{
+ checkElevation("MULTIPOINT Z ((5 5 5), (5 5 9))",
+ {
+ 0,9, 7, 5,9, 7, 9,9, 7,
+ 0,5, 7, 5,5, 7, 9,5, 7,
+ 0,0, 7, 5,0, 7, 9,0, 7
+ });
+}
+
+// testPopulateZLine
+template<>
+template<>
+void object::test<10> ()
+{
+ checkElevationPopulateZ("LINESTRING Z (0 0 0, 10 10 10)",
+ "LINESTRING (1 1, 9 9)",
+ "LINESTRING (1 1 0, 9 9 10)"
+ );
+}
+
+
+// testPopulateZBox
+template<>
+template<>
+void object::test<11> ()
+{
+ checkElevationPopulateZ("LINESTRING Z (0 0 0, 10 10 10)",
+ "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+ "POLYGON Z ((1 1 0, 1 9 5, 9 9 10, 9 1 5, 1 1 0))"
+ );
+}
+
+} // namespace tut
diff --git a/tests/unit/operation/overlayng/OverlayNGZTest.cpp b/tests/unit/operation/overlayng/OverlayNGZTest.cpp
index cbe4d8c..f56ab10 100644
--- a/tests/unit/operation/overlayng/OverlayNGZTest.cpp
+++ b/tests/unit/operation/overlayng/OverlayNGZTest.cpp
@@ -4,7 +4,7 @@
//
// Ported from JTS
// modules/core/src/test/java/org/locationtech/jts/operation/overlayng/OverlayNGZTest.java
-// c1a2e9a3cd723ebfb6c79e9daf5217f562f470e0
+// f3ebebb71ea838d00d29e564cf4de667a3b9a968
//
#include <tut/tut.hpp>
@@ -12,6 +12,8 @@
// geos
#include <geos/operation/overlayng/OverlayNG.h>
+#include <geos/io/WKTWriter.h>
+#include <geos/io/WKTReader.h>
// std
#include <memory>
@@ -19,6 +21,7 @@
using namespace geos::geom;
using namespace geos::operation::overlayng;
using geos::io::WKTReader;
+using geos::io::WKTWriter;
namespace tut {
//
@@ -29,6 +32,13 @@ namespace tut {
struct test_overlayngz_data {
WKTReader r;
+ WKTWriter w;
+
+ test_overlayngz_data()
+ {
+ w.setTrim(true);
+ w.setOutputDimension(3);
+ }
void checkOverlay(int opCode, const std::string& wktA, const std::string& wktB, const std::string& wktExpected)
{
@@ -36,7 +46,11 @@ struct test_overlayngz_data {
std::unique_ptr<Geometry> b = r.read(wktB);
std::unique_ptr<Geometry> expected = r.read(wktExpected);
std::unique_ptr<Geometry> result = OverlayNG::overlay(a.get(), b.get(), opCode);
- ensure_equals_geometry_xyz(result.get(), expected.get());
+ expected->normalize();
+ result->normalize();
+ std::string obtainedWKT = w.write(result.get());
+ std::string expectedWKT = w.write(expected.get());
+ ensure_equals(obtainedWKT, expectedWKT);
}
void checkIntersection(const std::string& wktA, const std::string& wktB, const std::string& wktExpected)
@@ -65,100 +79,155 @@ group test_overlayngz_group("geos::operation::overlayng::OverlayNGZ");
// Test Cases
//
-// testLineIntersectionPointInterpolated
+// testPointXYPointDifference
template<>
template<>
void object::test<1> ()
{
- checkIntersection("LINESTRING (0 0 0, 10 10 10)", "LINESTRING (10 0 0, 0 10 10)", "POINT(5 5 5)");
+ checkDifference("MULTIPOINT ((1 1), (5 5))", "POINT Z (5 5 99)",
+ "POINT Z(1 1 99)");
}
-// testLineIntersectionPoint
+// checks that Point Z is preserved
+// testPointPolygonIntersection
template<>
template<>
void object::test<2> ()
{
- checkIntersection(
- "LINESTRING (0 0 0, 10 10 10)",
- "LINESTRING (10 0 0, 5 5 999, 0 10 10)",
- "POINT(5 5 999)"
- );
+ checkIntersection("POINT Z (5 5 99)", "POLYGON Z ((1 9 5, 9 9 9, 9 1 5, 1 1 1, 1 9 5))",
+ "POINT Z(5 5 99)");
}
-// testLineLineXYDifferenceLineInterpolate
+// testLineIntersectionPointZInterpolated
template<>
template<>
void object::test<3> ()
{
- checkDifference(
- "LINESTRING (0 0 0, 10 10 10)",
- "LINESTRING (5 5, 6 6)",
- "MULTILINESTRING ((0 0 0, 5 5 5), (6 6 6, 10 10 10))"
- );
+ checkIntersection("LINESTRING (0 0 0, 10 10 10)", "LINESTRING (10 0 0, 0 10 10)",
+ "POINT(5 5 5)");
}
-// testLineLineUnionOverlay
+// testLineIntersectionPointZValue
template<>
template<>
void object::test<4> ()
{
- checkUnion(
- "LINESTRING (0 0 0, 10 10 10)",
- "LINESTRING (5 5 990, 15 15 999)",
- "MULTILINESTRING Z((0 0 0, 5 5 990), (5 5 990, 10 10 10), (10 10 10, 15 15 999))"
- );
+ checkIntersection("LINESTRING (0 0 0, 10 10 10)", "LINESTRING (10 0 0, 5 5 999, 0 10 10)",
+ "POINT(5 5 999)");
}
-// testLinePolygonIntersectionLine
+// testLineOverlapUnion
template<>
template<>
void object::test<5> ()
{
- checkIntersection(
- "LINESTRING Z (0 0 0, 5 5 5)",
- "POLYGON Z ((1 9 5, 9 9 9, 9 1 5, 1 1 1, 1 9 5))",
- "LINESTRING Z (1 1 1, 5 5 5)"
- );
+ checkUnion("LINESTRING (0 0 0, 10 10 10)", "LINESTRING (5 5 990, 15 15 999)",
+ "MULTILINESTRING Z((0 0 0, 5 5 990), (5 5 990, 10 10 10), (10 10 10, 15 15 999))");
}
-// testLinePolygonDifferenceLine
+// testLineLineXYDifferenceLineInterpolated
template<>
template<>
void object::test<6> ()
{
- checkDifference(
- "LINESTRING Z (0 5 0, 10 5 10)",
- "POLYGON Z ((1 9 5, 9 9 9, 9 1 5, 1 1 1, 1 9 5))",
- "MULTILINESTRING Z((0 5 0, 1 5 2), (9 5 8, 10 5 10))"
- );
+ checkDifference("LINESTRING (0 0 0, 10 10 10)", "LINESTRING (5 5, 6 6)",
+ "MULTILINESTRING ((0 0 0, 5 5 5), (6 6 6, 10 10 10))");
}
-// testLinePolygonXYDifferenceLine
+// testLinePolygonIntersection
template<>
template<>
void object::test<7> ()
{
- checkDifference(
- "LINESTRING Z (0 5 0, 10 5 10)",
- "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
- "MULTILINESTRING Z((0 5 0, 1 5 1), (9 5 9, 10 5 10))"
- );
+ checkIntersection("LINESTRING Z (0 0 0, 5 5 5)", "POLYGON Z ((1 9 5, 9 9 9, 9 1 5, 1 1 1, 1 9 5))",
+ "LINESTRING Z (1 1 1, 5 5 5)");
+}
+
+// testLinePolygonDifference
+template<>
+template<>
+void object::test<8> ()
+{
+ checkDifference("LINESTRING Z (0 5 0, 10 5 10)", "POLYGON Z ((1 9 5, 9 9 9, 9 1 5, 1 1 1, 1 9 5))",
+ "MULTILINESTRING Z((0 5 0, 1 5 2), (9 5 8, 10 5 10))");
+}
+
+// testPointXYPolygonIntersection
+template<>
+template<>
+void object::test<9> ()
+{
+ checkIntersection("POINT (5 5)", "POLYGON Z ((1 9 50, 9 9 90, 9 1 50, 1 1 10, 1 9 50))",
+ "POINT Z(5 5 50)");
+}
+
+// XY Polygon gets Z value from Point
+// testPointPolygonXYUnionn
+template<>
+template<>
+void object::test<10> ()
+{
+ checkUnion("POINT Z (5 5 77)", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+ "POLYGON Z((1 1 77, 1 9 77, 9 9 77, 9 1 77, 1 1 77))");
+}
+
+// testLinePolygonXYDifference
+template<>
+template<>
+void object::test<11> ()
+{
+ checkDifference("LINESTRING Z (0 5 0, 10 5 10)", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+ "MULTILINESTRING Z((0 5 0, 1 5 1), (9 5 9, 10 5 10))");
}
-//// TODO: add Z population from model
-//// testLineXYPolygonDifferenceLine
-//template<>
-//template<>
-//void object::test<8> ()
-//{
-// checkDifference(
-// "LINESTRING (0 5, 10 5)",
-// "POLYGON Z ((1 9 5, 9 9 9, 9 1 5, 1 1 1, 1 9 5))",
-// "MULTILINESTRING Z((0 5 0, 1 5 2), (9 5 8, 10 5 10))"
-// );
-//}
+// testLineXYPolygonDifference
+template<>
+template<>
+void object::test<12> ()
+{
+ checkDifference("LINESTRING (0 5, 10 5)", "POLYGON Z ((1 9 50, 9 9 90, 9 1 50, 1 1 10, 1 9 50))",
+ "MULTILINESTRING Z((0 5 50, 1 5 30), (9 5 70, 10 5 50))");
+}
+// testPolygonXYPolygonIntersection
+template<>
+template<>
+void object::test<13> ()
+{
+ checkIntersection("POLYGON ((4 12, 2 6, 7 6, 11 4, 15 15, 4 12))", "POLYGON Z ((1 9 50, 9 9 90, 9 1 50, 1 1 10, 1 9 50))",
+ "POLYGON Z((2 6 50, 3 9 60, 9 9 90, 9 5 70, 7 6 90, 2 6 50))");
+}
+// testPolygonXYPolygonUnion
+template<>
+template<>
+void object::test<14> ()
+{
+ checkUnion("POLYGON ((0 3, 3 3, 3 0, 0 0, 0 3))", "POLYGON Z ((1 9 50, 9 9 90, 9 1 50, 1 1 10, 1 9 50))",
+ "POLYGON Z((0 0 10, 0 3 50, 1 3 20, 1 9 50, 9 9 90, 9 1 50, 3 1 20, 3 0 50, 0 0 10))");
+}
+
+// Test that operation on XY geoms produces XY (Z = NaN)
+// testPolygonXYPolygonXYIntersection
+template<>
+template<>
+void object::test<15> ()
+{
+ checkIntersection("POLYGON ((4 12, 2 6, 7 6, 11 4, 15 15, 4 12))", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+ "POLYGON ((2 6, 3 9, 9 9, 9 5, 7 6, 2 6))");
+}
+// from https://trac.osgeo.org/geos/ticket/435
+// testLineXYLineIntersection
+template<>
+template<>
+void object::test<16> ()
+{
+ checkIntersection(
+ "LINESTRING(0 0,0 10,10 10,10 0)",
+ "LINESTRING(10 10 4,10 0 5,0 0 5)",
+ "GEOMETRYCOLLECTION Z(POINT Z(0 0 5), LINESTRING Z(10 0 5, 10 10 4))"
+ );
+}
} // namespace tut
-----------------------------------------------------------------------
Summary of changes:
include/geos/geom/CoordinateSequence.h | 4 +
include/geos/operation/overlayng/ElevationModel.h | 170 +++++++++++++
include/geos/operation/overlayng/Makefile.am | 1 +
include/geos/operation/overlayng/OverlayNG.h | 5 +
src/operation/overlayng/ElevationModel.cpp | 277 +++++++++++++++++++++
src/operation/overlayng/Makefile.am | 1 +
src/operation/overlayng/OverlayNG.cpp | 66 ++++-
tests/unit/Makefile.am | 1 +
.../operation/overlayng/ElevationModelTest.cpp | 254 +++++++++++++++++++
tests/unit/operation/overlayng/OverlayNGZTest.cpp | 173 +++++++++----
10 files changed, 891 insertions(+), 61 deletions(-)
create mode 100644 include/geos/operation/overlayng/ElevationModel.h
create mode 100644 src/operation/overlayng/ElevationModel.cpp
create mode 100644 tests/unit/operation/overlayng/ElevationModelTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list