Resolving triangulation issues

R/sparkR-mosaic/sparkrMosaic/tests/testthat/testRasterFunctions.R

@@ -8,7 +8,7 @@ generate_singleband_raster_df <- function() {
 
 test_that("mosaic can read single-band GeoTiff", {
   sdf <- generate_singleband_raster_df()
-  
+
   row <- first(sdf)
   expect_equal(row$length, 1067862L)
   expect_equal(row$x_size, 2400)
@@ -158,14 +158,16 @@ sdf <- createDataFrame(
 
 sdf <- agg(groupBy(sdf), masspoints = collect_list(column("wkt")))
 sdf <- withColumn(sdf, "breaklines", expr("array('LINESTRING EMPTY')"))
+sdf <- withColumn(sdf, "splitPointFinder", lit("NONENCROACHING"))
 sdf <- withColumn(sdf, "origin", st_geomfromwkt(lit("POINT (0.6 1.8)")))
 sdf <- withColumn(sdf, "xWidth", lit(12L))
 sdf <- withColumn(sdf, "yWidth", lit(6L))
 sdf <- withColumn(sdf, "xSize", lit(0.1))
 sdf <- withColumn(sdf, "ySize", lit(0.1))
+sdf <- withColumn(sdf, "noData", lit(-9999.0))
 sdf <- withColumn(sdf, "tile", rst_dtmfromgeoms(
-column("masspoints"), column("breaklines"), lit(0.0), lit(0.01),
-column("origin"), column("xWidth"), column("yWidth"), column("xSize"), column("ySize"))
-)
+column("masspoints"), column("breaklines"), lit(0.0), lit(0.01), column("splitPointFinder"),
+column("origin"), column("xWidth"), column("yWidth"), column("xSize"), column("ySize"), column("noData")
+))
 expect_equal(SparkR::count(sdf), 1)
 })

R/sparkR-mosaic/sparkrMosaic/tests/testthat/testVectorFunctions.R

@@ -109,23 +109,24 @@ sdf <- createDataFrame(
 
 sdf <- agg(groupBy(sdf), masspoints = collect_list(column("wkt")))
 sdf <- withColumn(sdf, "breaklines", expr("array('LINESTRING EMPTY')"))
-triangulation_sdf <- withColumn(sdf, "triangles", st_triangulate(column("masspoints"), column("breaklines"), lit(0.0), lit(0.01)))
+triangulation_sdf <- withColumn(sdf, "triangles", st_triangulate(column("masspoints"), column("breaklines"), lit(0.0), lit(0.01), lit("NONENCROACHING")))
 cache(triangulation_sdf)
 expect_equal(SparkR::count(triangulation_sdf), 2)
 expected <- c("POLYGON Z((0 2 2, 2 1 0, 1 3 3, 0 2 2))", "POLYGON Z((1 3 3, 2 1 0, 3 2 1, 1 3 3))")
 expect_contains(expected, first(triangulation_sdf)$triangles)
 
 interpolation_sdf <- sdf
-interpolation_sdf <- withColumn(interpolation_sdf, "origin", st_geomfromwkt(lit("POINT (0.6 1.8)")))
+interpolation_sdf <- withColumn(interpolation_sdf, "origin", st_geomfromwkt(lit("POINT (0.55 1.75)")))
 interpolation_sdf <- withColumn(interpolation_sdf, "xWidth", lit(12L))
 interpolation_sdf <- withColumn(interpolation_sdf, "yWidth", lit(6L))
 interpolation_sdf <- withColumn(interpolation_sdf, "xSize", lit(0.1))
 interpolation_sdf <- withColumn(interpolation_sdf, "ySize", lit(0.1))
 interpolation_sdf <- withColumn(interpolation_sdf, "interpolated", st_interpolateelevation(
   column("masspoints"),
   column("breaklines"),
-  lit(0.0),
   lit(0.01),
+  lit(0.01),
+  lit("NONENCROACHING"),
   column("origin"),
   column("xWidth"),
   column("yWidth"),
@@ -134,5 +135,5 @@ interpolation_sdf <- withColumn(interpolation_sdf, "interpolated", st_interpolat
 ))
 cache(interpolation_sdf)
 expect_equal(SparkR::count(interpolation_sdf), 6 * 12)
-expect_contains(collect(interpolation_sdf)$interpolated, "POINT Z(0.6 2 1.8)")
+expect_contains(collect(interpolation_sdf)$interpolated, "POINT Z(1.6 1.8 1.2)")
 })

