"Have You Never Been Mellow" 是已故英国女爵 Olivia Newton-John 的一首歌曲，收录于她的同名专辑《Have You Never Been Mellow》中。这首歌曲于 1975 年 1 月 21 日发行。而这首歌是德国制作人 Keanu Silva 在 2020 年于 Hexagon 厂牌发布的 Future House 重制作品。

前言

昨天有同学来问问题，问题很有趣，也就正如标题所说的，如何将 WGS84 坐标系矢量图转换为国测局 02 坐标系矢量图。Jim 虽然之前没做过，但是和那个同学讨论了一番，半个小时就给搞定了 hhh。本文就来讲讲这个问题的解决思路。

解决逻辑

我们知道在 shapefile 矢量图中，子图形都是以 WKT（Well-Known Text）格式存储的，其可以表示如下几种几何对象：

其中，WKT 文本中每一个数对都是一个点，每个点都是一个二维坐标，坐标的顺序是 x，y。

为了简化问题，我们假设我们的矢量图中只有一个多边形，那么我们的问题就变成了如何将 WGS84 坐标系的多边形顶点坐标转换为国测局 02 坐标系的多边形顶点坐标。所以问题就变成了，把一个 WGS84 坐标系下的多边形拆成多个顶点，将他们的坐标进行转换，再把转换后的顶点连起来。这就是我们解决这个问题的思路。

第一步：获取多边形顶点坐标

我们知道，一个多边形的顶点坐标是以 WKT 格式存储的，所以我们只需要读取 WKT 文本，然后解析出多边形的顶点坐标即可。这里我们使用 Python 的 Geopandas 库来读取 shapefile 文件，使用 Shapely 库来解析 WKT 文本。

我们这里以一个中国分省地图举例子，名字就叫 china_province.shp：

from shapely.geometry import Polygon
import geopandas as gpd
import numpy as np

gdf = gpd.read_file('china_province.shp')

# 以第 0 行为例，第 0 行是黑龙江省
lngs, lats = gdf.iloc[0, :].geometry.exterior.coords.xy
coordinates = np.vstack([lngs, lats]).T

这样，我们就可以得到一个多边形的顶点坐标了，我们来看看输出：

array([[121.48844147,  53.33264923],
       [121.49954224,  53.33600616],
       [121.51840973,  53.33919144],
       ...,
       [121.50035095,  53.31388855],
       [121.49738312,  53.32104492],
       [121.48844147,  53.33264923]])

可以发现，每一行都是 [longitude, latitude] 的数对，并且最后一行和第一行是重复的，这是因为多边形的最后一个顶点和第一个顶点是重合的。

第二步：坐标转换

为了方便，这里我们直接使用懒人高德地图 API 进行坐标转换，这里我们使用 requests 库来进行网络请求：

import requests
import os

def translate_cord(lng: float, lon: float) -> (float, float):
    # GPS坐标转换为高德坐标
    url = 'https://restapi.amap.com/v3/assistant/coordinate/convert'
    params = {
        'key': os.getenv('AMAP_KEY'), # 这里需要你自己申请一个高德地图的 API key，这里写到了环境变量里
        'locations': f'{lng},{lon}',
        'coordsys': 'gps'
    }
    res = requests.get(url, params=params)
    lng, lon = res.json()['locations'].split(',')
    return lng, lon

这样，我们就可以将 WGS84 坐标系下的经纬度坐标转换为国测局 02 坐标系下的经纬度坐标了。另外，高德地图的用量限制是免费版 5000 次/天，如果用量大得考虑送钱了。

第三步：将转换后的坐标连起来

昨天看到一个很有帮助的文章：Create Polygon from Points，看起来是对一个列表或者 numpy.ndarray 套上一个 Polygon() 就可以了，我们也来试试：

from shapely.geometry import Polygon

l_translated = []
for lng, lat in tqdm(coordinates):
    lng, lat = translate_cord(lng, lat)
    l_translated.append([lng, lat])
translated = np.array(l_translated)

heilongjiang = Polygon(translated)

这样，变量 heilongjiang 就是一个 Polygon 对象了，然后我们把它存储成一个 shapefile 文件：

gdf_hlj = gpd.GeoDataFrame({'geometry': [heilongjiang]})
gdf_hlj.to_file('heilongjiang.shp')

最后：效果

我们来到 QGIS 里面看看效果（投影：EPSG 4500: CGCS2000 / Gauss-Kruger zone 22）：

乍一看还差不多，我们拉得近一点：

在五常市这里，我们就能看到明显的偏移了。所以在 1km 级分辨率情况下，我们这个工作就很有必要进行了。

总结

本文的工作看似很简单，但其实是一个开创性的工作（至少 Jim 和同学都没有找到网上现成的基于国测局 02 坐标的中国地图），希望大家能用得上。