小白学习cartopy画地图的第六天

jl2587t

小白学习cartopy画地图的第六天从开始学习画地图开始，就是打算用cartopy的，但是不知不觉跑偏了，跑去了Basemap，最近感觉这样不对，因为Basemap已经停止更新了，主要还是个人觉得cartopy逼格高一点，所以今天再次用cartopy来画图。
本帖主要内容如下：
1.介绍maskout.py文件，用于白化、掩膜
2.用cartopy画一个湖南省某个时刻的24小时降水量的分部
介绍maskout.py

很多朋友可能会在一些文章里面看到以下代码：

import maskout

然后以为这是一个库，然后就去pip install maskout了，然后发现装不了，嘿嘿，本人也是这么认为的，后来发现不是。这个其实就是某位大神写的一个掩膜py文件，其代码核心在我之前的帖子里面就有，就是利用matplotlib的PathPatch来做路径边缘白化的操作，大家可以去看我之前的帖子，当然也可以去网站上找一找，有人分享了。我下载下列以后做了一些修改，因为发现写的不够普适性，容易报错，对小白不友好，有需要的也可以留言找我要。

为了文章不至于太长，这里就不贴出来了。本质就是一个py文件，不是一个复杂的库。学会了我之前的帖子，可以忽略这个东西。

我建议大家拿到maskout.py以后也不要放到库的文件夹下面，这样可以时刻提醒自己这不是库，防止换电脑，换环境以后又去pip，还找不到

画cartopy地图，降水量导入库
from read_mdfs import MDFS_Station
import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import maskout
画布和绘图区设置

plt.rcParams['font.sans-serif'] = ['SimHei']#正常显示中文
plt.rcParams['axes.unicode_minus'] = False#显示负号
fig = plt.figure(figsize=(16,9.6),facecolor='#666666', edgecolor='Blue', frameon=False)
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())


读取micaps数据
x = MDFS_Station('r20210210100000.000')
lon = x.data['Lon']
lat = x.data['Lat']
var = x.data[1011]#24小时降水量

插值成格点数据
#对经纬度进行插值，这里作者讲经纬度插值成了0.05°*0.05°的网格
olon = np.linspace(108,115,140)
olat = np.linspace(24,31,140)
olon,olat = np.meshgrid(olon,olat)#组合成一个二维的网格坐标矩阵
#站点数据插值
func = Rbf(lon,lat,var,function='cubic')
var_new1 = func(olon,olat)
#筛选数据清洗
var_new1[var_new1 < 0] = 0
var_new1[var_new1 > 116] = 116


色彩定制：自己画色标（不均匀色标）
clevs = [0.1,10.,25.,50.,100.,250.,500.]#自定义颜色列表
cdict = ['#A6F28F','#3DBA3D','#61B8FF','#0000FF','#FA00FA','#800040']#自定义颜色列表
my_cmap = colors.ListedColormap(cdict)#自定义色板
norm = colors.BoundaryNorm(clevs,my_cmap.N)#归一化