R/sparklyr-mosaic/sparklyrMosaic/tests/testthat/testRasterFunctions.R

@@ -209,11 +209,13 @@ test_that ("a terrain model can be produced from point geometries", {
         breaklines,
         as.double(0.0),
         as.double(0.01),
+        "NONENCROACHING",
         origin,
         xWidth,
         yWidth,
         xSize,
-        ySize
+        ySize,
+        as.double(-9999.0)
       )
     )
   expect_equal(sdf_nrow(sdf), 1)

R/sparklyr-mosaic/sparklyrMosaic/tests/testthat/testVectorFunctions.R

@@ -134,7 +134,7 @@ test_that ("triangulation and interpolation functions behave as intended", {
     mutate(breaklines = array("LINESTRING EMPTY"))
 
   triangulation_sdf <- sdf %>%
-    mutate(triangles = st_triangulate(masspoints, breaklines, as.double(0.00), as.double(0.01)))
+    mutate(triangles = st_triangulate(masspoints, breaklines, as.double(0.00), as.double(0.01), "NONENCROACHING"))
 
   expect_equal(sdf_nrow(triangulation_sdf), 2)
 
@@ -144,16 +144,17 @@ test_that ("triangulation and interpolation functions behave as intended", {
 
   interpolation_sdf <- sdf %>%
     mutate(
-      origin = st_geomfromwkt("POINT (0.6 1.8)"),
+      origin = st_geomfromwkt("POINT (0.55 1.75)"),
       xWidth = 12L,
       yWidth = 6L,
       xSize = as.double(0.1),
       ySize = as.double(0.1),
       interpolated = st_interpolateelevation(
         masspoints,
         breaklines,
-        as.double(0.0),
         as.double(0.01),
+        as.double(0.01),
+        "NONENCROACHING",
         origin,
         xWidth,
         yWidth,
@@ -163,5 +164,5 @@ test_that ("triangulation and interpolation functions behave as intended", {
     )
   expect_equal(sdf_nrow(interpolation_sdf), 6 * 12)
   expect_contains(sdf_collect(interpolation_sdf)$interpolated,
-                  "POINT Z(0.6 2 1.8)")
+                  "POINT Z(1.6 1.8 1.2)")
 })

docs/docs-requirements.txt

@@ -1,4 +1,4 @@
-setuptools==68.1.2
+setuptools>=70.0.0
 Sphinx==6.1.3
 sphinx-material==0.0.36
 nbsphinx>=0.8.8

docs/source/api/raster-functions.rst

@@ -491,7 +491,7 @@ rst_derivedband
 rst_dtmfromgeoms
 ****************
 
-.. function:: rst_dtmfromgeoms(pointsArray, linesArray, mergeTolerance, snapTolerance, origin, xWidth, yWidth, xSize, ySize)
+.. function:: rst_dtmfromgeoms(pointsArray, linesArray, mergeTolerance, snapTolerance, splitPointFinder, origin, xWidth, yWidth, xSize, ySize, noData)
 
     Generate a raster with interpolated elevations across a grid of points described by:
 
@@ -518,6 +518,13 @@ rst_dtmfromgeoms
       Setting this value to zero may result in the output triangle vertices being assigned a null Z value.
     Both tolerance parameters are expressed in the same units as the projection of the input point geometries.
 
+    Additionally, you have control over the algorithm used to find split points on the constraint lines. The recommended
+    default option here is the "NONENCROACHING" algorithm. You can also use the "MIDPOINT" algorithm if you find the
+    constraint fitting process fails to converge. For full details of these options see the JTS reference
+    `here <https://locationtech.github.io/jts/javadoc/org/locationtech/jts/triangulate/ConstraintSplitPointFinder.html>`__.
+
+    The :code:`noData` value of the output raster can be set using the :code:`noData` parameter.
+
     This is a generator expression and the resulting DataFrame will contain one row per point of the grid.
 
     :param pointsArray: Array of geometries respresenting the points to be triangulated
@@ -528,6 +535,8 @@ rst_dtmfromgeoms
     :type mergeTolerance: Column (DoubleType)
     :param snapTolerance: A snapping tolerance used to relate created points to their corresponding lines for elevation interpolation.
     :type snapTolerance: Column (DoubleType)
+    :param splitPointFinder: Algorithm used for finding split points on constraint lines. Options are "NONENCROACHING" and "MIDPOINT".
+    :type splitPointFinder: Column (StringType)
     :param origin: A point geometry describing the bottom-left corner of the grid.
     :type origin: Column (Geometry)
     :param xWidth: The number of points in the grid in x direction.
@@ -538,6 +547,8 @@ rst_dtmfromgeoms
     :type xSize: Column (DoubleType)
     :param ySize: The spacing between each point on the grid's y-axis.
     :type ySize: Column (DoubleType)
+    :param noData: The no-data value of the output raster.
+    :type noData: Column (DoubleType)
     :rtype: Column (RasterTileType)
 
     :example:
@@ -567,7 +578,8 @@ rst_dtmfromgeoms
     df.select(
         rst_dtmfromgeoms(
             "masspoints", "breaklines", lit(0.0), lit(0.01),
-            "origin", "xWidth", "yWidth", "xSize", "ySize"
+            "origin", "xWidth", "yWidth", "xSize", "ySize",
+            split_point_finder="NONENCROACHING", no_data_value=-9999.0
         )
     ).show(truncate=False)
     +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
@@ -594,8 +606,10 @@ rst_dtmfromgeoms
 
     df.select(
       rst_dtmfromgeoms(
-        $"masspoints", $"breaklines", lit(0.0), lit(0.01),
-        $"origin", $"xWidth", $"yWidth", $"xSize", $"ySize"
+        $"masspoints", $"breaklines",
+        lit(0.0), lit(0.01), lit("NONENCROACHING")
+        $"origin", $"xWidth", $"yWidth",
+        $"xSize", $"ySize", lit(-9999.0)
       )
     ).show(1, false)
     +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
@@ -615,8 +629,9 @@ rst_dtmfromgeoms
           "POINT Z (0 2 2)"
         ),
         ARRAY("LINESTRING EMPTY"),
-        DOUBLE(0.0), DOUBLE(0.01),
-        "POINT (0.6 1.8)", 12, 6, DOUBLE(0.1), DOUBLE(0.1)
+        DOUBLE(0.0), DOUBLE(0.01), "NONENCROACHING",
+        "POINT (0.6 1.8)", 12, 6,
+        DOUBLE(0.1), DOUBLE(0.1), DOUBLE(-9999.0)
       ) AS tile
     +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
     |rst_dtmfromgeoms(masspoints, breaklines, 0.0, 0.01, origin, xWidth, yWidth, xSize, ySize)                                                                                                                                                                                                                                                                                                                                                              |
@@ -638,8 +653,9 @@ rst_dtmfromgeoms
     sdf <- withColumn(sdf, "breaklines", expr("array('LINESTRING EMPTY')"))
     sdf <- select(sdf, rst_dtmfromgeoms(
       column("masspoints"), column("breaklines"),
-      lit(0.0), lit(0.01),
-      lit("POINT (0.6 1.8)"), lit(12L), lit(6L), lit(0.1), lit(0.1)
+      lit(0.0), lit(0.01), lit("NONENCROACHING"),
+      lit("POINT (0.6 1.8)"), lit(12L), lit(6L),
+      lit(0.1), lit(0.1), lit(-9999.0)
       )
     )
     showDF(sdf, n=1, truncate=F)

docs/source/api/spatial-functions.rst

@@ -881,7 +881,7 @@ st_haversine
 st_interpolateelevation
 ***********************
 
-.. function:: st_interpolateelevation(pointsArray, linesArray, mergeTolerance, snapTolerance, origin, xWidth, yWidth, xSize, ySize)
+.. function:: st_interpolateelevation(pointsArray, linesArray, mergeTolerance, snapTolerance, splitPointFinder, origin, xWidth, yWidth, xSize, ySize)
 
     Compute interpolated elevations across a grid of points described by:
 
@@ -908,6 +908,11 @@ st_interpolateelevation
       Setting this value to zero may result in the output triangle vertices being assigned a null Z value.
     Both tolerance parameters are expressed in the same units as the projection of the input point geometries.
 
+    Additionally, you have control over the algorithm used to find split points on the constraint lines. The recommended
+    default option here is the "NONENCROACHING" algorithm. You can also use the "MIDPOINT" algorithm if you find the
+    constraint fitting process fails to converge. For full details of these options see the JTS reference
+    `here <https://locationtech.github.io/jts/javadoc/org/locationtech/jts/triangulate/ConstraintSplitPointFinder.html>`__.
+
     This is a generator expression and the resulting DataFrame will contain one row per point of the grid.
 
     :param pointsArray: Array of geometries respresenting the points to be triangulated
@@ -919,6 +924,8 @@ st_interpolateelevation
     :param snapTolerance: A snapping tolerance used to relate created points to their corresponding lines for elevation interpolation.
     :type snapTolerance: Column (DoubleType)
     :param origin: A point geometry describing the bottom-left corner of the grid.
+    :param splitPointFinder: Algorithm used for finding split points on constraint lines. Options are "NONENCROACHING" and "MIDPOINT".
+    :type splitPointFinder: Column (StringType)
     :type origin: Column (Geometry)
     :param xWidth: The number of points in the grid in x direction.
     :type xWidth: Column (IntegerType)
@@ -957,7 +964,8 @@ st_interpolateelevation
     df.select(
         st_interpolateelevation(
             "masspoints", "breaklines", lit(0.0), lit(0.01),
-            "origin", "xWidth", "yWidth", "xSize", "ySize"
+            "origin", "xWidth", "yWidth", "xSize", "ySize",
+            split_point_finder="NONENCROACHING"
         )
     ).show(4, truncate=False)
     +--------------------------------------------------+
@@ -987,7 +995,8 @@ st_interpolateelevation
 
     df.select(
       st_interpolateelevation(
-        $"masspoints", $"breaklines", lit(0.0), lit(0.01),
+        $"masspoints", $"breaklines",
+        lit(0.0), lit(0.01), lit("NONENCROACHING"),
         $"origin", $"xWidth", $"yWidth", $"xSize", $"ySize"
       )
     ).show(4, false)
@@ -1011,7 +1020,7 @@ st_interpolateelevation
           "POINT Z (0 2 2)"
         ),
         ARRAY("LINESTRING EMPTY"),
-        DOUBLE(0.0), DOUBLE(0.01),
+        DOUBLE(0.0), DOUBLE(0.01), "NONENCROACHING",
         "POINT (0.6 1.8)", 12, 6, DOUBLE(0.1), DOUBLE(0.1)
       )
     +--------------------------------------------------+
