Py学习  »  Python

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

气象学家 • 1 年前 • 416 次点击  

项目地址: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=(1010))
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[0for s in data]
y = [s[1for s in data]
z = [s[2for s in data]

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

fig = plt.figure(figsize=(1010))
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(601501000)
lat = np.linspace(0601000)
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(601501000)
lat = np.linspace(0601000)

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=(1010))
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)

cnmaps==1.0.1 下运行耗时达到了惊人的 902 秒(15分钟,emmmmmmm),而在 cnmaps==1.1.0 下运行耗时为 52 秒(虽然可能还是有点慢,但没有15分钟那么离谱),当然后续还会继续优化。

再次感谢炸鸡人在新版本 cnmaps 性能优化上的贡献。


声明:欢迎转载、转发本号原创内容,可留言区留言或者后台联系小编(微信:gavin7675)进行授权。气象学家公众号转载信息旨在传播交流,其内容由作者负责,不代表本号观点。文中部分图片来源于网络,如涉及作品内容、版权和其他问题,请后台联系小编处理。






往期推荐

 ERA5-Land陆面高分辨率再分析数据(~32TB)

★ ERA5常用变量再分析数据(~18TB)

 TRMM 3B42降水数据(Daily/3h)

 科研数据免费共享: GPM卫星降水数据

 气象圈子有人就有江湖,不要德不配位!

 请某气象公众号不要 “以小人之心,度君子之腹”!

 EC数据商店推出Python在线处理工具箱

★ EC打造实用气象Python工具Metview

★ 机器学习简介及在短临天气预警中的应用

★ AMS推荐|气象学家-海洋学家的Python教程

★ Nature-地球系统科学领域的深度学习及理解

★ 采用神经网络与深度学习来预报降水、温度


   欢迎加入气象学家交流群   

请备注: 姓名/昵称-单位/学校-研究方向

未备注的不通过申请



❤️ 「气象学家」 点赞

Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/146512
 
416 次点击