画图设置
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xticks(np.arange(108,115,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(24,31,2),crs=ccrs.PlateCarree())
ax.gridlines()

到这为止，都是之前帖子的知识点，就不做讲解了。

画等值线和填色填色
cf = ax.contourf(olon,olat,var_new1,clevs,transform=ccrs.PlateCarree(),cmap=my_cmap,norm =norm)


等值线及标签
ct = ax.contour(olon,olat,var_new1,clevs)
clabel = ax.clabel(ct,fmt='%i')


色标
position = fig.add_axes([0.74, 0.15, 0.05, 0.2])
plt.colorbar(cf, cax=position)
position.set_yticklabels((0.1,10.,25.,50.,100.,250.,500.))

掩膜
clip = maskout.shp2clip(cf,ax,'Hunan_Province')
clip = maskout.shp2clip(ct,ax,'Hunan_Province')

图中发现两个问题：
1.等值线数值标签没有白化
2.省界、地级市界没有画出

等值线数值标签修改

我的解决方法是直接修改标签坐标，定死标签位置，经过尝试，发现其实你的坐标不在的等值线上，他依旧会去找最近的等值线，很智能的一个设计。当然还有其他的方法，我只是觉得这样最简单

clabel = ax.clabel(ct,fmt='%i',manual=[(110.6,28.6),(110.6,28.3),(111.7,27.6),(112.9,26.4),(112.5,26.5),(112.7,27.3)])


画省界方法一
shpname = r'Hunan_Province.shp'
shapes = list(shapereader.Reader(shpname).geometries())
ax.add_geometries(shapes,crs=ccrs.PlateCarree(),edgecolor = 'k', facecolor = '', linewidth = 0.3)

这一种方法可以画省界，但是我的地级市SHP文件就会出错，所以我不推荐，感觉普适性有问题

方法二

直接读取shp文件用plot函数画，但是由于shp坐标点排列的问题，图片有问题

sf = Reader('Hunan_Province.shp', encoding='utf-8')
for shape_rec in sf.shapeRecords():
    pts = shape_rec.shape.points#边界点
    x,y = zip(*pts)#把经纬度分别给到x,y
    ax.plot(x, y, '-', lw=1, color='k',transform=ccrs.PlateCarree())

画地级市也是一样的

方法三

但是我发现，用同样的地图文件，用Basemap读取就没问题，所以仿照Basemap，我在maskout.py内部写了一个一样的

maskout.readshapefile('Hunan_city.shp',linewidth=1,ax=ax)


拓展

那么问题来了，如果我想单独凸显某个地级市怎么做呢？
比如：

sf = Reader('Hunan_city.shp', encoding='utf-8')

pts = sf.shapeRecords()[3].shape.points#边界点

x,y = zip(*pts)#把经纬度分别给到x,y

ax.plot(x, y, '-', lw=2, color='r',transform=ccrs.PlateCarree())

大家可以自己去理解，也可以改写mask.py，写出一个自己的库

下一个问题，如果不画单独一个省呢，比如画京津冀、珠三角？

本帖的主要知识点是介绍maskout.py，同时甩开了Basemap，新的知识点不多，注释写的少了，大家有什么不理解的可以直接留言给我
我一定给解答
记住，maskout.py是一个你自己可以修改的py文件。

小其其格

赞！
python画图上手简单，就是出图分辨率让人捉急

jl2587t

小其其格 发表于 2021-3-8 20:18
赞！
python画图上手简单，就是出图分辨率让人捉急

那是不是你要求很高呀，我记得分辨率可以调的呀  比如：plt.figure(dpi=300)

小其其格

jl2587t 发表于 2021-3-8 20:31
那是不是你要求很高呀，我记得分辨率可以调的呀  比如：plt.figure(dpi=300)

还是太低感觉……很多期刊都要矢量图，python画eps简直灾难……

VV09

你好，读取micaps数据画图挺不错的，请问能不能出一期读取NC数据的？感谢！特别是有很多个NC文件需要全部遍历的，一直没搞懂怎么读取。请指教！数据是FY-4A的L2级产品LMI数据，一直没搞懂！

jl2587t

VV09 发表于 2021-03-09 08:53
你好，读取micaps数据画图挺不错的，请问能不能出一期读取NC数据的？感谢！特别是有很多个NC文件需要全部遍历的，一直没搞懂怎么读取。请指教！数据是FY-4A的L2级产品LMI数据，一直没搞懂！

要不你邮箱给我发一个，我画画看

来自安卓APP客户端

VV09

好的！还有就是，上面的代码，read_mdfs.py模块里 buf = open('data_table.pickle', 'rb')这个打开的文件是什么？找不到这个文件，整个代码跑不起来~~可否发我邮箱一份，感谢

jl2587t

VV09 发表于 2021-03-09 11:20
好的！还有就是，上面的代码，read_mdfs.py模块里 buf = open('data_table.pickle', 'rb')这个打开的文件是什么？找不到这个文件，整个代码跑不起来~~可否发我邮箱一份，感谢

之前没有发给过你吗？

来自安卓APP客户端

VV09

没有，'data_table.pickle'这个文件麻烦发给我一份，谢谢！
QQ：497456613@qq.com

jl2587t

VV09 发表于 2021-3-9 12:40
没有，'data_table.pickle'这个文件麻烦发给我一份，谢谢！
QQ：

你的NC文件发我一份，我读读看