@@ -1037,7 +1046,7 @@ st_interpolateelevation
     sdf <- withColumn(sdf, "breaklines", expr("array('LINESTRING EMPTY')"))
     sdf <- select(sdf, st_interpolateelevation(
       column("masspoints"), column("breaklines"),
-      lit(0.0), lit(0.01),
+      lit(0.0), lit(0.01), lit("NONENCROACHING"),
       lit("POINT (0.6 1.8)"), lit(12L), lit(6L), lit(0.1), lit(0.1)
       )
     )
@@ -1757,7 +1766,7 @@ st_transform
 st_triangulate
 **************
 
-.. function:: st_triangulate(pointsArray, linesArray, mergeTolerance, snapTolerance)
+.. function:: st_triangulate(pointsArray, linesArray, mergeTolerance, snapTolerance, splitPointFinder)
 
     Performs a conforming Delaunay triangulation using the points in :code:`pointsArray` including :code:`linesArray` as constraint / break lines.
 
@@ -1773,6 +1782,11 @@ st_triangulate
       Setting this value to zero may result in the output triangle vertices being assigned a null Z value.
     Both tolerance parameters are expressed in the same units as the projection of the input point geometries.
 
+    Additionally, you have control over the algorithm used to find split points on the constraint lines. The recommended
+    default option here is the "NONENCROACHING" algorithm. You can also use the "MIDPOINT" algorithm if you find the
+    constraint fitting process fails to converge. For full details of these options see the JTS reference
+    `here <https://locationtech.github.io/jts/javadoc/org/locationtech/jts/triangulate/ConstraintSplitPointFinder.html>`__.
+
     This is a generator expression and the resulting DataFrame will contain one row per triangle returned by the algorithm.
 
     :param pointsArray: Array of geometries respresenting the points to be triangulated
