Segmentize geography using equal length segments, from

https://github.com/postgis/postgis/pull/90
Hugo Mercier, @mhugo

git-svn-id: http://svn.osgeo.org/postgis/trunk@14706 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 33a7f14c9d9f855c8490f155064f4002a507e335..85434652409a81a30f1a9bbdc337cb6cbf112d9c 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1550,7 +1550,6 @@ ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length)
        int hasm = ptarray_has_m(pa_in);
        int pa_in_offset = 0; /* input point offset */
        POINT4D p1, p2, p;
-       POINT3D q1, q2, q, qn;
        GEOGRAPHIC_POINT g1, g2, g;
        double d;
        
@@ -1592,32 +1591,21 @@ ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length)
                {
                        int nsegs = 1 + d / max_seg_length;
                        int i;
-                       double dx, dy, dz, dzz = 0, dmm = 0;
-                       
-                       geog2cart(&g1, &q1);
-                       geog2cart(&g2, &q2);
-                       
-                       dx = (q2.x - q1.x) / nsegs;
-                       dy = (q2.y - q1.y) / nsegs;
-                       dz = (q2.z - q1.z) / nsegs;
-                       
+                       double dx, dy, dzz = 0, dmm = 0;
+                       g = g1;
+                       dx = (g2.lat - g1.lat) / nsegs;
+                       dy = (g2.lon - g1.lon) / nsegs;
+
                        /* The independent Z/M values on the ptarray */
                        if ( hasz ) dzz = (p2.z - p1.z) / nsegs;
                        if ( hasm ) dmm = (p2.m - p1.m) / nsegs;
                        
-                       q = q1;
                        p = p1;
-                       
                        for ( i = 0; i < nsegs - 1; i++ )
                        {
                                /* Move one increment forwards */
-                               q.x += dx; q.y += dy; q.z += dz;
-                               qn = q;
-                               normalize(&qn);
-                               
-                               /* Back to spherical coordinates */
-                               cart2geog(&qn, &g);
-                               /* Back to lon/lat */
+                               g.lat += dx;
+                               g.lon += dy;
                                p.x = rad2deg(g.lon);
                                p.y = rad2deg(g.lat);
                                if ( hasz )

diff --git a/regress/geography.sql b/regress/geography.sql
index f9a769afffb3eacd0cda4d0d588182d914469178..bca28b0946a5efa916deac3a9c52f339e2f4570c 100644 (file)
--- a/regress/geography.sql
+++ b/regress/geography.sql
@@ -107,6 +107,13 @@ CROSS JOIN (VALUES (1609),(1600),(1068)) AS t (radius)
 ORDER BY distance_t;
 DROP TABLE dwithgeogbug;
 
+-- Test segmentize on geography
+-- Check that the maximum sphere distance over all the segments is lesser than the max distance used for segmentize
+WITH
+seg as (select ST_Segmentize('LINESTRING(0 0 10,0 90 20)'::geography, 50000)::geometry as geom),
+dumped as (SELECT (st_dumppoints(geom)).path[1] as id, (st_dumppoints(geom)).geom from seg)
+SELECT 'segmentize_geography', max(st_distance(d1.geom::geography, d2.geom::geography, false))::int FROM dumped as d1, dumped as d2 where d2.id = d1.id + 1;
+
 -- Clean up spatial_ref_sys
 DELETE FROM spatial_ref_sys WHERE srid IN (4269,4326);
-    
\ No newline at end of file
+    

diff --git a/regress/geography_expected b/regress/geography_expected
index dc9f2be96c4a679240567dac0e8d0040d71fe5d5..577e53b1e1ceb5e9841c8ab00f906e558814be88 100644 (file)
--- a/regress/geography_expected
+++ b/regress/geography_expected
@@ -27,3 +27,4 @@ geog_precision_pazafir|0|0
 #2422|1|1609|t|t|1400.230|1396.816|1400.230|1400.230
 #2422|1|1600|t|t|1400.230|1396.816|1400.230|1400.230
 #2422|1|1068|f|f|1400.230|1396.816|1400.230|1400.230
+segmentize_geography|49789

diff --git a/regress/tickets_expected b/regress/tickets_expected
index 8a61b1371a2119355c9553fc7f8ffe3a3b867349..5c60626dea8e418cc3f938af649a0f6028abbe34 100644 (file)
--- a/regress/tickets_expected
+++ b/regress/tickets_expected
@@ -236,7 +236,7 @@ ERROR:  invalid GML representation
 #2110.1|f
 #2110.2|t
 #2110.3|t
-#2145|6792004
+#2145|6792007
 #2232|M 0 0 l 0 0 1 0 0 0 1 0 0 0
 #2307|MULTIPOLYGON(((-41.1932 -7.3257,-41.1616 -7.3257,-41.1569 -7.3257,-41.1569 -7.3483,-41.1932 -7.3483,-41.1932 -7.3257),(-41.1616 -7.3257,-41.1879 -7.3257,-41.1879 -7.3425,-41.1616 -7.3425,-41.1616 -7.3257)))
 #2409|GeometryCollection[B] with 2 elements