I want to manually convert coordinates (latitude-longitude in degree decimals) data into Lambert Conformal Conic (more precisely, EPSG:3347). I already know how to use GIS software functions for this (st_transform in R, ogr2ogr in gdal, using directly QGIS, etc.). All my search online points to those specific functions.
What I need is the actual equation used to do the conversion so I can program my own function in another software (SAS, Excel, etc) that is less GIS friendly.
I've found this link that refers to this manual which I tried to reproduce in R, with no success. I'm not even sure it's to good formula.
## The coordinates to reproject coo = data.frame(lat=c(46.313503, 44.336374, 50.174041, 50.955616, 48.221700), lon=c(-72.663290, -77.936728, -100.085165, -124.870320, -56.842978)) lat = coo$lat *pi / 180 lon = coo$lon *pi / 180 Then extract the information about the projection I want.
sf::st_crs(3347) Coordinate Reference System: User input: EPSG:3347 wkt: PROJCRS["NAD83 / Statistics Canada Lambert", BASEGEOGCRS["NAD83", DATUM["North American Datum 1983", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4269]], CONVERSION["Statistics Canada Lambert", METHOD["Lambert Conic Conformal (2SP)", ID["EPSG",9802]], PARAMETER["Latitude of false origin",63.390675, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8821]], PARAMETER["Longitude of false origin",-91.8666666666667, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8822]], PARAMETER["Latitude of 1st standard parallel",49, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8823]], PARAMETER["Latitude of 2nd standard parallel",77, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8824]], PARAMETER["Easting at false origin",6200000, LENGTHUNIT["metre",1], ID["EPSG",8826]], PARAMETER["Northing at false origin",3000000, LENGTHUNIT["metre",1], ID["EPSG",8827]]], CS[Cartesian,2], AXIS["(E)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["(N)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["Topographic mapping (small scale)."], AREA["Canada - onshore and offshore - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon."], BBOX[38.21,-141.01,86.46,-40.73]], ID["EPSG",3347]] Which we transform into R variables and transform in Radian instead of degree
lat1 = 49 * pi / 180 lat2 = 77 * pi / 180 latF = 63.390675 * pi / 180 # is it the latitude of the false origin? not sure. lon0 = -91.8666666666667 * pi / 180 # Not sure about this one. EF = 6200000 NF = 3000000 Other variables I need to extract. Their description was at page 6 of the manual cited above.
a = 6378137 f = 1 / 298.257222101 e = sqrt(2*f - f^2) Then, we adapt the formula in the website
m1 = cos(lat1)/(1 - e^2 * sin(lat1)^2)^0.5 # for m1, lat1, and m2, lat2 where lat1 and lat2 are the latitudes of the two standard parallels. m2 = cos(lat2)/(1 - e^2 * sin(lat2)^2)^0.5 t1 = tan(pi/4 - lat1/2)/((1 - e * sin(lat1))/(1 + e * sin(lat1)))^(e^2) # for t1, t2, tF and t using lat1, lat2, latF and lat respectively. t2 = tan(pi/4 - lat2/2)/((1 - e * sin(lat2))/(1 + e * sin(lat2)))^(e^2) tF = tan(pi/4 - latF/2)/((1 - e * sin(latF))/(1 + e * sin(latF)))^(e^2) t = tan(pi/4 - lat/2)/((1 - e * sin(lat))/(1 + e * sin(lat)))^(e^2) I'm not sure I'm using the right latF and lat
Then:
n = (log(m1) - log(m2))/(log(t1) - log(t2)) FF = m1/(n * t1^n) r = a * FF * t^n # for rF and r, where rF is the radius of the parallel of latitude of the false origin. rF = a * FF * tF^n theta = n * (lon - lon0) Here, I'm confuse about lon0 is it the Longitude of false origin? And lon is it the lon I want to convert?
Finally:
E = EF + r * sin(theta) N = NF + rF - r * cos(theta) Running this produce:
cbind(E, N) E N [1,] 7673386 1349263 [2,] 7325110 1025063 [3,] 5617116 1594762 [4,] 3999380 2227177 [5,] 8680877 2039659 If I compare to what st_transform gives me:
sf::st_as_sf(coo, coords = c("lon", "lat"), remove=F, crs = st_crs(4326)) |> sf::st_transform(st_crs(3347)) |> sf::st_coordinates() X Y [1,] 7673199 1352744 [2,] 7324899 1029094 [3,] 5617097 1596987 [4,] 3999454 2229785 [5,] 8680546 2043038 I'm not quite there yet. It's close, but not the same.
In summary, what is the mathematical procedure to convert degree-decimals coordinates to lcc. why doesn't the process I did is not producing similar values than st_transform?