@@ -1783,6 +1797,8 @@ st_triangulate
     :type mergeTolerance: Column (DoubleType)
     :param snapTolerance: A snapping tolerance used to relate created points to their corresponding lines for elevation interpolation.
     :type snapTolerance: Column (DoubleType)
+    :param splitPointFinder: Algorithm used for finding split points on constraint lines. Options are "NONENCROACHING" and "MIDPOINT".
+    :type splitPointFinder: Column (StringType)
     :rtype: Column (Geometry)
 
     :example:
@@ -1803,7 +1819,7 @@ st_triangulate
       .groupBy()
       .agg(collect_list("wkt").alias("masspoints"))
       .withColumn("breaklines", array(lit("LINESTRING EMPTY")))
-      .withColumn("triangles", st_triangulate("masspoints", "breaklines", lit(0.0), lit(0.01))
+      .withColumn("triangles", st_triangulate("masspoints", "breaklines", lit(0.0), lit(0.01), "NONENCROACHING"))
     )
     df.show(2, False)
     +---------------------------------------+
@@ -1826,7 +1842,7 @@ st_triangulate
     .withColumn("triangles",
       st_triangulate(
         $"masspoints", $"breaklines",
-        lit(0.0), lit(0.01)
+        lit(0.0), lit(0.01), lit("NONENCROACHING")
       )
     )
 
