geopandas中一个非常重要,之前有让我非常头大的一个内容,就是坐标系。Dos中的语句很清晰,但是缺少坐标系的背景知识,看着每一个都认识的字母凑在一起完全不认识的语句,一个头两个大。让我看不懂我就更不想学下去了。所以梨准备从基础开始学起来,理清楚坐标系的门门道道,不至于混乱不清。混乱不清的逻辑什么的最令人讨厌了。
借鉴了许多前人的文档,我只是针对我的逻辑点总结了一下坐标系的基础知识,觉得写的不好的盆友们可以看看这些大大们的文档(他们说的超好的❤❤( ̄︶ ̄)↗ ):
链接: (数据科学学习手札75)基于geopandas的空间数据分析——坐标参考系篇.
链接: 地球坐标系与投影方式的理解(关于北京54,西安80,WGS84;高斯,兰勃特,墨卡托投影)
1. 坐标系基础知识
先上图,预告一下逻辑框架!!
这啥?别问,问就是参考椭球之间的区别!!
1.1 大地坐标系
- 定义:大地测量中以参考椭球面为基准面建立起来的坐标系。地面点的位置用大地经度、大地纬度和大地高度表示:(L, B, H)。
大家都知道地球它就是个椭球🌍,地球表面坑坑挖挖,高低不平,高有喜马拉雅山⛰,低有马里亚纳海沟。所以这个地球是无法用完美的数学椭球进行描述的,只能力求误差分布均衡,地势高低相补,总误差最小,尽可能拟合。
- WGS84坐标系: 椭球若以地球为椭球的球心,最大程度地拟合坑坑挖挖的整个地球🌍,就形成了WGS84坐标系。而且,这个椭球公公正正不偏不倚,没有偏向任何一个地区,当椭球体模型的几何中心与地球质心重合时,模型就会最接近整个地球,这叫做质心坐标系。
但是呢,因为它尽量逼近整个地球表面,优点是范围大,缺点是局部不够准确。你就这么一小个球就想高精度描述每一个地区??痴心妄想,痴人说梦!!━┳━ ━┳━
所以呢,每个国家根据自己的国情,又建立对本国精度更高的椭球(你国外变形,谁爱管谁管ㄟ( ▔, ▔ )ㄏ)。 - Beijing54坐标系:参考前苏联坐标系,它在前苏联区域这个曲面尽量逼近,而其它国家地区偏多少它不管。它以苏联的普尔科沃为中心,离那越远,误差就越大。
- Xian80坐标系:西安80是面向中国的,所以它在中国区域这个曲面尽量逼近,而其它国家地区偏多少它不管,而且这个逼近是以西安附近的大地原点为中心的。也就是说,在西安附近,模型和真实地表参考海平面重合,误差趋近于0,而离大地原点越远的地方,误差越大。
不管是北京54还是西安80,它们的椭球是偏向中国的,几何中心与地球质心不重合,而是参考了地球上的某一个点,所以叫做参心坐标系。
国际上多用WGS84,但是中国的还是Xian80坐标系,Beijing54被后浪冲在了沙滩上。
1.2 空间直角坐标系(暂时不重要)
- 定义:大地测量中以参考椭球中心为原点,以原点到0度经线与赤道交点的射线为x轴,原点到90度经线与赤道交点的射线为y轴,以地球旋转轴向北为z轴:(x, y, z)。
- 大地坐标系与地理坐标系都依赖与一个参考椭球。不同的是,大地坐标系以基准面为基准(标准海平面),空间直角坐标系则以一个点为基准,需要确定中心点。
只要确定了椭球基本参数,则大地坐标系和空间直角坐标系就相对确定了,只是两种不同的表达而矣,这两个坐标系的点是一一对应的。
不同椭球的基本参数,看这里👇:
链接: 地理坐标系、大地坐标系与地图投影与重投影详解
1.3 投影坐标系
地理坐标系 + 投影方法 = 投影坐标系
最常用的投影方法是横轴墨卡托, 在横轴墨卡托的基础上,形成了UTM和高斯-克吕格两种投影方法;当然,还有很多都其他的投影方法,比如👆,上面那个链接超详细:
·兰伯特Lambert投影(正轴等角割圆锥投影)
·阿伯斯Albers投影(正轴等积割圆锥投影)
·…
从3D→2D总有失真与变形的发生,不同的投影方法有不同的用途与适用范围,存在即合理。
2. Geopandas坐标系操作
背景知识讲完了,深刻体会思维导图对我个人的重要意义。无逻辑就丧失了我学下去的动力。之前会让我打退堂鼓的东西,现在学来也不是很难,真是奇了怪了。
接下来是Geopandas Document大法~~
2.1 CRS语句
任何能够被pyproj.CRS.from_user_input()语句接受的字符串,都能被geopandas正确使用。
- EPSG分区:官方推荐使用
GeoDataFrame(..., crs = 'EPSG:4326')
gdf.crs = 'EPSG:4326'
gdf.to_csr('EPSG:4326')
也就是说,推荐使用EPSG分区编号,中国国内的分区编号可以查看一下链接:
链接:Spatial Reference
- proj4语句:+号声明每一个元素,建议转换为EPSG;
‘+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0’
proj = utm: 投影方法为墨卡托
zone=11: 分区为北11区;南11区(11s)
datum=WGS84: 基准面为WGS84椭球面
units=m: 单位为米
ellps=WGS84: 如何计算地球的圆度
如果输入非得是proj4形式,修改为epsg:
import pyproj
crs = pyproj.CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")
crs.to_epsg()
2163
2.2 添加CRS
当通过geopandas.read_file()读取shp或者其他地理文件时,即使是不传入crs参数,geopandas也会自动读取文件的crs.
Geopandas.crs
Geoseris.crs
#查看geopandas/Geoseris的crs
当geopandas缺少crs信息时,可以通过set_crs()设置crs
my_geoseries = my_geoseries.set_crs("EPSG:4326")
my_geoseries = my_geoseries.set_crs(epsg=4326)
#也可以传出proj4语句(一定要是字符串,不接受dict)
但是呢,我用set_crs语句设置geoDataFrame的crs,返回值是空,也就是说没有设置上
用geoDataFrame.crs = {'init' :'epsg:4326'}就成功了
2.3 重设CRS
world = world.to_crs("EPSG:3395")
#将一个crs转换到另一个crs
#前提:要有crs,不然该语句失效报错
2.4 geopanas坐标输入顺序
对于geopanas 中的每一个坐标点,输入的形式一定要是:
(lon, lat) = (经度,纬度)
虽然在crs中,标准的坐标点的顺序应该是(lat,lon);但是,Geopandas自行处理二者之间的差距,咱们不同担心~~
就记着geopandas需要的输入形式就好啦~~