From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Fri, 10 Aug 2012 16:00:31 +0000 (+0000)
Subject: Switch from gnomic to LAEA for the projection for the custom zones. Less perfect... 
X-Git-Tag: 2.1.0beta2~699
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=f9f7133a1b4599ff707b86a3fd5a8533db08e008;p=postgis

Switch from gnomic to LAEA for the projection for the custom zones. Less perfect intersections, but measure metric fidelity for distances, areas, etc. (#1610)


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

diff --git a/libpgcommon/lwgeom_transform.c b/libpgcommon/lwgeom_transform.c
index 238803644..dd26c4528 100644
--- a/libpgcommon/lwgeom_transform.c
+++ b/libpgcommon/lwgeom_transform.c
@@ -407,14 +407,18 @@ 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 )
+		/* Lambert zones (about 30x30, larger in higher latitudes) */
+		/* There are three latitude zones, divided at -90,-60,-30,0,30,60,90. */
+		/* In yzones 2,3 (equator) zones, the longitudinal zones are divided every 30 degrees (12 of them) */
+		/* In yzones 1,4 (temperate) zones, the longitudinal zones are every 45 degrees (8 of them) */
+		/* In yzones 0,5 (polar) zones, the longitudinal zones are ever 90 degrees (4 of them) */
+		else if ( id >= SRID_LAEA_START && id <= SRID_LAEA_END )
 		{
-			int zone = id - SRID_GNOMIC_START;
+			int zone = id - SRID_LAEA_START;
 			int xzone = zone % 20;
 			int yzone = zone / 20;
 			double lat_0 = 30.0 * (yzone - 3) + 15.0;
-			double lon_0;
+			double lon_0 = 0.0;
 			
 			/* The number of xzones is variable depending on yzone */
 			if  ( yzone == 2 || yzone == 3 )
@@ -426,7 +430,7 @@ static char* GetProj4String(int srid)
 			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);
+			snprintf(proj_str, maxproj4len, "+proj=laea +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 11574df35..ee79179cc 100644
--- a/libpgcommon/lwgeom_transform.h
+++ b/libpgcommon/lwgeom_transform.h
@@ -52,7 +52,7 @@ void srid_is_latlong(FunctionCallInfo fcinfo, int srid);
 /**  End of UTM North zone, equivalent to EPSG:32660 */
 #define SRID_NORTH_UTM_END   999060
 
-/** Lambert Azimuthal Equal Area North Pole, equivalent to EPSG:3574 */
+/** Lambert Azimuthal Equal Area (LAEA) North Pole, equivalent to EPSG:3574 */
 #define SRID_NORTH_LAMBERT   999061
 
 /** PolarSteregraphic North, equivalent to EPSG:3995 */
@@ -64,17 +64,17 @@ void srid_is_latlong(FunctionCallInfo fcinfo, int srid);
 /**  Start of UTM South zone, equivalent to EPSG:32760 */
 #define SRID_SOUTH_UTM_END   999160
 
-/** Lambert Azimuthal Equal Area South Pole, equivalent to EPSG:3409 */
+/** Lambert Azimuthal Equal Area (LAEA) South Pole, equivalent to EPSG:3409 */
 #define SRID_SOUTH_LAMBERT   999161
 
 /** 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
+/** LAEA zones start (6 latitude bands x up to 20 longitude bands) */
+#define SRID_LAEA_START 999163
 
-/** Gnomic zones end (6 latitude bands x up to 20 latitude bands) */
-#define SRID_GNOMIC_END 999283
+/** LAEA zones end (6 latitude bands x up to 20 longitude bands) */
+#define SRID_LAEA_END 999283
 
 
 /** @} */
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 57a6d8e3e..76ac3bdd7 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -774,9 +774,9 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 	}
 
 	/*
-	** Can we fit into a custom Gnomic area? (30 degrees high, variable width) 
+	** Can we fit into a custom LAEA 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.
+	** and minimize the worst case.
 	** Again, we are hoping the dateline doesn't trip us up much
 	*/
 	if ( ywidth < 25.0 )
@@ -803,7 +803,7 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
 		/* Did we fit into an appropriate xzone? */
 		if ( xzone != -1 )
 		{
-			PG_RETURN_INT32(SRID_GNOMIC_START + 20 * yzone + xzone);
+			PG_RETURN_INT32(SRID_LAEA_START + 20 * yzone + xzone);
 		}
 	}