请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 42378|回复: 16

[经验总结] 小白学习cartopy画地图的第六天

[复制链接]

新浪微博达人勋

发表于 2021-3-8 19:47:14 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
小白学习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)
20210307215624924.png

等值线及标签
ct = ax.contour(olon,olat,var_new1,clevs)
clabel = ax.clabel(ct,fmt='%i')
2021-03-08_193442.png

色标
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')

2021-03-08_193610.png

图中发现两个问题:
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)])

2021-03-08_193808.png
画省界方法一
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())
2021-03-08_193908.png

画地级市也是一样的

方法三

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


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

2021-03-08_193947.png
拓展

那么问题来了,如果我想单独凸显某个地级市怎么做呢?
比如:

2021-03-08_194016.png

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文件。







地图初体验七.rar

303.11 KB, 下载次数: 81, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-3-8 20:18:07 | 显示全部楼层
赞!
python画图上手简单,就是出图分辨率让人捉急
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-3-8 20:31:46 | 显示全部楼层
小其其格 发表于 2021-3-8 20:18
赞!
python画图上手简单,就是出图分辨率让人捉急

那是不是你要求很高呀,我记得分辨率可以调的呀  比如:plt.figure(dpi=300)
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-3-8 20:59:56 | 显示全部楼层
jl2587t 发表于 2021-3-8 20:31
那是不是你要求很高呀,我记得分辨率可以调的呀  比如:plt.figure(dpi=300)

还是太低感觉……很多期刊都要矢量图,python画eps简直灾难……
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-3-9 08:53:48 | 显示全部楼层
你好,读取micaps数据画图挺不错的,请问能不能出一期读取NC数据的?感谢!特别是有很多个NC文件需要全部遍历的,一直没搞懂怎么读取。请指教!数据是FY-4A的L2级产品LMI数据,一直没搞懂!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-3-9 10:34:16 来自手机 | 显示全部楼层
VV09 发表于 2021-03-09 08:53
你好,读取micaps数据画图挺不错的,请问能不能出一期读取NC数据的?感谢!特别是有很多个NC文件需要全部遍历的,一直没搞懂怎么读取。请指教!数据是FY-4A的L2级产品LMI数据,一直没搞懂!

要不你邮箱给我发一个,我画画看
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-3-9 11:20:52 | 显示全部楼层
好的!还有就是,上面的代码,read_mdfs.py模块里 buf = open('data_table.pickle', 'rb')这个打开的文件是什么?找不到这个文件,整个代码跑不起来~~可否发我邮箱一份,感谢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-3-9 11:45:38 来自手机 | 显示全部楼层
VV09 发表于 2021-03-09 11:20
好的!还有就是,上面的代码,read_mdfs.py模块里 buf = open('data_table.pickle', 'rb')这个打开的文件是什么?找不到这个文件,整个代码跑不起来~~可否发我邮箱一份,感谢

之前没有发给过你吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-3-9 12:40:22 | 显示全部楼层
没有,'data_table.pickle'这个文件麻烦发给我一份,谢谢!
QQ:497456613@qq.com
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-3-9 13:24:00 | 显示全部楼层
VV09 发表于 2021-3-9 12:40
没有,'data_table.pickle'这个文件麻烦发给我一份,谢谢!
QQ:

你的NC文件发我一份,我读读看
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表