专栏首页气象学家Python气象绘图教程(十五)—Cartopy_5

Python气象绘图教程(十五)—Cartopy_5

本节提要:仿制中央气象台气象服务图片、关于cartopy里的投影与转换、cartopy中extent与boundary。

一、仿制中央气象台图片

从鄙人高三填报了南信的志愿开始,就一直持续的关注中央气象台,也算是一个老看客了。期间还经历了比较大型的改版(我感觉改版后变丑了O(∩_∩)O哈哈~),我觉得央台的气象服务图片做得比较好,基本上以填色图和散点图为主,辅以少量雷达图,绝大多数图片都能很好的向各层次的读者展现天气预报与实况。

今天就接着上次谈论的图例和色条,来谈谈py下仿制央台图片流程。

我仿制的图片如下:

A、仿制的高温图

恩施想要这样的高温还是比较困难的(恩施算得上矮高原了,鄂西凉都),所以修改成了地质灾害的预报。

首先涉及到资料的问题,地质灾害不在常规预报里,但是气象局也必须发这项预警。这里只能用事先做好的实验数据,预报在清江两岸有一定的山洪泥石流风险。使用的仍然是不均匀的站点资料,所以先要将站点资料网格化,变为格点资料后再用等值线填色的方式画出危险区域。最后,通过前面几节提到的添加图例的方法,完善图形。关键代码如下:

olon=np.linspace(108,111,90)
olat=np.linspace(29,32,90)
olon,olat=np.meshgrid(olon,olat)
func=Rbf(lon,lat,danger,function='linear')
danger_new=func(olon,olat)
colorlevel=[0,60,80,100]#危险等级
colordict=['white','orange','red']#颜色列表
danger_map=mcolors.ListedColormap(colordict)#产生颜色映射
norm=mcolors.BoundaryNorm(colorlevel,danger_map.N)#生成索引
ax.contourf(olon,olat,danger_new,levels=colorlevel,cmap=danger_map,norm=norm)#填色图
clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000)#白化
ax.set_title('地质灾害风险(sample data)',fontsize=12)
danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")
low_danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="orange")
labels = ['泥石流高风险','内涝积水风险']
ax.legend([danger, low_danger], labels,fontsize=7,frameon=False,title='图例',
              loc='lower left', bbox_to_anchor=(-0.01, 0), fancybox=False)

B、降水量图的仿制

使用的是最开始实验自定义colorbar时的那张图的数据,但是当时用的是色条来表示降水量,这次我们用图例的方式表示降水量,前面的步骤和A中的类似。最后添加图例时候,我用的是逐个添加,实际上可以用循环的方式添加。代码如下:

filename=r'C:\Users\lenovo\Desktop\ex1.xlsx'
df=pd.read_excel(filename)
lon=df['lon']#经度
lat=df['lat']#纬度
olon=np.linspace(108,111,90)
olat=np.linspace(29,32,90)
olon,olat=np.meshgrid(olon,olat)
rain=df['rain_new']
func=Rbf(lon,lat,rain,function='linear')
rain_new=func(olon,olat)
colorlevel=[0.1,10.0,25.0,50.0,100.0,250.0,500.0]#雨量等级
colordict=['#A6F28F','#3DBA3D','#61BBFF','#0000FF','#FA00FA','#800040']#颜色列表
rain_map=mcolors.ListedColormap(colordict)#产生颜色映射
norm=mcolors.BoundaryNorm(colorlevel,rain_map.N)#生成索引
cs= ax.contourf(olon,olat,rain_new,levels=colorlevel,cmap=rain_map,norm=norm)
clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000)#白化
############添加图例#######################
larger1 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#A6F28F")
larger2 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#3DBA3D")
larger3 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#61BBFF")
larger4 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#0000FF")
larger5 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#FA00FA")
labels = ['0~10','10~25','25~50','50~100','100~250']
ax.legend([larger1,larger2,larger3,larger4,larger5], labels,fontsize=12,frameon=False,title='图例(mm)',
              loc='lower left', bbox_to_anchor=(-0.01, -0.03), fancybox=False)

其中,mpatches.Rectangle这个语句就是添加图例上的色条,rectangle就是矩形的意思,延伸的,你也可以使用其他形状添加图例(包括自定义形状)。如何提供更丰富的标记,只能烦请各位在官网文档上查找了。

二、Cartopy里的投影与转换

我在两个月前经常碰到这个问题,有两个关于投影的—crs、transform。前一个设定投影方式,后一个涉及投影与数据转换。

首先说crs,这个是GeoAxes的基础,只有在projection设置了投影之后,才能添加地图。如果没有一步转化,有一些专属于GeoAxes的命令就无法生效:

 'AxesSubplot' object has no attribute 'add_feature'

比较典型的如add_feature、add_geometries添加地图和地理特征的命令不能使用。

crs还控制着数据绘图与边界的裁剪,比如set_extent(crs=ccrs.PlateCarree()) ,就使裁剪的方式按照PlateCarree()的方式进行边界的裁剪,一个经典的案例即兰勃脱下的使用extent进行裁剪,那么将会在地图中出现范围之外的地图:

