对 draw_map 函数的改进，以及 1.03 版本下 cartopy <= 0.19 时的报错

[ZhaJiMan]

改动主要有三个：

1. 1.03 版本里用到了pyproj 的 Transformer 类对每个点进行坐标变换，加快了运行速度。其实 shapely 的文档介绍过，利用 shapely.ops.transform 函数，将 Transformer.transform 方法作为可调用对象传入，可以将整个几何对象从源坐标系变换到新坐标系中。变换失败的数值会以 inf 填充，而不会产生报错。
2. 加入了 ax.projection 是否等价于地图数据坐标系的判断。这样一来，如果 GeoAxes 基于无参数的等经纬度投影，就能直接跳过坐标变换的步骤。
3. 考虑到 map_polygon 参数只可能是 MapPolygon（继承自 sgeom.MultiPoygon）或 sgeom.MultiLineString 对象，简化了绘图步骤。

import numpy as np
import shapely.ops as sops

def draw_map(map_polygon: Union[MapPolygon, sgeom.MultiLineString], **kwargs):
    ax = plt.gca()
    crs_from = ccrs.PlateCarree()
    crs_to = ax.projection
    if crs_from != crs_to:
        transformer = Transformer.from_crs(crs_from, crs_to, always_xy=True)
        map_polygon = sops.transform(transformer.transform, map_polygon)

    if 'color' not in kwargs and 'c' not in kwargs:
        kwargs['color'] = 'k'

    if isinstance(map_polygon, sgeom.MultiPolygon):
        for polygon in map_polygon.geoms:
            for ring in [polygon.exterior] + polygon.interiors[:]:
                coords = np.array(ring.coords)
                ax.plot(coords[:, 0], coords[:, 1], **kwargs)
    elif isinstance(map_polygon, sgeom.MultiLineString):
        for line in map_polygon.geoms:
            coords = np.array(line.coords)
            ax.plot(coords[:, 0], coords[:, 1], **kwargs)

测试代码如下，画出兰伯特投影的中国省界

from cnmaps import get_adm_maps, draw_maps
provinces = get_adm_maps(level='省')
crs = ccrs.LambertConformal(central_longitude=105, standard_parallels=(25, 47))
fig = plt.figure()
ax = fig.add_subplot(111, projection=crs)
draw_maps(provinces, lw=0.5)
ax.set_extent([78, 128, 15, 55], crs=ccrs.PlateCarree())

fig.savefig('draw_maps.png')
plt.close(fig)

结果图片为

基于同样的测试代码，对 1.03 的 draw_map 函数应用 line_profiler 模块，耗时 14 s。

Total time: 14.3249 s
File: .\revised.py
Function: draw_map at line 236

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   236                                           @profile
   237                                           def draw_map(map_polygon, **kwargs):
   238                                               """
   239                                               绘制单个地图边界线
   240                                           
   241                                               参数:
   242                                                   map_polygon (Union[MapPolygon,  sgeom.MultiLineString]): 地图边界线对象
   243                                           
   244                                               示例:
   245                                                   >>> beijing = get_adm_maps(city='北京市', level='市')[0]
   246                                                   >>> beijing
   247                                                   {'国家': '中华人民共和国', '省/直辖市': '北京市', '市': '北京市', '区/县': None, '级别': '市', '来源': '高德', '类型': '陆地', 'geometry': <cnmaps.maps.MapPolygon object at 0x7f92fe514490>}
   248                                           
   249                                                   >>> beijing_polygon = get_adm_maps(city='北京市', level='市')[0]['geometry']
   250                                                   >>> beijing_polygon
   251                                                   <cnmaps.maps.MapPolygon object at 0x7f92fe4fe410>
   252                                           
   253                                                   >>> draw_map(beijing_polygon)
   254                                           
   255                                                   >>> get_adm_maps(city='北京市', level='市')[0]['geometry'] == get_adm_maps(city='北京市', level='市', only_polygon=True, record='first')
   256                                                   True
   257                                                   >>> draw_map(get_adm_maps(city='北京市', level='市', only_polygon=True, record='first'))
   258                                               """
   259        34        746.3     22.0      0.0      ax = plt.gca()
   260        34      14779.4    434.7      0.1      crs = ccrs.PlateCarree()
   261        34     255073.7   7502.2      1.8      transformer = Transformer.from_crs(crs, ax.projection, always_xy=True)
   262      2036      51407.3     25.2      0.4      for geometry in map_polygon.geoms:
   263      2002       1967.1      1.0      0.0          if isinstance(map_polygon, sgeom.MultiPolygon):
   264      2002        877.3      0.4      0.0              try:
   265      2002      81090.5     40.5      0.6                  coords = geometry.boundary.coords
   266                                                       except NotImplementedError:
   267                                                           exterior_coords = geometry.exterior.coords
   268                                                           interiors = geometry.interiors
   269                                                           xs = []
   270                                                           ys = []
   271                                                           for coord in exterior_coords:
   272                                                               try:
   273                                                                   x, y = transformer.transform(*coord)
   274                                                               except AttributeError:
   275                                                                   x, y = coord
   276                                                               xs.append(x)
   277                                                               ys.append(y)
   278                                           
   279                                                           ax.plot(xs, ys, **kwargs)
   280                                                           for interior in interiors:
   281                                                               interior_coords = interior.coords
   282                                                               xs = []
   283                                                               ys = []
   284                                                               for coord in interior_coords:
   285                                                                   try:
   286                                                                       x, y = transformer.transform(*coord)
   287                                                                   except AttributeError:
   288                                                                       x, y = coord
   289                                                                   xs.append(x)
   290                                                                   ys.append(y)
   291                                                               ax.plot(xs, ys, **kwargs)
   292                                                           continue
   293                                                   elif isinstance(map_polygon, sgeom.MultiLineString):
   294                                                       coords = geometry.coords
   295      2002       4784.1      2.4      0.0          xs = []
   296      2002       4023.5      2.0      0.0          ys = []
   297    788498     416001.5      0.5      2.9          for coord in coords:
   298    786496     311393.0      0.4      2.2              try:
   299    786496   11158540.7     14.2     77.9                  x, y = transformer.transform(*coord)
   300                                                       except AttributeError:
   301                                                           x, y = coord
   302    786496     395319.6      0.5      2.8              xs.append(x)
   303    786496     406209.2      0.5      2.8              ys.append(y)
   304      2002       1573.8      0.8      0.0          if kwargs.get("color"):
   305      2002    1221139.4    610.0      8.5              ax.plot(xs, ys, **kwargs)
   306                                                   else:
   307                                                       ax.plot(xs, ys, color="k", **kwargs)

