Within the infrastructure where I work there are several different applications. One application was mapping via lat/long and the other by state plane. In order to properly pinpoint locations on a map I needed a quick way to go from state plane and back again. This was a trying effort. There are numerous map projections and coordinate systems and to someone who has never delved into this area it can be quite daunting.

Furthermore there is a great deal of information stretching back literally hundreds of years dealing with algorithms for finding locations and translating projections and coordinate systems. I found all of these to be either indecipherable, incomplete, incorrect, or too inaccurate/with too wide a margin of error. While the USGS site is a wonderful resource, the documentation available on the site is multitudinous but ungraded. A document stretching back to the early 20th century would be listed beside a document dealing with satellite GPS tracking.

Another challenge was finding the right language which could be easily integrated into systems I was using. Initially I tried running the algorithms in pure TSQL thinking the decimal type would suffice but TSQL is not designed for precision to the number of decimal places I needed. I ended up using a Visual C# SQL CLR Database Project since both the GIS database and the application database were on the same server. Microsoft SQL server also contains functions for dealing with spatial data which would be important when I tried to match coordinates to parcels.

For testing my results I used the Earth Point website.
The basis of the final algorithm was extracted from the SPCS83 program provided by the National Geodetic Survey .
To test and complie the FORTRAN code I used MinGW.

Function COSStatePlaneToLatorLong

This is the basic algorith extracted from various sources. Ultimately I tried to follow the SPCS83 program as closely as possible. It takes in the easting, northing, latitude or longitude depending on what we are trying to get, and the unit since a survey foot differs from an international foot.

using System;
using System.Data.SqlTypes;
using System.Globalization;


namespace COSCLRFunctions
{
    public partial class UserDefinedFunctions
    {         
        [Microsoft.SqlServer.Server.SqlFunction]
        public static String COSStatePlaneToLatorLong(SqlString Easting, SqlString Northing, SqlString latorlong, SqlString unit) //in survey feet
        {
            decimal eastingInMeters;
            decimal northingInMeters;

            if (unit == "sfeet")
            {
                eastingInMeters = (1200.0m / 3937.0m) * (decimal)Easting.ToSqlDecimal();
                northingInMeters = (1200.0m / 3937.0m) * (decimal)Northing.ToSqlDecimal();
            }
            else
            {
                eastingInMeters = (decimal)Easting.ToSqlDecimal();
                northingInMeters = (decimal)Northing.ToSqlDecimal();
            }

            const decimal rad = 57.295779513082323m; //RadianToDegree
            const decimal er = 6378137.0000000000m;  //equatorial radius of the ellipsoid (= major semiaxis).
            const decimal esq = 6.6943800229007869E-003m; //the square of first eccentricity of the ellipsoid.
            decimal e = (decimal)Math.Sqrt((double)esq);
            const decimal cm = -71.5m / rad;                   //Central Meridian
            const decimal sp1 = 41.716666666666669m; //Standard Parallel 1
            const decimal sp2 = 42.683333333333330m; //Standard Parallel 2
            const decimal nb = 750000;               //False Northing
            const decimal eo = 200000;               //False Easting
            const decimal rl = 41.000000000000000m;  //Reference Latitude

            decimal fis = sp1 / rad;
            decimal sinFs = (decimal)Math.Sin((double)fis);
            decimal cosFs = (decimal)Math.Cos((double)fis);

            decimal fin = sp2 / rad;
            decimal sinFn = (decimal)Math.Sin((double)fin);
            decimal cosFn = (decimal)Math.Cos((double)fin);

            decimal fib = rl / rad;
            decimal sinFb = (decimal)Math.Sin((double)fib);

            decimal qs = GetQ(e, sinFs);
            decimal qn = GetQ(e, sinFn);
            decimal qb = GetQ(e, sinFb);

            decimal w1 = (decimal)Math.Sqrt(1 - (double)esq * (double)(sinFs * sinFs));

            decimal w2 = (decimal)Math.Sqrt(1 - (double)esq * (double)(sinFn * sinFn));

            decimal sinFo = (decimal)Math.Log((double)w2 * (double)cosFs / (double)(w1 * cosFn)) / (qn - qs);

            decimal k = (decimal)((double)er * (double)cosFs * Math.Exp((double)qs * (double)sinFo) / (double)(w1 * sinFo));

            decimal rb = (decimal)((double)k / Math.Exp((double)qb * (double)sinFo));

            decimal npr = rb - northingInMeters + nb;

            decimal epr = eastingInMeters - eo;

            decimal gam = (decimal)Math.Atan((double)(epr / npr));

            decimal lon = (decimal)(Math.Abs((double)cm) - (double)gam / (double)sinFo);


            decimal rpt = (decimal)Math.Sqrt((double)(npr * npr + epr * epr));

            decimal q = (decimal)(Math.Log((double)(k / rpt)) / (double)sinFo);

            decimal temp = (decimal)Math.Exp((double)(q + q));

            decimal sine = (temp - 1) / (temp + 1);

            double i = 1;

            while (i <= 3)
            {
                var f1 = (decimal)(Math.Log((double)((1 + sine) / (1 - sine))) - (double)e * Math.Log((double)((1 + e * sine) / (1 - e * sine)))) / 2 - q;

                var f2 = 1 / (1 - sine * sine) - esq / (1 - esq * sine * sine);

                sine = sine - f1 / f2;

                i++;
            }

            decimal lat = (decimal)Math.Asin((double)sine);
            decimal fi = lat;

            fi = ConvertRadiansToDegrees(fi);

            if (latorlong == "latitude")
            {
                return fi.ToString(CultureInfo.CurrentCulture);
            }
            else
            {
                return ConvertRadiansToDegrees(-1 * lon).ToString(CultureInfo.CurrentCulture); 
            }
        }

Function COSLatLongToStatePlane

This does the reverse of the above function.

[Microsoft.SqlServer.Server.SqlFunction]
        public static SqlString COSLatLongToStatePlane(SqlString latitude, SqlString longitude, SqlString eastingornorthing, SqlString units) //in survey feet
        {
            const decimal rad = 57.295779513082323m;                    //RadianToDegree
            const decimal er = 6378137.0000000000m;                     //equatorial radius of the ellipsoid (= major semiaxis).
            const decimal esq = 6.6943800229007869E-003m;               //the square of first eccentricity of the ellipsoid.
            const decimal cm = 71.5m / rad;                             //Central Meridian
            const decimal sp1 = 41.716666666666669m;                    //Standard Parallel 1
            const decimal sp2 = 42.683333333333330m;                    //Standard Parallel 2
            const decimal nb = 750000;                                  //False Northing
            const decimal eo = 200000;                                  //False Easting
            const decimal rl = 41.000000000000000m;                     //Reference Latitude
            decimal e = (decimal)Math.Sqrt(6.6943800229007869E-003);


            decimal fis = sp1 / rad;
            decimal sinFs = (decimal)Math.Sin((double)fis);
            decimal cosFs = (decimal)Math.Cos((double)fis);

            decimal fin = sp2 / rad;
            decimal sinFn = (decimal)Math.Sin((double)fin);
            decimal cosFn = (decimal)Math.Cos((double)fin);

            decimal fib = rl / rad;
            decimal sinFb = (decimal)Math.Sin((double)fib);

            decimal qs = GetQ(e, sinFs);
            decimal qn = GetQ(e, sinFn);
            decimal qb = GetQ(e, sinFb);

            decimal w1 = (decimal)Math.Sqrt(1 - (double)esq * (double)(sinFs * sinFs));
            decimal w2 = (decimal)Math.Sqrt(1 - (double)esq * (double)(sinFn * sinFn));

            decimal sinFo = (decimal)Math.Log((double)w2 * (double)cosFs / (double)(w1 * cosFn)) / (qn - qs);

            decimal k = (decimal)((double)er * (double)cosFs * Math.Exp((double)qs * (double)sinFo) / (double)(w1 * sinFo));

            decimal rb = (decimal)((double)k / Math.Exp((double)qb * (double)sinFo));

            decimal fi = ConvertDegreesToRadians((decimal)latitude.ToSqlDecimal());
            decimal lam = ConvertDegreesToRadians(Math.Abs((decimal)longitude.ToSqlDecimal()));

            decimal sinLat = (decimal)Math.Sin((double)fi);
            decimal conv = (decimal)(((double)cm - (double)lam) * (double)sinFo);

            var q = (decimal)((Math.Log((1 + (double)sinLat) / (1 - (double)sinLat)) -
                               (double)e * Math.Log((1 + (double)e * (double)sinLat) / (1 - (double)e * (double)sinLat))) / 2);
            var rpt = (decimal)((double)k / Math.Exp((double)sinFo * (double)q));
            var north = (decimal)((double)nb + (double)rb - (double)rpt * Math.Cos((double)conv));
            var east = (decimal)((double)eo + (double)rpt * Math.Sin((double)conv));

            if (units == "sfeet")
            {
                east =  east / (1200.0m / 3937.0m);
                north = north / (1200.0m / 3937.0m);
            }

            if (eastingornorthing.ToString() == "easting")
            {
                return east.ToString(CultureInfo.CurrentCulture);
            }
            else
            {
                return north.ToString(CultureInfo.CurrentCulture);
            }

        }

Helper Functions

Various helper functions.

        private static decimal GetQ(decimal E, decimal S)
        {
            decimal Q = (decimal)(Math.Log((1 + (double)S) / (1 - (double)S)) -
                                  (double)E * Math.Log((1 + (double)E * (double)S) / (1 - (double)E * (double)S))) / 2;
            return Q;
        }
        public static decimal ConvertRadiansToDegrees(decimal radians)
        {
            decimal degrees = (decimal)((180 / Math.PI) * (double)radians);
            return ((decimal)degrees);
        }   
        public static decimal ConvertDegreesToRadians(decimal degrees)
        {
            var radians = 0.01745329252 * (double)degrees;
            return (decimal)radians;
        }
		public static Boolean COSIsValidPolygon(SqlString sqlgeometrystring)
        {
            Boolean result = false;
            sqlgeometrystring = sqlgeometrystring.ToString().Replace("(", "").Replace(")", "").Replace("POLYGON ", "").Replace(", ", ",");
            string firstentry = sqlgeometrystring.ToString().Substring(0, sqlgeometrystring.ToString().IndexOf(','));
            string lastentry = ReverseString(ReverseString(sqlgeometrystring.ToString()).Substring(0, sqlgeometrystring.ToString().IndexOf(',')));
            if (firstentry.Equals(lastentry))
            {
                result = true;
            }
            else
            {
                result = false;
            }
            return result;
        }
        public static string ReverseString(string s)
        {
            char[] arr = s.ToCharArray();
            Array.Reverse(arr);
            return new string(arr);
        }
[Microsoft.SqlServer.Server.SqlFunction]
        public static SqlString COSStatePlaneGeometrytoLatLongGeometry(SqlString sqlgeometrystring) //in survey feet
        {


            Boolean torf = true;
            String[] temparray = sqlgeometrystring.ToString().Split(new[] { "), (" }, StringSplitOptions.None);
            String[] result = new String[temparray.Length];

            for (int i = 0; i < temparray.Length; i++)
            {
				//Some of the polygons in the database do not have closing elements so are not real polygons
                if (COSIsValidPolygon(temparray[i]) != true)
                {
                    torf = false;
                    break;
                }
            }

			//if the polygon is good
            if (torf == true)
            {
                for (int z = 0; z < temparray.Length; z++)
                {
                    temparray[z] = temparray[z].ToString().Replace("(", "").Replace(")", "")
                        .Replace("POLYGON ", "")
                        .Replace(", ", ",");
                    String[] sqlgeometry_c_stringarray_source = temparray[z].ToString().Split(',');
                    result[z] = "(";
                    String[] sqlgeometry_c_stringarray_final = new String[sqlgeometry_c_stringarray_source.Length];
                    for (int i = 0; i < sqlgeometry_c_stringarray_source.Length; i++)
                    {

                        String temp1 = COSSataePlaneToLatorLong((String)sqlgeometry_c_stringarray_source[i].Split(' ')[0].Trim(),
                            (String)sqlgeometry_c_stringarray_source[i].Split(' ')[1].Trim(), (String)"easting",
                            (String)"sfeet");
                        String temp2 = COSSataePlaneToLatorLong((String)sqlgeometry_c_stringarray_source[i].Split(' ')[0].Trim(),
                            (String)sqlgeometry_c_stringarray_source[i].Split(' ')[1], (String)"northing",
                            (String)"sfeet");

                        result[z] = result[z] + temp1 + " " + temp2 + ", ";
                    }

                    char[] trimchars = { ',', ' ' };
                    result[z] = result[z].TrimEnd(trimchars);
                    result[z] = result[z] + ")";
                    String temp = sqlgeometry_c_stringarray_final.ToString();
                    //Console.WriteLine("result: " + result[z]);
                }
            }
        }

Usage

the basic usage for the functions follows. If all is correct the first two function calls should put out results that can be plugged into the next two function calls which will return the original lat/long.

SELECT top 1 
	dbo.COSLatLongToStatePlane(42.1013889,-72.5902778,'easting','') AS easting
	,dbo.COSLatLongToStatePlane(42.1013889,-72.5902778,'','') AS northing
	,dbo.COSSataePlaneToLatorLong(109818.706390295,872906.990453886,'latitude','') AS lat
	,dbo.COSSataePlaneToLatorLong(109818.706390295,872906.990453886,'','') AS long
FROM SOMETABLE A

Parcel from Coordinates

Finally, just for completion, the Scalar-valued functions to match to parcels.

USE [SOMEDATABASE]
GO

SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO

-- =============================================
-- Author:		<Author,,Name>
-- Create date: <Create Date, ,>
-- Description:	<Description, ,>
-- =============================================
ALTER FUNCTION [dbo].[COSSVFParcelFromLatLong] 
(
	@lat VARCHAR(50) = NULL, 
	@LONG VARCHAR(50) = NULL
)
RETURNS [nvarchar](max)
AS
BEGIN

	DECLARE @return_value varchar(50)

	SELECT @return_value = SP --easting + ',' + northing
	FROM [springfield].[gis].[PARCELS]

	OUTER APPLY(
		SELECT dbo.COSLatLongToStatePlane(convert(varchar(20),@lat),convert(varchar(20),@LONG),'easting','sfeet') as easting
			  ,dbo.COSLatLongToStatePlane(convert(varchar(20),@lat),convert(varchar(20),@LONG),'northing','sfeet') as northing
	) A

	OUTER APPLY(
		SELECT geometry::Point(easting,northing,2249) as geopoint
	) B

	WHERE [PARCELS].[Shape].STContains(geopoint) <> 0


RETURN @return_value
END
USE [SOMEDATABASE]
GO

SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO


-- =============================================
-- Author:		<Author,,Name>
-- Create date: <Create Date, ,>
-- Description:	<Description, ,>
-- =============================================
ALTER FUNCTION [dbo].[COSSVFParcelFromStatePlane] 
(
	@lat VARCHAR(50) = NULL, 
	@LONG VARCHAR(50) = NULL
)
RETURNS [nvarchar](max)
AS
BEGIN

	DECLARE @return_value varchar(50)

	SELECT @return_value = SP --easting + ',' + northing
	FROM [springfield].[gis].[PARCELS]

	OUTER APPLY(
		SELECT geometry::Point(@lat,@LONG,2249) as geopoint
	) B

	WHERE [PARCELS].[Shape].STContains(geopoint) <> 0


RETURN @return_value
END

 365 total views,  1 views today