geography_bestsrid()

Datum geography_bestsrid

(

PG_FUNCTION_ARGS

)

Definition at line 786 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().

{
        GBOX gbox, gbox1, gbox2;
        GSERIALIZED *g1 = NULL;
        GSERIALIZED *g2 = NULL;
        int empty1 = LW_FALSE;
        int empty2 = LW_FALSE;
        double xwidth, ywidth;
        POINT2D center;

        Datum d1 = PG_GETARG_DATUM(0);
        Datum d2 = PG_GETARG_DATUM(1);

        /* Get our geometry objects loaded into memory. */
        g1 = (GSERIALIZED*)PG_DETOAST_DATUM(d1);
        /* Synchronize our box types */
        gbox1.flags = g1->flags;
        /* Calculate if the geometry is empty. */
        empty1 = gserialized_is_empty(g1);
        /* Calculate a geocentric bounds for the objects */
        if ( ! empty1 && gserialized_get_gbox_p(g1, &gbox1) == LW_FAILURE )
                elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g1, &gbox1)");

        POSTGIS_DEBUGF(4, "calculated gbox = %s", gbox_to_string(&gbox1));

        /* If we have a unique second argument, fill in all the necessary variables. */
        if ( d1 != d2 )
        {
                g2 = (GSERIALIZED*)PG_DETOAST_DATUM(d2);
                gbox2.flags = g2->flags;
                empty2 = gserialized_is_empty(g2);
                if ( ! empty2 && gserialized_get_gbox_p(g2, &gbox2) == LW_FAILURE )
                        elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g2, &gbox2)");
        }
        /*
        ** If no unique second argument, copying the box for the first
        ** argument will give us the right answer for all subsequent tests.
        */
        else
        {
                gbox = gbox2 = gbox1;
        }

        /* Both empty? We don't have an answer. */
        if ( empty1 && empty2 )
                PG_RETURN_NULL();

        /* One empty? We can use the other argument values as infill. Otherwise merge the boxen */
        if ( empty1 )
                gbox = gbox2;
        else if ( empty2 )
                gbox = gbox1;
        else
                gbox_union(&gbox1, &gbox2, &gbox);

        gbox_centroid(&gbox, &center);

        /* Width and height in degrees */
        xwidth = 180.0 * gbox_angular_width(&gbox)  / M_PI;
        ywidth = 180.0 * gbox_angular_height(&gbox) / M_PI;

        POSTGIS_DEBUGF(2, "xwidth %g", xwidth);
        POSTGIS_DEBUGF(2, "ywidth %g", ywidth);
        POSTGIS_DEBUGF(2, "center POINT(%g %g)", center.x, center.y);

        /* Are these data arctic? Lambert Azimuthal Equal Area North. */
        if ( center.y > 70.0 && ywidth < 45.0 )
        {
                PG_RETURN_INT32(SRID_NORTH_LAMBERT);
        }

        /* Are these data antarctic? Lambert Azimuthal Equal Area South. */
        if ( center.y < -70.0 && ywidth < 45.0 )
        {
                PG_RETURN_INT32(SRID_SOUTH_LAMBERT);
        }

        /*
        ** Can we fit these data into one UTM zone?
        ** We will assume we can push things as
        ** far as a half zone past a zone boundary.
        ** Note we have no handling for the date line in here.
        */
        if ( xwidth < 6.0 )
        {
                int zone = floor((center.x + 180.0) / 6.0);

                if ( zone > 59 ) zone = 59;

                /* Are these data below the equator? UTM South. */
                if ( center.y < 0.0 )
                {
                        PG_RETURN_INT32( SRID_SOUTH_UTM_START + zone );
                }
                /* Are these data above the equator? UTM North. */
                else
                {
                        PG_RETURN_INT32( SRID_NORTH_UTM_START + zone );
                }
        }

        /*
        ** 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 minimize 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(center.y / 30.0); /* (range of 0-5) */

                /* Equatorial band, 12 zones, 30 degrees wide */
                if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 )
                {
                        xzone = 6 + floor(center.x / 30.0);
                }
                /* Temperate band, 8 zones, 45 degrees wide */
                else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 )
                {
                        xzone = 4 + floor(center.x / 45.0);
                }
                /* Arctic band, 4 zones, 90 degrees wide */
                else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 )
                {
                        xzone = 2 + floor(center.x / 90.0);
                }

                /* Did we fit into an appropriate xzone? */
                if ( xzone != -1 )
                {
                        PG_RETURN_INT32(SRID_LAEA_START + 20 * yzone + xzone);
                }
        }

        /*
        ** Running out of options... fall-back to Mercator
        ** and hope for the best.
        */
        PG_RETURN_INT32(SRID_WORLD_MERCATOR);

}

Here is the call graph for this function:

Here is the caller graph for this function: