# Mathematical formula to convert lat-lon to Lambert Conformal Conic

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:
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
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.
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`?

So I've figured it out.

I had the right formula. To validate that, the projection wkt gives:

``````METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
``````

Which fits with the manual at page 19

Than, I validated my equations with their example (which is projection EPSG:32040) to figure that I used `e^2` instead of `e/2` in my `t` calculation (dumb, I know). Changing this fixed my problem.

The final script in case somebody wants to reproduce this:

``````## 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

## The information about the wanted projection
sf::st_crs(3347)

## Which we transform into the variables and transform in Radian instead of degree
lat1 = 49 * pi / 180
lat2 = 77 * pi / 180
latF = 63.390675 * pi / 180
lon0 = -91.8666666666667 * pi / 180
EF = 6200000
NF = 3000000

a = 6378137
f = 1 / 298.257222101
e = sqrt(2*f - f^2)

## Adaption of 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)

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)

E = EF + r * sin(theta)
N = NF + rF - r * cos(theta)

cbind(E, N)

sf::st_as_sf(coo, coords = c("lon", "lat"), remove=F, crs = st_crs(4326)) |>
sf::st_transform(st_crs(3347)) |>
sf::st_coordinates()
``````