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?
