pyproj 将 UTM 转换为 lat/long 3 度

Question

我正在从形状文件中读取 UTM 点数据。与形状文件关联的 geopandas CRS 字符串是：

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]]

我使用 pyproj:

将点数据转换为 lat/long

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)

这给了我正确的纬度，但经度刚好偏出 3 度。我将 UTM 区域更改为 utm=+20，但现在向另一个方向倾斜了 3 度。我还尝试用 x_0=500000 设置东偏移量，用 lon_0=-60 设置中心经度，但这没有区别。最后，我尝试使用 CRS 字符串中的 EPSG 设置之一设置投影系统，例如

projn = pyproj.CRS.from_epsg(6326)

但是给出了错误信息 CRSError: Invalid projection: epsg:6326: (Internal Proj Error: proj_create: crs not found)。感谢任何建议，因为我是 GIS 的新手并且发现很难理解投影。问题的说明如下所示：

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("")

编辑：经过进一步调查，我发现我应该使用横向墨卡托投影，而不是通用横向墨卡托投影。如果我使用：

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')

然后点绘制在正确的位置。

Answer 1

建议：

使用Transformer (https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1)
WKT可以直接作为输入

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))

输出：

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

pyproj 将 UTM 转换为 lat/long 3 度

pyproj conversion UTM to lat/long out by 3 deg

python

gis

pyproj