PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ geography_bestsrid()

Datum geography_bestsrid ( PG_FUNCTION_ARGS  )

Definition at line 815 of file geography_measurement.c.

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

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: