[postgis-users] function to find UTM zone SRID from POINT

[Rubén Ruiz]

>From BBN's OpenMap

/**
     * Converts a set of Longitude and Latitude co-ordinates to UTM
     * given an ellipsoid
     *
     * @param ellip an ellipsoid definition.
     * @param llpoint the coordinate to be converted
     * @param utmpoint A UTMPoint instance to put the results in. If
     *        null, a new UTMPoint will be allocated.
     * @return A UTM class instance containing the value of
     *         <code>null</code> if conversion failed. If you pass
     *         in a UTMPoint, it will be returned as well if
     *         successful.
     */
    public static UTMPoint LLtoUTM(LatLonPoint llpoint, Ellipsoid ellip,
                                   UTMPoint utmpoint) {

        double Lat = llpoint.getLatitude();
        double Long = llpoint.getLongitude();
        double a = ellip.radius;
        double eccSquared = ellip.eccsq;
        double k0 = 0.9996;

        double LongOrigin;
        double eccPrimeSquared;
        double N, T, C, A, M;

        double LatRad = llpoint.radlat_;
        double LongRad = llpoint.radlon_;
        double LongOriginRad;
        int ZoneNumber;

        ZoneNumber = (int) ((Long + 180) / 6) + 1;

        //Make sure the longitude 180.00 is in Zone 60
        if (Long == 180) {
            ZoneNumber = 60;
        }

        // Special zone for Norway
        if (Lat >= 56.0f && Lat < 64.0f && Long >= 3.0f && Long < 12.0f) {
            ZoneNumber = 32;
        }

        // Special zones for Svalbard
        if (Lat >= 72.0f && Lat < 84.0f) {
            if (Long >= 0.0f && Long < 9.0f)
                ZoneNumber = 31;
            else if (Long >= 9.0f && Long < 21.0f)
                ZoneNumber = 33;
            else if (Long >= 21.0f && Long < 33.0f)
                ZoneNumber = 35;
            else if (Long >= 33.0f && Long < 42.0f)
                ZoneNumber = 37;
        }
        LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3; //+3 puts origin
                                                     // in middle of
                                                     // zone
        LongOriginRad = ProjMath.degToRad(LongOrigin);

        eccPrimeSquared = (eccSquared) / (1 - eccSquared);

        N = a / Math.sqrt(1 - eccSquared * Math.sin(LatRad) *
Math.sin(LatRad));
        T = Math.tan(LatRad) * Math.tan(LatRad);
        C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad);
        A = Math.cos(LatRad) * (LongRad - LongOriginRad);

        M = a
                * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 -
5
                        * eccSquared * eccSquared * eccSquared / 256)
                        * LatRad
                        - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared
                                / 32 + 45 * eccSquared * eccSquared
                                * eccSquared / 1024)
                        * Math.sin(2 * LatRad)
                        + (15 * eccSquared * eccSquared / 256 + 45 *
eccSquared
                                * eccSquared * eccSquared / 1024)
                        * Math.sin(4 * LatRad) - (35 * eccSquared *
eccSquared
                        * eccSquared / 3072)
                        * Math.sin(6 * LatRad));

        double UTMEasting = (double) (k0
                * N
                * (A + (1 - T + C) * A * A * A / 6.0d + (5 - 18 * T + T * T
                        + 72 * C - 58 * eccPrimeSquared)
                        * A * A * A * A * A / 120.0d) + 500000.0d);

        double UTMNorthing = (double) (k0 * (M + N
                * Math.tan(LatRad)
                * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A
                        / 24.0d + (61 - 58 * T + T * T + 600 * C - 330 *
eccPrimeSquared)
                        * A * A * A * A * A * A / 720.0d)));
        if (Lat < 0.0f) {
            UTMNorthing += 10000000.0f; //10000000 meter offset for
                                        // southern hemisphere
        }

        if (utmpoint == null) {
            utmpoint = new UTMPoint();
        }

        utmpoint.northing = (double) Math.rint(UTMNorthing);
        utmpoint.easting = (double) Math.rint(UTMEasting);
        utmpoint.zone_number = ZoneNumber;
        utmpoint.zone_letter = utmpoint.getLetterDesignator(Lat);

        return utmpoint;
    }

-----Mensaje original-----
De: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] En nombre de Alex
Mayrhofer
Enviado el: sábado, 03 de diciembre de 2005 14:09
Para: PostGIS Users Discussion
Asunto: [postgis-users] function to find UTM zone SRID from POINT

All,

Is there a function/table to find out the SRID of the UTM zone where a given

point falls into? This should be a fairly frequent problem - i'd like to 
avoid reinventing the wheel ;)

thanks,

Alex Mayrhofer
_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users