@@ -1849,7 +1865,8 @@ st_triangulate
           "POINT Z (0 2 2)"
         ),
         ARRAY("LINESTRING EMPTY"),
-        DOUBLE(0.0), DOUBLE(0.01)
+        DOUBLE(0.0), DOUBLE(0.01),
+        "NONENCROACHING"
       )
     +---------------------------------------+
     |triangles                              |
@@ -1872,7 +1889,7 @@ st_triangulate
     sdf <- withColumn(sdf, "breaklines", expr("array('LINESTRING EMPTY')"))
     result <- select(sdf, st_triangulate(
       column("masspoints"), column("breaklines"),
-      lit(0.0), lit(0.01))
+      lit(0.0), lit(0.01), lit("NONENCROACHING")
       )
     showDF(result, truncate=F)
     +---------------------------------------+

pom.xml

@@ -98,12 +98,12 @@
         <dependency>
             <groupId>org.locationtech.jts</groupId>
             <artifactId>jts-core</artifactId>
-            <version>1.19.0</version>
+            <version>1.20.0</version>
         </dependency>
         <dependency>
             <groupId>org.locationtech.jts.io</groupId>
             <artifactId>jts-io-common</artifactId>
-            <version>1.19.0</version>
+            <version>1.20.0</version>
         </dependency>
         <dependency>
             <groupId>org.locationtech.proj4j</groupId>