[postgis-tickets] r14706 - Segmentize geography using equal length segments, from

[Paul Ramsey]

Author: pramsey
Date: 2016-02-26 06:25:03 -0800 (Fri, 26 Feb 2016)
New Revision: 14706

Modified:
   trunk/liblwgeom/lwgeodetic.c
   trunk/regress/geography.sql
   trunk/regress/geography_expected
   trunk/regress/tickets_expected
Log:
Segmentize geography using equal length segments, from
https://github.com/postgis/postgis/pull/90
Hugo Mercier, @mhugo


Modified: trunk/liblwgeom/lwgeodetic.c
===================================================================
--- trunk/liblwgeom/lwgeodetic.c	2016-02-26 11:37:39 UTC (rev 14705)
+++ trunk/liblwgeom/lwgeodetic.c	2016-02-26 14:25:03 UTC (rev 14706)
@@ -1550,7 +1550,6 @@
 	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 @@
 		{
 			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 )

Modified: trunk/regress/geography.sql
===================================================================
--- trunk/regress/geography.sql	2016-02-26 11:37:39 UTC (rev 14705)
+++ trunk/regress/geography.sql	2016-02-26 14:25:03 UTC (rev 14706)
@@ -107,6 +107,13 @@
 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
+    

Modified: trunk/regress/geography_expected
===================================================================
--- trunk/regress/geography_expected	2016-02-26 11:37:39 UTC (rev 14705)
+++ trunk/regress/geography_expected	2016-02-26 14:25:03 UTC (rev 14706)
@@ -27,3 +27,4 @@
 #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

Modified: trunk/regress/tickets_expected
===================================================================
--- trunk/regress/tickets_expected	2016-02-26 11:37:39 UTC (rev 14705)
+++ trunk/regress/tickets_expected	2016-02-26 14:25:03 UTC (rev 14706)
@@ -236,7 +236,7 @@
 #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