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.
Advertisement
Answer
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