Python -- 地理/投影坐标系转换

当我们要在ArcGis中绘制一个图形的时候，我们需要确认三点
1、坐标系参考规则
2、要绘制的点集合
3、点的顺序
有了以上三点，我们就可以在指定的坐标系中，按照点的顺序将点连起来，绘制成一个图形。
在实际开发过程中，我们的的数据来源可能是来源于其他网站，我们所获取的数据可能不能直接满足展示需求，例如我们爬取的网站绘图的坐标系与我们使用的坐标系不同，展示出来的效果肯定是有差异的，所以在这里记录一下坐标系转换的方式以及需要用到各类资源。

Geojson格式

在网页开发过程中，当我们要在地图上展示某个区域的轮廓时，我们一般会返回给前端一段 Geojson 数据，格式如下：

{
    "type": "FeatureCollection",
    "features": [{
    		//点
            "type": "Feature",
            "geometry": {
                "type": "Point",
                "coordinates": [102.0, 0.5]
            },
            "properties": {
                "prop0": "value0"
            }
        }, {
        	//线
            "type": "Feature",
            "geometry": {
                "type": "LineString",
                "coordinates": [[102.0, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]]
            },
            "properties": {
                "prop0": "value0",
                "prop1": 0.0
            }
        }, {
        	//面
            "type": "Feature",
            "geometry": {
                "type": "Polygon",
                "coordinates": [[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]]
            },
            "properties": {
                "prop0": "value0",
                "prop1": {
                    "this": "that"
                }
            }
        }
    ]
}

其中 coordinates 属性中就是坐标信息。

坐标转换

从某网站上爬取到的地图数据是基于 Arcgis 绘图的 wkid=‘4528’ 坐标系的坐标信息
其中 wkid=‘4528’ 表示坐标系是 CGCS2000_3_Degree_GK_Zone_40 投影坐标系
然而我们需要比较通用的坐标系，一般我们需要转换为 wkid=’4326‘ 的 GCS_WGS_1984 地理坐标系

投影坐标系和地理坐标系

其实简单来讲地理坐标系就是球面坐标，参考平面是椭球面，坐标单位是经纬度；
投影坐标系是平面坐标系，参考平面是水平面，坐标单位是米、千米等。
地理坐标系WKID详情
投影坐标系WKID详情

使用CRS库转换坐标系（以CGCS2000_3_Degree_GK_Zone_40转WGS84为例）

from pyproj import CRS
from pyproj import Transformer
def marcator(x,y):
    lonlat = {}
    # CGCS2000_3_Degree_GK_Zone_40 的wkid
    crs_CGCS2000 = CRS.from_epsg(4528)
    # crs_CGCS2000 = CRS.from_wkt('PROJCS["CGCS2000_3_Degree_GK_Zone_40",GEOGCS["GCS_China_Geodetic_Coordinate_System_2000",DATUM["D_China_2000",'
    #                             'SPHEROID["CGCS2000",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],'
    #                             'PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",40500000.0],PARAMETER["False_Northing",0.0],'
    #                             'PARAMETER["Central_Meridian",120.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],'
    #                             'AUTHORITY["EPSG",4528]]')
    #我们需要转换的目标坐标系的WKID
    crs_WGS84 = CRS.from_epsg(4326)
    from_crs = crs_CGCS2000
    to_crs = crs_WGS84
    transformer = Transformer.from_crs(from_crs, to_crs)
    # 即为转换后的坐标，也可以分别使⽤数组
    new_x, new_y = transformer.transform(x, y)
    lonlat['x'] = new_x
    lonlat['y'] = new_y
    return lonlat

代码很简单，最简单的是我们知道两个坐标系都是标准坐标系，那我们只要能够获取到他们的 WKID 就可以转换了。

Attention：

CRS.from_epsg 貌似仅支持地理坐标系的 WKID

使用投影坐标系的 WKID（EPSG） 生成的 crs 属性中的 WKID 与实际使用的 WKID 不一致，有些会直接报错提示 crs not found
所以上述代码实际上运行是有问题的，因为 4528 对应的是投影坐标系
那我们可以使用 CRS.from_wkt 去定义坐标系，就像上面被屏蔽的代码，CRS.from_wkt 需要将坐标系配置文件中的内容放进去。

大地两千的坐标系配置文件可以从以下链接下载

大地两千坐标系配置文件

大地两千坐标系配置文件