pyproj conversion UTM to lat/long out by 3 deg

Question

I am reading UTM point data from a shape file. The geopandas CRS string associated with the shape file is:

PROJCS["WGS 84 / Falk",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-60],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

I convert the point data to lat/long using pyproj:

projn = pyproj.Proj('+proj=utm +zone=21 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +lon_0=-60')
ylon, ylat = projn(utm_e, utm_n, inverse=True)

This gives me the correct latitude but the longitude is exactly 3 degrees out. I changed the UTM zone to utm=+20 but now am 3 degress out in the other direction. I also tried setting the false easting with x_0=500000 and the central longitude with lon_0=-60 but that made no difference. Finally, I tried setting a projection system using one of the EPSG settings in the CRS string eg

projn = pyproj.CRS.from_epsg(6326)

but that gave the error message CRSError: Invalid projection: epsg:6326: (Internal Proj Error: proj_create: crs not found). Would appreciate any suggestions as I am new to GIS and finding it difficult to understand projections. An illustration of the problem is shown below:

import pyproj

# Define points to process

well_dict = {}
well_dict['14/09-1']  = [-59.384869, -49.319319, 544706.1998681703, 4536872.299629836]
well_dict['14/09-2']  = [-59.349633, -49.247411, 547331.1995800878, 4544831.399531693]
well_dict['14/10-1']  = [-59.176736, -49.275033, 559882.9998544991, 4541663.299837932]
well_dict['14/13-1']  = [-59.496483, -49.417692, 536519.4998223917, 4525987.899903225]
well_dict['14/05-1A'] = [-59.177950, -49.162275, 559930.8995227005, 4554179.299983081]
well_dict['14/24-1']  = [-59.297611, -49.812692, 550533.2498328319, 4481958.129964017]

# Define projections

projns = {}
projns['base'] = pyproj.Proj('+proj=utm +zone=21 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +lon_0=-60')
projns['6326'] = pyproj.CRS.from_epsg(6326) #FIXME: gives CRSerror
# projns['6326'] = pyproj.Proj('epsg:6326')
# projns['7030'] = pyproj.Proj('epsg:7030')
# projns['9001'] = pyproj.Proj('epsg:9001')

# Convert UTMs for each well to lat/long using each projection

for well in well_dict:
    xlon, xlat, utm_e, utm_n = well_dict[well]
    print("%-9s  %10.2f %10.2f  %7.3f %7.3f" % (well, utm_e, utm_n, xlat, xlon), end='')    
    
    for pname in projns:
        projn = projns[pname]
        ylon, ylat = projn(utm_e, utm_n, inverse=True)
        print("  %7.3f %7.3f" % (ylat, ylon), end='')
        
    print("")

EDIT: After further investigation I found that I should have been using the Transverse Mercator projection, not Universal Transverse Mercator. If I use:

projns['tmerc'] = pyproj.Proj('+proj=tmerc +lat_0=0 +lon_0=-60 +k=0.9996 +x_0=500000 +y_0=10000000 +datum=WGS84 +units=m +no_defs')

then the points plot in the correct position.

6326 is the id of your Datum, not the CRS. Your proj string doesn't match either of the UTM zones. You'll need to go back to the data supplier and ask them for more projection details - or try opening the shapefile in QGIS and exporting it from there in 4326 — Ian Turton, Apr 22 '22 at 14:39
Hi Ian thanks. Sorry for sounding dumb but I don't understand the difference between an EPGS ID and a datum. Also don't understand why my proj string does not match either of the UTM zones. The data came from the public domain so it might be difficult to get more projection details - what extra info do I need that is not in the WKT string in my post? Thanks for you patience — PetGriffin, Apr 24 '22 at 08:24

score 1 · Accepted Answer · answered Apr 23 '22 at 01:10

Recommendations:

Use Transformer (https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1)
The WKT can be used as input directly

import pyproj


wkt = 'PROJCS["WGS 84 / Falk",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-60],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
transformer = pyproj.Transformer.from_crs(wkt, "EPSG:4326", always_xy=True)
# Define points to process

well_dict = {}
well_dict['14/09-1']  = [-59.384869, -49.319319, 544706.1998681703, 4536872.299629836]
well_dict['14/09-2']  = [-59.349633, -49.247411, 547331.1995800878, 4544831.399531693]
well_dict['14/10-1']  = [-59.176736, -49.275033, 559882.9998544991, 4541663.299837932]
well_dict['14/13-1']  = [-59.496483, -49.417692, 536519.4998223917, 4525987.899903225]
well_dict['14/05-1A'] = [-59.177950, -49.162275, 559930.8995227005, 4554179.299983081]
well_dict['14/24-1']  = [-59.297611, -49.812692, 550533.2498328319, 4481958.129964017]

# Define projections

# Convert UTMs for each well to lat/long using each projection

for well in well_dict:
    xlon, xlat, utm_e, utm_n = well_dict[well]
    print("%-9s  %10.2f %10.2f  %7.3f %7.3f" % (well, utm_e, utm_n, xlat, xlon), end='')    
    ylon, ylat = transformer.transform(utm_e, utm_n)
    print("  %7.3f %7.3f" % (ylat, ylon))

Output:

14/09-1     544706.20 4536872.30  -49.319 -59.385  -49.319 -59.385
14/09-2     547331.20 4544831.40  -49.247 -59.350  -49.247 -59.350
14/10-1     559883.00 4541663.30  -49.275 -59.177  -49.275 -59.177
14/13-1     536519.50 4525987.90  -49.418 -59.496  -49.418 -59.496
14/05-1A    559930.90 4554179.30  -49.162 -59.178  -49.162 -59.178
14/24-1     550533.25 4481958.13  -49.813 -59.298  -49.813 -59.298

Thanks @snowman2. Can I use the transform on the data read from the original shape file. I am using `cartopy` to plot the data and read the shape file using `reader=cartopy.io.shapereader.Reader(shp_file)` — PetGriffin, Apr 24 '22 at 08:33
You should be able to do it natively with geopandas if you use the `.to_crs("EPSG:4326")` — snowman2, Apr 24 '22 at 14:25

pyproj conversion UTM to lat/long out by 3 deg

1 Answers1