而改动后的 draw_map 耗时 2.8 s。

Timer unit: 1e-06 s

Total time: 2.84613 s
File: .\revised.py
Function: draw_map at line 236

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   236                                           @profile
   237                                           def draw_map(map_polygon, **kwargs):
   238        34        571.5     16.8      0.0      ax = plt.gca()
   239        34      14208.4    417.9      0.5      crs_from = ccrs.PlateCarree()
   240        34         39.5      1.2      0.0      crs_to = ax.projection
   241        34      42037.9   1236.4      1.5      if crs_from != crs_to:
   242        34     233624.2   6871.3      8.2          transformer = Transformer.from_crs(crs_from, crs_to, always_xy=True)
   243        34    1272800.0  37435.3     44.7          map_polygon = sops.transform(transformer.transform, map_polygon)
   244                                           
   245        34         68.9      2.0      0.0      if 'color' not in kwargs and 'c' not in kwargs:
   246        34         26.3      0.8      0.0          kwargs['color'] = 'k'
   247                                           
   248        34         49.3      1.5      0.0      if isinstance(map_polygon, sgeom.MultiPolygon):
   249      2036      36299.7     17.8      1.3          for polygon in map_polygon.geoms:
   250      4004      87826.7     21.9      3.1              for ring in [polygon.exterior] + polygon.interiors[:]:
   251      2002     104124.4     52.0      3.7                  coords = np.array(ring.coords)
   252      2002    1054452.6    526.7     37.0                  ax.plot(coords[:, 0], coords[:, 1], **kwargs)
   253                                               elif isinstance(map_polygon, sgeom.MultiLineString):
   254                                                   for line in map_polygon.geoms:
   255                                                       coords = np.array(line.coords)
   256                                                       ax.plot(coords[:, 0], coords[:, 1], **kwargs)

个人感觉 shapely 的 transform 函数有助于加快运行速度，同时将坐标变换的部分移出循环，从而简化代码。

另外测试发现，当 cartopy <= 0.19 时，cartopy 的 CRS 类并非简单继承自 pyproj 的 CRS 类，而是与 _crs.pxd、_crs.pyx 和 _proj4.pxd 等文件有关。具体我也不懂，反正 0.19 的 cartopy 在执行 Transformer.from_crs(crs_from, crs_to, always_xy=True) 时会产生如下错误

CRSError                                  Traceback (most recent call last)
<ipython-input-5-e063796a802e> in <module>
----> 1 t = Transformer.from_crs(crs_from, crs_to, always_xy=True)

/data/anaconda3/envs/pp/lib/python3.7/site-packages/pyproj/transformer.py in from_crs(crs_from, crs_to, skip_equivalent, always_xy)
    151         transformer = Transformer(
    152             _Transformer.from_crs(
--> 153                 CRS.from_user_input(crs_from),
    154                 CRS.from_user_input(crs_to),
    155                 skip_equivalent=skip_equivalent,

/data/anaconda3/envs/pp/lib/python3.7/site-packages/pyproj/crs.py in from_user_input(cls, value)
    427         if isinstance(value, _CRS):
    428             return value
--> 429         return cls(value)
    430 
    431     def get_geod(self):

/data/anaconda3/envs/pp/lib/python3.7/site-packages/pyproj/crs.py in __init__(self, projparams, **kwargs)
    294             projstring = projparams.to_wkt()
    295         else:
--> 296             raise CRSError("Invalid CRS input: {!r}".format(projparams))
    297 
    298         super(CRS, self).__init__(projstring)

CRSError: Invalid CRS input: <cartopy.crs.PlateCarree object at 0x7fa8fc554b30>

对于老版本的 cartopy，也许可以继续沿用 ax.projection.transform_point。

[Clarmy]

@ZhaJiMan Z 对于 cartopy <= 0.19 的问题其实我也发现了，我准备做的是在下一个版本中，在 requirements.txt 里对 cartopy 的版本限制在 0.20 及以上，这样安装的时候如果 cartopy 的版本不满足就无法安装。

[Clarmy]

@ZhaJiMan 对于函数优化的部分，如果代码已经测试过了，可以在本地跑一遍单元测试，如果单元测试没有报错，可以直接提交PR。

跑单元测试的方法：

1. 在cnmaps的环境中安装 pytest: pip install pytest
2. 进入 tests 目录： cd tests
3. 执行 pytest --verbose -p no:warnings . 即可