In [2]:
from pyproj import Proj
使用关键词参数定义投影：
In [3]:
proj1 = Proj(proj='utm',zone=10,ellps='WGS84')
使用Proj字符串定义投影：
In [4]:
proj2 = Proj('+proj=utm +zone=10 +ellps=WGS84')
使用 EPSG 编码定义投影：
注意以为是弃用的方式：
In [5]:
proj3 = Proj("+init=epsg:32667")
/opt/conda/lib/python3.12/site-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
In [6]:
proj4 = Proj("epsg:32667",preserve_units=True)
转换不同的格式。
In [7]:
proj4.definition_string()
Out[7]:
'proj=tmerc lat_0=0 lon_0=-81 k=0.9996 x_0=500000.001016002 y_0=0 datum=WGS84 units=us-ft no_defs ellps=WGS84 towgs84=0,0,0'
In [8]:
proj4.srs
Out[8]:
'+proj=tmerc +lat_0=0 +lon_0=-81 +k=0.9996 +x_0=500000.001016002 +y_0=0 +datum=WGS84 +units=us-ft +no_defs'
In [9]:
proj4.to_wkt()
Out[9]:
'CONVERSION["PROJ-based coordinate operation",METHOD["PROJ-based operation method: +proj=tmerc +lat_0=0 +lon_0=-81 +k=0.9996 +x_0=500000.001016002 +y_0=0 +datum=WGS84 +units=us-ft +no_defs"]]'
In [10]:
proj4.to_proj4()
Out[10]:
'+proj=tmerc +lat_0=0 +lon_0=-81 +k=0.9996 +x_0=500000.001016002 +y_0=0 +datum=WGS84 +units=us-ft +no_defs'
In [11]:
proj4.to_json_dict()
Out[11]:
{'$schema': 'https://proj.org/schemas/v0.7/projjson.schema.json',
'type': 'Conversion',
'name': 'PROJ-based coordinate operation',
'method': {'name': 'PROJ-based operation method: +proj=tmerc +lat_0=0 +lon_0=-81 +k=0.9996 +x_0=500000.001016002 +y_0=0 +datum=WGS84 +units=us-ft +no_defs'},
'parameters': []}
In [12]:
p=Proj('+proj=aea +lon_0=105 +lat_1=25 +lat_2=47 +ellps=krass')
x,y=p(105,36)
print('%.3f,%.3f' %(x,y))
0.000,3847866.973
Proj4 投影控制参数必须在字典 projparams 或关键字参数给定。
查看页面（ https://proj4.org/usage/projections.html ）了解有关指定投影参数的更多信息。
In [13]:
lon,lat=p(x,y,inverse=True)
print('%.3f,%.3f' %(lon,lat))
105.000,36.000
In [14]:
import math
使用 math 模块将角度转换为弧度作为参数传递，使用关键字 radians=True 来声明：
In [15]:
x,y=p(math.radians(105),math.radians(36),radians=True)
print( '%.3f,%.3f' %(x,y) )
0.000,3847866.973
In [16]:
lons=(105,106,104)
lats=(36,35,34)
x,y=p(lons,lats) #将经纬度放入元组中
print('%.3f,%.3f,%.3f' %x) #普通打印
0.000,89660.498,-90797.784
In [17]:
print('%.3f,%.3f,%.3f' %y)
3847866.973,3735328.476,3622421.811
In [18]:
type(x) #输出的类型为元组
Out[18]:
tuple
In [19]:
zip(x,y) #用 zip函数包装
Out[19]:
<zip at 0x7f902578cd40>
In [20]:
utm=Proj(proj='utm',zone=48,ellps='WGS84')
x,y=utm(105,36) #用关键字定义一个投影。
x,y
Out[20]:
(499999.99999999773, 3983948.4533356656)
投影转换
transform() 函数是在pyproj 库下， 可以进行两个不同投影的转换。
相当于 proj 程序里的 cs2cs 子程序。用法如下：
transform(p1, p2, x, y, z=None, radians=False)
x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False)这里p1与p2都是Proj对象。
在 p1，p2两个投影之间进行投影转换。
把 p1坐标系下的点(x1,y1,z1)转换到 p2所定义的投影中去。
z1是可选的，如果没有设 z1的话，将假定为 0， 并返回 x2，y2。
使用这个函数的时候要注意不要进行基准面的变换（datum）。
关键字radians只有在 p1，p2中有一个为地理坐标系时才起作用，
并且是把地理坐标系的投影的值当作弧度。 判断是否为地理坐标可以用
p1.is_latlong() 和 p2.is_latlong()函数。输入的
(x,y,z)可以分别是数组或序列的某一种形式。
例如：
In [21]:
from pyproj import Proj
albers=Proj('+proj=aea +lon_0=105 +lat_1=25 +lat_2=47 +ellps=krass')
utm=Proj(proj='utm',zone=48,ellps='krass')
albers_x,albers_y=albers(105,36)
albers_x,albers_y
Out[21]:
(0.0, 3847866.972516728)
In [22]:
utm_x,utm_y=utm(105,36)
print(utm_x,utm_y )
499999.99999999773 3984019.058813517
下边直接从 albers转为 utm坐标
In [23]:
from pyproj import transform
to_utm_x,to_utm_y = transform(albers,utm,albers_x ,albers_y)
print(to_utm_x,to_utm_y )
499999.99999999773 3984019.058813516
/tmp/ipykernel_2392/1887088792.py:2: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 to_utm_x,to_utm_y = transform(albers,utm,albers_x ,albers_y)
Geod类
Geod类主要用来求算地球大圆两点间的距离和相对方位，以及相反的操作。 同时也可以在两点间插入等分点。
该类主要包括三个函数：
正转换
fwd()函数可以进行正转换，返回经纬度。也可以用 NumPy
与常规的Python数组对象，Python序列和标量。如果
弧度为真，把输入的经纬度单位当作弧度，而不是度。距离单位为米。
fwd(self,lons,lats,az,dist,radians=False)参数为经度，纬度，相对方位，距离。
逆变换
inv()函数可以进行逆变换，已知两点经纬度，返回其前方位解，后方位角，以及距离。
inv(self, lons1, lats1, lons2, lats2, radians=False)弧线等分
npts()函数可以进行弧线等分。给出一个始点（lon1，lat1）和终点(lon2，lat2)以及等分点数目
npts。
npts(self, lon1, lat1, lon2, lat2, npts, radians=False)实例：进行两个点的等分
例如：
In [24]:
from pyproj import Geod
g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
specify the lat/lons of Boston and Portland.
In [25]:
boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
find ten equally spaced points between Boston and Portland.
In [26]:
lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
for lon,lat in lonlats: print('%6.3f %7.3f' % (lat, lon))
43.528 -75.414 44.637 -79.883 45.565 -84.512 46.299 -89.279 46.830 -94.156 47.149 -99.112 47.251 -104.106 47.136 -109.100 46.805 -114.051 46.262 -118.924
In [27]:
from pyproj import Geod
g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
london_lat = 51.+(32./60.); london_lon = -(5./60.)
compute forward and back azimuths, plus distance between Boston and Portland.
In [28]:
az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print ("%7.3f %6.3f %12.3f" % (az12,az21,dist))
-66.531 75.654 4164192.708
compute latitude, longitude and back azimuth of Portland, given Boston lat/lon, forward azimuth and distance to Portland.
In [29]:
endlon, endlat, backaz = g.fwd(boston_lon,boston_lat, az12, dist)
print("%6.3f %6.3f %13.3f" % (endlat,endlon,backaz))
45.517 -123.683 75.654
compute the azimuths, distances from New York to several cities (pass a list)
In [30]:
lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
lons2 = [boston_lon, portland_lon, london_lon]
lats2 = [boston_lat, portland_lat, london_lat]
az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
for faz,baz,d in zip(az12,az21,dist): print("%7.3f %7.3f %9.3f" % (faz,baz,d))
54.663 -123.448 288303.720 -65.463 79.342 4013037.318 51.254 -71.576 5579916.651