Convert Feet to Decimal Degrees

Question

I need to convert a bunch of coordinates that were stored in Feet to Decimal Degrees. I found a Function in VB.Net, though it's results are not coming back correctly. When I input 3106972, 1417936 into this function I get: -80.7964701731349, 32.8465039721407 where they should be: -105.115621, 40.480104. Where am I going wrong here?

Here is the function:

Private Function ConvertFeetToDecimalDegrees(lat As Integer, lng As Integer) As String
' This function converts coordinates from Florida State Plane North Zone NAD83 feet  
' to Geographic coordinates in decimal degrees  
' Set up the coordinate system parameters.  
Dim A = 20925604.47                         ' major radius of GRS 1980 ellipsoid, feet  
Dim Ec = 0.0818191910435                    ' eccentricity of GRD 1980 ellipsoid  
Dim Ec2 = Ec * Ec                           ' eccentricity squared  
Dim AngRad = 0.0174532925199433             ' number of radians in a degree  
Dim Pi4 = 3.141592653582 / 4                ' Pi / 4  
Dim P1 = 29.583333333333329 * AngRad        ' latitude of first standard parallel  
Dim P2 = 30.75 * AngRad                     ' latitude of second standard parallel  
Dim P0 = 29.0 * AngRad                      ' latitude of origin  
Dim M0 = -84.5 * AngRad                     ' central meridian  
Dim X0 = 1968500.0                          ' False easting of central meridian, map units  

' Calculate the coordinate system constants.  
Dim m1 = Math.Cos(P1) / Math.Sqrt(1 - (Ec2 * Math.Pow((Math.Sin(P1)), 2)))
Dim m2 = Math.Cos(P2) / Math.Sqrt(1 - (Ec2 * Math.Pow((Math.Sin(P2)), 2)))
Dim t1 = Math.Tan(Pi4 - (P1 / 2)) / Math.Pow((1 - Ec * Math.Sin(P1)) / (1 + Ec * Math.Sin(P1)), (Ec / 2))
Dim t2 = Math.Tan(Pi4 - (P2 / 2)) / Math.Pow((1 - Ec * Math.Sin(P2)) / (1 + Ec * Math.Sin(P2)), (Ec / 2))
Dim t0 = Math.Tan(Pi4 - (P0 / 2)) / Math.Pow((1 - Ec * Math.Sin(P0)) / (1 + Ec * Math.Sin(P0)), (Ec / 2))
Dim n = Math.Log(m1 / m2) / Math.Log(t1 / t2)
Dim F = m1 / (n * Math.Pow(t1, n))
Dim rho0 = A * F * Math.Pow(t0, n)
' Convert the coordinate to Latitude/Longitude.  

' Calculate the Longitude.  
lat = lat - X0
Dim Pi2 = Pi4 * 2
Dim rho = Math.Sqrt(Math.Pow(lat, 2) + (Math.Pow(rho0 - lng, 2)))
Dim theta = Math.Atan(lat / (rho0 - lng))
Dim t = Math.Pow(rho / (A * F), (1 / n))
Dim LonR = theta / n + M0
lat = lat + X0

' Estimate the Latitude  
Dim Lat0 = Pi2 - (2 * Math.Atan(t))

' Substitute the estimate into the iterative calculation that  
' converges on the correct Latitude value.  
Dim part1 = (1 - (Ec * Math.Sin(Lat0))) / (1 + (Ec * Math.Sin(Lat0)))
Dim LatR = Pi2 - (2 * Math.Atan(t * Math.Pow(part1, (Ec / 2))))

Lat0 = LatR
part1 = (1 - (Ec * Math.Sin(Lat0))) / (1 + (Ec * Math.Sin(Lat0)))
LatR = Pi2 - (2 * Math.Atan(t * Math.Pow(part1, (Ec / 2))))
Debug.WriteLine(LatR & "|" & LonR)

If (Math.Abs(LatR - Lat0) > 0.000000002) Then
  ' Convert from radians to degrees.  
  Dim _Lat = LatR / AngRad
  Dim _Lng = LonR / AngRad
  Return _Lng.ToString + "|" + _Lat.ToString
End If
Return ""
End Function

Answer 1

There is a similar formula described in this Esri forum post (for Washington State rather than Florida), so it is certainly a problem that's been addressed for a range of locations. The post refers to Map Projections, A Working Manual for the mathematical basics to modify the script for a different location.

However, that would require some familiarity with projections, a fair amount of reading, and knowing the specific projection of the original data. If you don't know it, the question becomes much trickier and there is a chance of significant errors being introduced. The error will be less than "somewhere in Florida when it's supposed to be in Colorado" of course, but whether that's acceptable depends on the application (are you locating cities or utility poles?)...

Another option would be to explore scripting through a GIS program (e.g., QGIS which is free) and using their pre-existing transformation functions. Again, this depends on knowing the original projection of your data. If you choose to go that route, then scripting the function would need to be a new Question (because it's pretty different, and because I don't know how to answer it.)

Answer 2

The feet and degree coordinates fit to each other if you use EPSG:2231 NAD83 / Colorado North (ftUS) with the following parameters:

+proj=lcc +lat_1=40.78333333333333 +lat_2=39.71666666666667 +lat_0=39.33333333333334 +lon_0=-105.5 +x_0=914401.8288036576 +y_0=304800.6096012192 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

Using QGIS measurement tool, the offset is about 21m. I added the point for a NAD27 CRS as well, that is 85m West of NAD83/WGS84.