[GIS] Converting Lambert Conformal (x, y) of spheroid to (lon, lat) of NAD83 with proj/cs2cs

cs2csellipsoidprojspherical-geometry

Revised Trial: Correct?

According to Andre Joost advice, I tried with "+towgs84=0,0,0", and the both two ways seem to get the same answer that mkennedy posted below. So, just to make sure, anybody can tell me whether this is correct? What would be an easy way to verify my results?

And the original two ways were done as follows:

The first way:

$ cs2cs -v -f '%.16f' +proj=lcc +no_defs +a=6370000.0m
+b=6370000.0m +lon_0=97w +lat_0=40n +lat_1=33n +lat_2=45n +\ x_0=900000.0m +y_0=1620000.0m +towgs84=0,0,0 +to +proj=lonlat
+datum=NAD83 +no_defs<<EOF
2340000.0 1800000.0 EOF
# ---- From Coordinate System ----
#Lambert Conformal Conic
#       Conic, Sph&Ell
#       lat_1= and lat_2= or lat_0
# +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n
# +lat_1=33n +lat_2=45n +x_0=900000.0m +y_0=1620000.0m +towgs84=0,0,0
# ---- To Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=lonlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
-79.8022823549574412    40.5818946695673759 869.0948340790346265

The second way step (1):

$ proj -I -v -f '%.16f' +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n +lat_1=33n +lat_2=45n +x_0=900\
000.0m +y_0=1620000.0m +towgs84=0,0,0<<EOF
2340000.0 1800000.0 EOF
#Lambert Conformal Conic
#       Conic, Sph&Ell
#       lat_1= and lat_2= or lat_0
# +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n
# +lat_1=33n +lat_2=45n +x_0=900000.0m +y_0=1620000.0m +towgs84=0,0,0
-79.8022823549574412    40.3918791792399077

The second way step (2):

$ cs2cs -v -f '%.16f' +proj=lonlat +a=6370000.0m +b=6370000.0m +towgs84=0,0,0 +no_defs +to +proj=lonlat +datum=NAD83 +no_de\ fs<<EOF
-79.8022823549574412    40.3918791792399077 EOF
# ---- From Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=lonlat +a=6370000.0m +b=6370000.0m +towgs84=0,0,0 +no_defs
# ---- To Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=lonlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
-79.8022823549574554    40.5818946695673759 869.0948340790346265

Original Question (unmodified)

I want to convert Lambert Conformal (x, y) of a spheroid datum to (lon, lat) of an ellipsoid datum (NAD83). I thought this could be done in two ways.

The first one is to use cs2cs in one step:

$ cs2cs -v -f '%.16f' +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n +lat_1=33n +lat_2=45n +x_0=900000.0m +y_0=1620000.0m +to +proj=lonlat +datum=NAD83 +no_defs<<EOF
2340000.0 1800000.0
EOF
# ---- From Coordinate System ----
#Lambert Conformal Conic
#       Conic, Sph&Ell
#       lat_1= and lat_2= or lat_0
# +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n
# +lat_1=33n +lat_2=45n +x_0=900000.0m +y_0=1620000.0m
# ---- To Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=lonlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
-79.8022823549574412    40.3918791792399077 0.0000000000000000

The second is to do in two steps. (1) to use proj to convert (x, y) to (lon, lat):

$ proj -I -v -f '%.16f' +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n +lat_1=33n +lat_2=45n +x_0=900000.0m +y_0=1620000.0m<<EOF
2340000.0 1800000.0
EOF
#Lambert Conformal Conic
#       Conic, Sph&Ell
#       lat_1= and lat_2= or lat_0
# +proj=lcc +no_defs +a=6370000.0m +b=6370000.0m +lon_0=97w +lat_0=40n
# +lat_1=33n +lat_2=45n +x_0=900000.0m +y_0=1620000.0m
-79.8022823549574412    40.3918791792399077

(2) to use cs2cs convert its datum.

$ cs2cs -v -f '%.16f' +proj=lonlat +a=6370000.0m +b=6370000.0m +no_defs +to +proj=lonlat +datum=NAD83 +no_defs<<EOF
-79.8022823549574412    40.3918791792399077
EOF
# ---- From Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=lonlat +a=6370000.0m +b=6370000.0m +no_defs
# ---- To Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=lonlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
-79.8022823549574554    40.3918791792399148 0.0000000000000000

But as you can see, the first one's results are same with the second one's first step. And, the two ways produce different answers in the end. What is wrong here?

Best Answer

The first workflow and the first step of the second workflow are doing the same thing. They're unprojecting the data from the LCC-based projected coordinate reference system (CRS) to its sphere-based geographic CRS.

...Originally I thought the second step of the second workflow was doing a datum transformation between the sphere and the NAD83 datum (based on GRS80 ellipsoid) but states that the transformation parameters are zeroes. Using zero parameters means that the sphere and the ellipsoid are assumed to have the same center point.

But the longitude and latitude are honestly, almost the same, so it's possible you're seeing 'space dust' due to internal number conversions. If a geocentric translation (3 parameter) in XYZ space had occurred even with zero parameters, I would expect to see a larger difference in the latitude value and a new ellipsoidal height value (not zero).

Using the Esri projection engine, I set up a custom ProjCRS and GeoCRS using your definition, then added a geocentric translation from the custom GeoCRS to NAD83. The results were:

Unprojection to sphere (your result): 
       -79.8022823549574412    40.3918791792399077
Transformation: 
       -79.80228235495744      40.58189466956738        869.0948340798497

You can see that there are significant differences in the latitude and height values.

Related Question