From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Thu, 5 Jul 2012 16:50:44 +0000 (+0000)
Subject: Support for larger objects in ST_Intersection(geography) (#1610) This adds a set... 
X-Git-Tag: 2.1.0beta2~818
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=6ba12da35aee3a70dd2d81fae18b8648323732a7;p=postgis

Support for larger objects in ST_Intersection(geography) (#1610) This adds a set of larger regions that are handled with a gnomic projection. It could also use an orthographic. There's still a case to be handled for objects that cross the dateline or poles.


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

diff --git a/libpgcommon/lwgeom_transform.c b/libpgcommon/lwgeom_transform.c
index 07b0e5e5d..238803644 100644
--- a/libpgcommon/lwgeom_transform.c
+++ b/libpgcommon/lwgeom_transform.c
@@ -407,6 +407,27 @@ static char* GetProj4String(int srid)
 		{
 			snprintf(proj_str, maxproj4len, "+proj=utm +zone=%d +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs", id - SRID_SOUTH_UTM_START + 1);
 		}
+		/* Gnomic zones (about 30x30, larger in higher latitudes) */
+		else if ( id >= SRID_GNOMIC_START && id <= SRID_GNOMIC_END )
+		{
+			int zone = id - SRID_GNOMIC_START;
+			int xzone = zone % 20;
+			int yzone = zone / 20;
+			double lat_0 = 30.0 * (yzone - 3) + 15.0;
+			double lon_0;
+			
+			/* The number of xzones is variable depending on yzone */
+			if  ( yzone == 2 || yzone == 3 )
+				lon_0 = 30.0 * (xzone - 6) + 15.0;
+			else if ( yzone == 1 || yzone == 4 )
+				lon_0 = 45.0 * (xzone - 4) + 22.5;
+			else if ( yzone == 0 || yzone == 5 )
+				lon_0 = 90.0 * (xzone - 2) + 45.0;
+			else
+				lwerror("Unknown yzone encountered!");
+			
+			snprintf(proj_str, maxproj4len, "+proj=gnom +ellps=WGS84 +datum=WGS84 +lat_0=%g +lon_0=%g +units=m +no_defs", lat_0, lon_0);
+		}
 		/* Lambert Azimuthal Equal Area South Pole */
 		else if ( id == SRID_SOUTH_LAMBERT )
 		{
diff --git a/libpgcommon/lwgeom_transform.h b/libpgcommon/lwgeom_transform.h
index 235151f57..11574df35 100644
--- a/libpgcommon/lwgeom_transform.h
+++ b/libpgcommon/lwgeom_transform.h
@@ -70,5 +70,12 @@ void srid_is_latlong(FunctionCallInfo fcinfo, int srid);
 /** PolarSteregraphic South, equivalent to EPSG:3031 */
 #define SRID_SOUTH_STEREO    999162
 
+/** Gnomic zones start (6 latitude bands x up to 20 latitude bands) */
+#define SRID_GNOMIC_START 999163
+
+/** Gnomic zones end (6 latitude bands x up to 20 latitude bands) */
+#define SRID_GNOMIC_END 999283
+
+
 /** @} */
 
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index ae54df620..c224a9b76 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -696,6 +696,7 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 	int type1, type2;
 	int empty1 = LW_FALSE;
 	int empty2 = LW_FALSE;
+	double ycenter, xcenter, xwidth, ywidth;
 
 	Datum d1 = PG_GETARG_DATUM(0);
 	Datum d2 = PG_GETARG_DATUM(1);
@@ -747,6 +748,11 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 	if ( empty2 )
 		gbox2 = gbox1;
 
+	ycenter = (gbox1.ymin + gbox1.ymax + gbox2.ymin + gbox2.ymax) / 4.0;
+	xcenter = (gbox1.xmin + gbox1.xmax + gbox2.xmin + gbox2.xmax) / 4.0;
+
+	xwidth = fabs(FP_MAX(gbox1.xmax, gbox2.xmax) - FP_MIN(gbox1.xmin, gbox2.xmin));
+	ywidth = fabs(FP_MAX(gbox1.ymax, gbox2.ymax) - FP_MIN(gbox1.ymin, gbox2.ymin));
 
 	/* Are these data arctic? Lambert Azimuthal Equal Area North. */
 	if ( gbox1.ymin > 65.0 && gbox2.ymin > 65.0 )
@@ -766,7 +772,7 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 	** far as a half zone past a zone boundary.
 	** Note we have no handling for the date line in here.
 	*/
-	if ( fabs(FP_MAX(gbox1.xmax, gbox2.xmax) - FP_MIN(gbox1.xmin, gbox2.xmin)) < 6.0 )
+	if ( xwidth < 6.0 )
 	{
 		/* Cheap hack to pick a zone. Average of the box center points. */
 		double dzone = (gbox1.xmin + gbox1.xmax + gbox2.xmin + gbox2.xmax) / 4.0;
@@ -785,6 +791,40 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 		}
 	}
 
+	/*
+	** Can we fit into a custom Gnomic area? (30 degrees high, variable width) 
+	** We will allow overlap into adjoining areas, but use a slightly narrower test (25) to try
+	** and lower the worst case.
+	** Again, we are hoping the dateline doesn't trip us up much
+	*/
+	if ( ywidth < 25.0 )
+	{
+		int xzone = -1;
+		int yzone = 3 + floor(ycenter / 30.0); /* (range of 0-5) */
+		
+		/* Equatorial band, 12 zones, 30 degrees wide */
+		if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 )
+		{
+			xzone = 6 + floor(xcenter / 30.0);
+		}
+		/* Temperate band, 8 zones, 45 degrees wide */
+		else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 )
+		{
+			xzone = 4 + floor(xcenter / 45.0);
+		}
+		/* Arctic band, 4 zones, 90 degrees wide */
+		else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 )
+		{
+			xzone = 2 + floor(xcenter / 90.0);
+		}
+		
+		/* Did we fit into an appropriate xzone? */
+		if ( xzone != -1 )
+		{
+			PG_RETURN_INT32(SRID_GNOMIC_START + 20 * yzone + xzone);
+		}
+	}
+
 	/*
 	** Running out of options... fall-back to Mercator
 	** and hope for the best.