Support for larger objects in ST_Intersection(geography) (#1610) This adds a set...

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 07b0e5e5d6fa28f040083d4630eecaaebcb3ca45..238803644182c34c058abf278d523c709181e8b1 100644 (file)
--- 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 235151f57d47b3e5abcdb097f8a4217c0dd32fe0..11574df35fccfd0dcca3b66a90ce19ac5aa1de25 100644 (file)
--- 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 ae54df620c8cc4a29476e627c5df764326e17dc2..c224a9b7684b6fd4cd854ea7c9c4af1f0f5d98c5 100644 (file)
--- 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.