PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ geography_bestsrid()

Datum geography_bestsrid ( PG_FUNCTION_ARGS  )

Definition at line 782 of file geography_measurement.c.

References GBOX::flags, GSERIALIZED::flags, gbox_angular_height(), gbox_angular_width(), gbox_centroid(), gbox_to_string(), gbox_union(), geography_project(), gserialized_get_gbox_p(), gserialized_is_empty(), LW_FAILURE, LW_FALSE, PG_FUNCTION_INFO_V1(), POINT2D::x, and POINT2D::y.

Referenced by geography_covers().

783 {
784  GBOX gbox, gbox1, gbox2;
785  GSERIALIZED *g1 = NULL;
786  GSERIALIZED *g2 = NULL;
787  int empty1 = LW_FALSE;
788  int empty2 = LW_FALSE;
789  double xwidth, ywidth;
790  POINT2D center;
791 
792  Datum d1 = PG_GETARG_DATUM(0);
793  Datum d2 = PG_GETARG_DATUM(1);
794 
795  /* Get our geometry objects loaded into memory. */
796  g1 = (GSERIALIZED*)PG_DETOAST_DATUM(d1);
797  /* Synchronize our box types */
798  gbox1.flags = g1->flags;
799  /* Calculate if the geometry is empty. */
800  empty1 = gserialized_is_empty(g1);
801  /* Calculate a geocentric bounds for the objects */
802  if ( ! empty1 && gserialized_get_gbox_p(g1, &gbox1) == LW_FAILURE )
803  elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g1, &gbox1)");
804 
805  POSTGIS_DEBUGF(4, "calculated gbox = %s", gbox_to_string(&gbox1));
806 
807  /* If we have a unique second argument, fill in all the necessary variables. */
808  if ( d1 != d2 )
809  {
810  g2 = (GSERIALIZED*)PG_DETOAST_DATUM(d2);
811  gbox2.flags = g2->flags;
812  empty2 = gserialized_is_empty(g2);
813  if ( ! empty2 && gserialized_get_gbox_p(g2, &gbox2) == LW_FAILURE )
814  elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g2, &gbox2)");
815  }
816  /*
817  ** If no unique second argument, copying the box for the first
818  ** argument will give us the right answer for all subsequent tests.
819  */
820  else
821  {
822  gbox = gbox2 = gbox1;
823  }
824 
825  /* Both empty? We don't have an answer. */
826  if ( empty1 && empty2 )
827  PG_RETURN_NULL();
828 
829  /* One empty? We can use the other argument values as infill. Otherwise merge the boxen */
830  if ( empty1 )
831  gbox = gbox2;
832  else if ( empty2 )
833  gbox = gbox1;
834  else
835  gbox_union(&gbox1, &gbox2, &gbox);
836 
837  gbox_centroid(&gbox, &center);
838 
839  /* Width and height in degrees */
840  xwidth = 180.0 * gbox_angular_width(&gbox) / M_PI;
841  ywidth = 180.0 * gbox_angular_height(&gbox) / M_PI;
842 
843  POSTGIS_DEBUGF(2, "xwidth %g", xwidth);
844  POSTGIS_DEBUGF(2, "ywidth %g", ywidth);
845  POSTGIS_DEBUGF(2, "center POINT(%g %g)", center.x, center.y);
846 
847  /* Are these data arctic? Lambert Azimuthal Equal Area North. */
848  if ( center.y > 70.0 && ywidth < 45.0 )
849  {
850  PG_RETURN_INT32(SRID_NORTH_LAMBERT);
851  }
852 
853  /* Are these data antarctic? Lambert Azimuthal Equal Area South. */
854  if ( center.y < -70.0 && ywidth < 45.0 )
855  {
856  PG_RETURN_INT32(SRID_SOUTH_LAMBERT);
857  }
858 
859  /*
860  ** Can we fit these data into one UTM zone?
861  ** We will assume we can push things as
862  ** far as a half zone past a zone boundary.
863  ** Note we have no handling for the date line in here.
864  */
865  if ( xwidth < 6.0 )
866  {
867  int zone = floor((center.x + 180.0) / 6.0);
868 
869  if ( zone > 59 ) zone = 59;
870 
871  /* Are these data below the equator? UTM South. */
872  if ( center.y < 0.0 )
873  {
874  PG_RETURN_INT32( SRID_SOUTH_UTM_START + zone );
875  }
876  /* Are these data above the equator? UTM North. */
877  else
878  {
879  PG_RETURN_INT32( SRID_NORTH_UTM_START + zone );
880  }
881  }
882 
883  /*
884  ** Can we fit into a custom LAEA area? (30 degrees high, variable width)
885  ** We will allow overlap into adjoining areas, but use a slightly narrower test (25) to try
886  ** and minimize the worst case.
887  ** Again, we are hoping the dateline doesn't trip us up much
888  */
889  if ( ywidth < 25.0 )
890  {
891  int xzone = -1;
892  int yzone = 3 + floor(center.y / 30.0); /* (range of 0-5) */
893 
894  /* Equatorial band, 12 zones, 30 degrees wide */
895  if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 )
896  {
897  xzone = 6 + floor(center.x / 30.0);
898  }
899  /* Temperate band, 8 zones, 45 degrees wide */
900  else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 )
901  {
902  xzone = 4 + floor(center.x / 45.0);
903  }
904  /* Arctic band, 4 zones, 90 degrees wide */
905  else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 )
906  {
907  xzone = 2 + floor(center.x / 90.0);
908  }
909 
910  /* Did we fit into an appropriate xzone? */
911  if ( xzone != -1 )
912  {
913  PG_RETURN_INT32(SRID_LAEA_START + 20 * yzone + xzone);
914  }
915  }
916 
917  /*
918  ** Running out of options... fall-back to Mercator
919  ** and hope for the best.
920  */
921  PG_RETURN_INT32(SRID_WORLD_MERCATOR);
922 
923 }
int gserialized_get_gbox_p(const GSERIALIZED *g, GBOX *box)
Read the bounding box off a serialization and calculate one if it is not already there.
Definition: g_serialized.c:639
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: g_box.c:399
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:258
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: g_serialized.c:178
#define LW_FAILURE
Definition: liblwgeom.h:78
double x
Definition: liblwgeom.h:330
#define LW_FALSE
Definition: liblwgeom.h:76
double y
Definition: liblwgeom.h:330
int gbox_union(const GBOX *g1, const GBOX *g2, GBOX *gout)
Update the output GBOX to be large enough to include both inputs.
Definition: g_box.c:142
uint8_t flags
Definition: liblwgeom.h:293
double gbox_angular_height(const GBOX *gbox)
GBOX utility functions to figure out coverage/location on the globe.
Definition: lwgeodetic.c:179
uint8_t flags
Definition: liblwgeom.h:385
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:206
Here is the call graph for this function:
Here is the caller graph for this function: