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
Leave a Reply