PostGIS  3.7.0dev-r@@SVN_REVISION@@

◆ geography_bestsrid()

Datum geography_bestsrid ( PG_FUNCTION_ARGS  )

Definition at line 794 of file geography_measurement.c.

795 {
796  GBOX gbox, gbox1, gbox2;
797  GSERIALIZED *g1 = NULL;
798  GSERIALIZED *g2 = NULL;
799  int empty1 = LW_FALSE;
800  int empty2 = LW_FALSE;
801  double xwidth, ywidth;
802  POINT2D center;
803  LWGEOM *lwgeom;
804 
805  /* Get our geometry objects loaded into memory. */
806  g1 = PG_GETARG_GSERIALIZED_P(0);
807  /* Synchronize our box types */
808  gbox1.flags = gserialized_get_lwflags(g1);
809  /* Calculate if the geometry is empty. */
810  empty1 = gserialized_is_empty(g1);
811 
812  /* Convert g1 to LWGEOM type */
813  lwgeom = lwgeom_from_gserialized(g1);
814 
815  /* Calculate a geocentric bounds for the objects */
816  if ( ! empty1 && gserialized_get_gbox_p(g1, &gbox1) == LW_FAILURE )
817  elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g1, &gbox1)");
818 
819  POSTGIS_DEBUGF(4, "calculated gbox = %s", gbox_to_string(&gbox1));
820 
821  if ( !lwgeom_isfinite(lwgeom) ) {
822  elog(ERROR, "Error in geography_bestsrid calling with infinite coordinate geographies");
823  }
824  lwgeom_free(lwgeom);
825 
826  /* If we have a unique second argument, fill in all the necessary variables. */
827  if (PG_NARGS() > 1)
828  {
829  g2 = PG_GETARG_GSERIALIZED_P(1);
830  gbox2.flags = gserialized_get_lwflags(g2);
831  empty2 = gserialized_is_empty(g2);
832  if ( ! empty2 && gserialized_get_gbox_p(g2, &gbox2) == LW_FAILURE )
833  elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g2, &gbox2)");
834 
835  /* Convert g2 to LWGEOM type */
836  lwgeom = lwgeom_from_gserialized(g2);
837 
838  if ( !lwgeom_isfinite(lwgeom) ) {
839  elog(ERROR, "Error in geography_bestsrid calling with second arg infinite coordinate geographies");
840  }
841  lwgeom_free(lwgeom);
842  }
843  /*
844  ** If no unique second argument, copying the box for the first
845  ** argument will give us the right answer for all subsequent tests.
846  */
847  else
848  {
849  gbox = gbox2 = gbox1;
850  }
851 
852  /* Both empty? We don't have an answer. */
853  if ( empty1 && empty2 )
854  PG_RETURN_NULL();
855 
856  /* One empty? We can use the other argument values as infill. Otherwise merge the boxen */
857  if ( empty1 )
858  gbox = gbox2;
859  else if ( empty2 )
860  gbox = gbox1;
861  else
862  gbox_union(&gbox1, &gbox2, &gbox);
863 
864  gbox_centroid(&gbox, &center);
865 
866  /* Width and height in degrees */
867  xwidth = 180.0 * gbox_angular_width(&gbox) / M_PI;
868  ywidth = 180.0 * gbox_angular_height(&gbox) / M_PI;
869 
870  POSTGIS_DEBUGF(2, "xwidth %g", xwidth);
871  POSTGIS_DEBUGF(2, "ywidth %g", ywidth);
872  POSTGIS_DEBUGF(2, "center POINT(%g %g)", center.x, center.y);
873 
874  /* Are these data arctic? Lambert Azimuthal Equal Area North. */
875  if ( center.y > 70.0 && ywidth < 45.0 )
876  {
877  PG_RETURN_INT32(SRID_NORTH_LAMBERT);
878  }
879 
880  /* Are these data antarctic? Lambert Azimuthal Equal Area South. */
881  if ( center.y < -70.0 && ywidth < 45.0 )
882  {
883  PG_RETURN_INT32(SRID_SOUTH_LAMBERT);
884  }
885 
886  /*
887  ** Can we fit these data into one UTM zone?
888  ** We will assume we can push things as
889  ** far as a half zone past a zone boundary.
890  ** Note we have no handling for the date line in here.
891  */
892  if ( xwidth < 6.0 )
893  {
894  int zone = floor((center.x + 180.0) / 6.0);
895 
896  if ( zone > 59 ) zone = 59;
897 
898  /* Are these data below the equator? UTM South. */
899  if ( center.y < 0.0 )
900  {
901  PG_RETURN_INT32( SRID_SOUTH_UTM_START + zone );
902  }
903  /* Are these data above the equator? UTM North. */
904  else
905  {
906  PG_RETURN_INT32( SRID_NORTH_UTM_START + zone );
907  }
908  }
909 
910  /*
911  ** Can we fit into a custom LAEA area? (30 degrees high, variable width)
912  ** We will allow overlap into adjoining areas, but use a slightly narrower test (25) to try
913  ** and minimize the worst case.
914  ** Again, we are hoping the dateline doesn't trip us up much
915  */
916  if ( ywidth < 25.0 )
917  {
918  int xzone = -1;
919  int yzone = 3 + floor(center.y / 30.0); /* (range of 0-5) */
920 
921  /* Equatorial band, 12 zones, 30 degrees wide */
922  if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 )
923  {
924  xzone = 6 + floor(center.x / 30.0);
925  }
926  /* Temperate band, 8 zones, 45 degrees wide */
927  else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 )
928  {
929  xzone = 4 + floor(center.x / 45.0);
930  }
931  /* Arctic band, 4 zones, 90 degrees wide */
932  else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 )
933  {
934  xzone = 2 + floor(center.x / 90.0);
935  }
936 
937  /* Did we fit into an appropriate xzone? */
938  if ( xzone != -1 )
939  {
940  PG_RETURN_INT32(SRID_LAEA_START + 20 * yzone + xzone);
941  }
942  }
943 
944  /*
945  ** Running out of options... fall-back to Mercator
946  ** and hope for the best.
947  */
948  PG_RETURN_INT32(SRID_WORLD_MERCATOR);
949 
950 }
int gbox_union(const GBOX *g1, const GBOX *g2, GBOX *gout)
Update the output GBOX to be large enough to include both inputs.
Definition: gbox.c:135
char * gbox_to_string(const GBOX *gbox)
Allocate a string representation of the GBOX, based on dimensionality of flags.
Definition: gbox.c:404
int gserialized_get_gbox_p(const GSERIALIZED *g, GBOX *gbox)
Read the box from the GSERIALIZED or calculate it if necessary.
Definition: gserialized.c:94
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:268
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: gserialized.c:181
lwflags_t gserialized_get_lwflags(const GSERIALIZED *g)
Read the flags from a GSERIALIZED and return a standard lwflag integer.
Definition: gserialized.c:47
#define LW_FALSE
Definition: liblwgeom.h:94
#define LW_FAILURE
Definition: liblwgeom.h:96
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1218
int lwgeom_isfinite(const LWGEOM *lwgeom)
Check if a LWGEOM has any non-finite (NaN or Inf) coordinates.
Definition: lwgeom.c:2807
double gbox_angular_height(const GBOX *gbox)
GBOX utility functions to figure out coverage/location on the globe.
Definition: lwgeodetic.c:188
double gbox_angular_width(const GBOX *gbox)
Returns the angular width (longitudinal span) of the box in radians.
Definition: lwgeodetic.c:215
int gbox_centroid(const GBOX *gbox, POINT2D *out)
Computes the average(ish) center of the box and returns success.
Definition: lwgeodetic.c:267
lwflags_t flags
Definition: liblwgeom.h:353
double y
Definition: liblwgeom.h:390
double x
Definition: liblwgeom.h:390

References GBOX::flags, gbox_angular_height(), gbox_angular_width(), gbox_centroid(), gbox_to_string(), gbox_union(), gserialized_get_gbox_p(), gserialized_get_lwflags(), gserialized_is_empty(), LW_FAILURE, LW_FALSE, lwgeom_free(), lwgeom_from_gserialized(), lwgeom_isfinite(), POINT2D::x, and POINT2D::y.

Here is the call graph for this function: