From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Fri, 26 Feb 2016 14:25:03 +0000 (+0000)
Subject: Segmentize geography using equal length segments, from
X-Git-Tag: 2.3.0beta1~219
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=119c3ba64b91b7f05881c96c01bdd6162895b3e8;p=postgis

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 33a7f14c9..854346524 100644
--- 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 f9a769aff..bca28b094 100644
--- 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 dc9f2be96..577e53b1e 100644
--- 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 8a61b1371..5c60626de 100644
--- 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