如果主地图投影为默认投影方式(PlateCarree),那么,你在后续的绘制中是不需要进行投影转换的。

但是,在主地图不是默认投影方式时,需要进行转换。比如:

ax=fig.add_subplot(1,1,1,projection=ccrs.LambertConformal(central_longitude=110, central_latitude=35))
ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r')

这种情况下,将无法绘制出数据:

这是因为,数据默认的格式为PlateCarree下的投影,需要transform的介入:

ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r',transform=ccrs.PlateCarree())

transform=PlateCarree()的意思是,这个数据原来是PlateCarree()投影下的,转化为主图的格式后(即LambertConformal),再绘制等值线。

主图ax的投影方式

数据的投影方式

转换命令

LambertConformal

PlateCarree

transform=PlateCarree

LambertConformal

Mercator

transform=Mercator

PlateCarree

Mercator

transform=Mercator

三、通过set_boundary和set_extent绘制扇形图

cartopy官网实例有一个Custom Boundary Shape,通过该例子可以裁剪边框形状。现在尝试裁剪东亚地区的兰勃脱投影。

import numpy as np
from matplotlib.path import Path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat
vertices = [(65, 10), (65, 60), (145, 60), (145, 10), (65, 10)]#五个点,但是首尾是一样的,以连接为封闭的四边形
boundary = Path(vertices)#边界形状
fig=plt.figure(figsize=(5,5),dpi=500)
ax=plt.axes([0,0,1,1],projection=ccrs.LambertConformal(central_longitude=110, central_latitude=35))
ax.add_feature(cfeat.OCEAN)#添加海洋
ax.add_feature(cfeat.LAND)#添加陆地
ax.set_boundary(boundary, transform=ccrs.PlateCarree())#裁剪边界
ax.set_extent([65,145,0,60])#裁剪地图范围

这时,可能因为投影转换的原因,出现了几个明显的折痕,于是,进一步提高path上点的数量:

latmin = 10
latmax = 60
lonmin = 70
lonmax = 140
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
    [(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())
ax.set_extent(extent)

此时,因为vertices中的边界点的数量增多,边界连接更加圆滑。请注意,在这种方式下裁剪了边界,一定要使用set_extent来裁剪地图大小,不然就会出现地图蜷缩在图里这种情况:

在cartopy=0.17中,不能使用draw_labels=True来为除PlateCarree、Mercator之外的投影添加经纬标签,不过据说在0.18版本中已经优化,读者可以试试。通过前面提到的数据转换,绘制了一张气温图:

ax.contourf(data_lon,data_lat,temp,levels=np.arange(-40,40),cmap='Spectral_r',transform=ccrs.PlateCarree())

本文分享自微信公众号 - 气象学家(Meteorologist2019)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-06-11

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Python气象绘图教程(十六)—Cartopy_6

    本节提要:使用cartopy进行市县的色块填色、模仿geopandas绘制颜色图

    气象学家
  • Python-matplotlib 学术散点图完善

    上期的推文Python-matplotlib 学术型散点图绘制 推出后,很多小伙伴比较喜欢

    气象学家
  • Python高效批量绘图方法

    在数值预报后处理中经常需要批量出图,而基于matplotlib的图形渲染速度较慢,而提高出图的速度通常可通过两个方面来解决:

    气象学家
  • 详解 MNIST 数据集

    MNIST 数据集已经是一个被"嚼烂"了的数据集, 很多教程都会对它"下手", 几乎成为一个 "典范". 不过有些人可能对它还不是很了解, 下面来介绍一下.

    用户1558438
  • Python-EEG工具库MNE中文教程(10)-信号空间投影SSP数学原理

    projector(投影)(简称proj),也称为信号空间投影(SSP),定义了应用于空间上的EEG或MEG数据的线性操作。

    脑机接口社区
  • 用NW.js构建跨平台桌面应用(1)-入门案例

    NW.js 基于 Chromium 和 Node.js,从而可以在桌面app中使用浏览器开发技术并直接调用 Node.js 资源,甚至将一个web应用打包到本地...

    江米小枣
  • PIKOCUBE:带 LED、陀螺仪,WiFi 控制的可编程骰子

    今天给大家带来一个非常好玩的项目,带有 54 颗 LED、陀螺仪,支持 WiFi 控制的可编程骰子。

    MCU起航
  • 数据运营实战(三):用数据说话,从埋点开始

    埋点是 App 数据运营中很重要的一个环节。之前我们讨论过用户分群的方式、漏斗转化的改进,但所有 App 数据的来源是数据采集,很多时候就是 App 的埋点。

    腾讯大数据
  • 如何在 Windows 10 中移除 Internet Explorer 浏览器 如何通过控制面板删除 Internet Explorer 浏览器通过 PowerShell 删

    现在 Internet Explorer (IE)已经过时了,可以通过控制面板移除这个古老但是依然是一个伟大的浏览器

    林德熙
  • 微信、公网播放摄像机视频方案

    之前一直都是 描述easynvr的各种功能、应用场景。也介绍了几种可以将easynvr接入公网的方案。

    EasyNVR

扫码关注云+社区

领取腾讯云代金券