新功能速览 | cnmaps绘制中国地图提供一站式解决方案的Python扩展包

项目地址：https://github.com/cnmetlab/cnmaps

cnmaps经过很多小伙伴的反馈和协作开发，目前更新了一个小版本到1.1。增加了一些新的功能，同时做了一些优化，下面我们快速看一下都有哪些变化。

可以使用conda安装了

cnmaps 现在已经通过 conda-forge 发布，相比于 pip 来说，安装更为简单，只需要执行 conda install -c conda-forge cnmaps==1.1.0 即可， cartopy 等依赖会自动安装，但建议Python解释器在3.8及以上版本。

可以对散点和风矢做裁剪了

根据气象备忘录的反馈，增加了对风矢量图的裁剪功能。

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cnmaps import get_adm_maps, clip_quiver_by_map, clip_contours_by_map, draw_map
from cnmaps.sample import load_wind

lons, lats, u, v = load_wind()

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
map_polygon = get_adm_maps(country='中华人民共和国', record='first', only_polygon=True)

spd = (u ** 2 + v ** 2) ** 0.5

qv = ax.quiver(lons, lats, u, v,transform=ccrs.PlateCarree(), zorder=2)
cs = ax.contourf(lons, lats, spd, cmap=plt.cm.RdYlBu_r,
                levels=np.linspace(spd.min(), spd.max(), 50),
                transform=ccrs.PlateCarree(), zorder=1)

clip_contours_by_map(cs, map_polygon)
clip_quiver_by_map(qv, map_polygon)

draw_map(map_polygon, color='k', linewidth=1)

plt.show()

另外，也增加了裁剪散点图的功能。

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

from cnmaps import get_adm_maps, clip_scatter_by_map, draw_map
from cnmaps.sample import load_wind

lons, lats, u, v = load_wind()
spd = (u ** 2 + v ** 2) ** 0.5

data = list(zip(lons.flatten(), lats.flatten(), spd.flatten()))

x = [s[0] for s in data]
y = [s[1] for s in data]
z = [s[2] for s in data]

map_polygon = get_adm_maps(record='first', only_polygon=True)

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
scatter = ax.scatter(x, y, s=np.array(z)*10, c=z,
                     cmap=plt.cm.RdYlBu_r, transform=ccrs.PlateCarree())
clip_scatter_by_map(scatter, map_polygon)
draw_map(map_polygon, linewidth=1)
ax.set_extent(map_polygon.get_extent(buffer=1))

plt.show()

可以做栅格掩膜了

原本自己尝试实现过一套矢量地图做栅格掩膜的方案，但是效率很低。在炸鸡人的帮助下，使用他的递归方案大幅提高了效率，使掩膜方法的可用性大幅提高。炸鸡人介绍该方法的博客：https://zhajiman.github.io/post/cartopy_polygon_to_mask/，气象备忘录对该方法的介绍：Python 利用多边形生成掩膜数组

下面是cnmaps生成一套掩膜矩阵的例子。

import numpy as np
from cnmaps import get_adm_maps
import matplotlib.pyplot as plt

lon = np.linspace(60, 150, 1000)
lat = np.linspace(0, 60, 1000)
lons, lats = np.meshgrid(lon, lat)

china = get_adm_maps(level="国", record="first", only_polygon=True, wgs84=True)
china_maskout_array = china.make_mask_array(lons, lats)

plt.imshow(china_maskout_array, cmap='binary', origin='lower')

我们也可以直接对数据进行掩膜操作。

import numpy as np
from cnmaps import get_adm_maps
import matplotlib.pyplot as plt

lon = np.linspace(60, 150, 1000)
lat = np.linspace(0, 60, 1000)

lons, lats = np.meshgrid(lon, lat)
data = np.random.random(lons.shape)

china = get_adm_maps(level="国", record= "first", only_polygon=True, wgs84=True)

maskout_data = china.maskout(lons, lats, data)

plt.figure(figsize=(20,8))

plt.subplot(121)
plt.pcolormesh(lons, lats, data)
plt.title("no maskout")

plt.subplot(122)
plt.pcolormesh(lons, lats, maskout_data)
plt.title("maskout")
plt.show()

可以导出矢量文件了

cnmaps 支持将查询到的矢量边界输出为 GeoJSON 或 ESRI Shapefile 文件了。

from cnmaps import get_adm_maps

china = get_adm_maps(level="国", record= "first", only_polygon=True)

china.to_file('./china.geojson')  # 默认为 geojson 格式文件
china.to_file('./china.shp', engine='ESRI Shapefile')  # 也可以指定 shapefile 格式文件

可以按WGS84坐标导入了

由于原始数据来自于高德API，地图坐标系为火星坐标系（GCJ02），因此在新版本的 cnmaps 中，我们增加了对坐标转换的“开关”，在get_adm_maps函数中引入了 wgs84 的参数，例如：

from cnmaps import get_adm_maps

china = get_adm_maps(level="国", record= "first", only_polygon=True, wgs84=True)

当我们传入wgs84=True 时，加载的地图会自动从火星坐标系转为WGS84坐标系，若为 False 则按火星坐标系加载，默认为wgs84=True。

裁剪和绘图的效率提高了

新版本的 cnmaps 对裁剪和绘图的程序进行了性能优化（虽然可能还有很大的优化空间）

今天特意在我的Intel Mac上测试了一下以下面这段代码：

import time

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cnmaps import get_adm_maps, clip_contours_by_map, draw_map
from cnmaps.sample import load_dem

lons, lats, data = load_dem()

t0 = time.time()
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
map_polygon = get_adm_maps(country="中华人民共和国", record="first", only_polygon=True)

cs = ax.contourf(
    lons,
    lats,
    data,
    cmap=plt.cm.terrain,
    levels=np.linspace(-2800, data.max(), 10),
    transform=ccrs.PlateCarree(),
)

clip_contours_by_map(cs, map_polygon)
draw_map(map_polygon, color="k", linewidth=1)

plt.savefig('./china.png')

t1 = time.time()

time_delta = t1 - t0

print("time_delta", time_delta)