Python完美白化

平流层的萝卜 · 发表于 2016-3-4 20:29:02

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

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由平流层的萝卜于 2022-11-16 21:24 编辑

Python白化似乎是一个不那么容易的问题，不过这个问题相信现在已经得到了完美的解决。借助Basemap和pyshp这两个模块，实现对任意shp文件的区域进行白化，拒绝空白锯齿，甩开越界填图，让maskout后的数据完美吻合shp边界。

登录/注册后可看大图

具体的思路是，读取气象场数据，照常画图，得到matplotlib plot的一个instance，对应未白化的整图。然后进行关键的一步，读取shapefile（shp文件，比如meteoinfo的地图库文件country1.shp，是世界各国地图），提取shapefile里的感兴趣区（比如中国），形成一个clip对象，最后得到剪切图。

登录/注册后可看大图

效果图如下：一张海平面气压场的图，把中国外的地区maskout了。

放大之后来看，没有锯齿，完美吻合。

当然，也可以把县级地图叠加到上面去，只要保留clip对象不要动，同时叠加更加精细的县级底图就可以了。

第一张图的脚本如下（其他图大同小异）：

#coding=utf-8
import maskout
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import numpy
import netCDF4 as nc
font = FontProperties(fname=r"C:\\Windows\\Fonts\\impact.ttf")
font1 = FontProperties(fname=u"C:\\Windows\\Fonts\\simkai.ttf")
fig = plt.figure(figsize=(16,9.6))
ax = fig.add_subplot(111)
# ncdata=nc.Dataset(r'2004_500_grtuvw.nc')
ncdata=nc.Dataset('pres.mon.ltm.nc')
# data=ncdata.variables['t'][0,:,:]
data=ncdata.variables['pres'][0,:,:]
# exit()
# lat=ncdata.variables['latitude'][:]
# lon=ncdata.variables['longitude'][:]
lat=ncdata.variables['lat'][:]
lon=ncdata.variables['lon'][:]
lon1,lon2=lon[0],lon[-1]
lat1,lat2=lat[-1],lat[0]
nx=data.shape[1];ny=data.shape[0]
m = Basemap(llcrnrlon=lon1,llcrnrlat=lat1,urcrnrlon=lon2,urcrnrlat=lat2,projection = 'cyl')
xx,yy=m.makegrid(nx,ny)
yy=yy[::-1]
m.readshapefile('country1','whatevername',color='gray')
minval,maxval=int(numpy.amin(data)),int(numpy.amax(data))+1
cs = m.contourf(xx,yy,data,range(minval,maxval),cmap= plt.cm.get_cmap('jet') )
bar=m.colorbar(cs)
bar.set_ticks(range(minval-1,maxval,40))
bar.set_ticklabels(range(minval-1,maxval,40))
clip=maskout.shp2clip(cs,ax,'country1','China')
plt.title(u'Python Super Mask',fontproperties=font,fontsize=40)
lon1,lon2=ax.set_xlim(70,140);lat1,lat2=ax.set_ylim(15,55)
#signature=u"by 平流层的萝卜"
#plt.text(lon2-(lon2-lon1)*3.0/10,lat2+(lat2-lat1)*0.1/10,signature,fontproperties=font1,fontsize=25)
# plt.show()
plt.savefig('111.png')

复制代码

对于上面的clip=maskout.shp2clip(cs,ax,'country1','China') ，其中maskout是自己定义的一个模块，实现的功能就是最关键的对整图进行剪切，因为模块可以移植到其他程序中，方便以后使用。
把第一张图的测试数据、地图和脚本一块打包了，包括自定义maskout模块包，欢迎感兴趣的朋友下载：

Python完美白化_测试数据及脚本.rar (1.41 MB, 下载次数: 4552)

2016-3-4 20:24 上传

点击文件名下载附件

登录/注册后可看大图

shp的记录不尽相同，但只要找到你感兴趣区对应的特别标识符，然后用pyshp的shapefile模块提取即可。对于pyshp的安装，直接命令行pip install pyshp即可。
比如把除了辽宁之外的地区maskout了。

这里需要用到bou2_4p这个shp，它包含中国各省地图，通过shp中的BOU2_4M的数字找到辽宁。
QQ截图20160304190534.png

表格是从meteoinfo中打开，获取图层的Attribute Layer,找到辽宁的id，在各列中，只有BOU2_4M_这一列和省份一一对应，（其他“辽宁省”不是辽宁省主体，第一个出现的应该是本省主体）。

登录/注册后可看大图

file:///C:/Users/Administrator/AppData/Local/YNote/data/qq816ADF5E3ABFA94344F59384C4F422EE/046034d9e41f424d86f9584196d42bb1/clipboard.png
所以，只要你有感兴趣区的shp文件，那么对任意国家、任意地区，都可以Python进行完美的白化。
如有改进，实时更新。欢迎交流。

登录/注册后可看大图

2016-3-11更新：
修改几行代码,只显示几个地区/省份，其他地区/省份被白化，且大大小小所有岛屿全部画上数据。
之前看过一篇气象公众号的文章，也是关于python白化的，讲到set_clip_path不能同时遮掉两个不相连多边形之外的部分，如果同时遮蔽两个地区，需要使用多个绘图对象，即用到多个matplotlib instance。事实上，这和几个多边形无关，也无需用多个绘图对象。之所以有这样的错觉，是因为没注意到path中的codes这个命令。
如果我们打开path的函数，可以看到codes的用法：它是一套动作命令的集合，它的数量是和多边形的点的数量是一致的。codes的元素=1时，从对应的点开始起画线，这个点也是这条线的第一个点，即MOVETO，codes的元素=2时，从上一个点到这个点画线，即LINETO，codes的元素=79时，这条线达到了终点，不再画了，即CLOSEPOLY。详见下图Path中codes的用法截图（大家可以自行help（path)查看）：
QQ截图20160311160111.png

也就是说如果从一个多边形对应的岛屿A跳到另一个多边形对应的岛屿B时，只需要把A的最后一个点codes元素设成79，B的第一个点的codes元素设成1，就可以了。这样，就把所有大大小小的多边形用一系列动作串联起来，画完岛屿A画岛屿B，画完B画C。。。。直到把所有的岛屿同陆地包含起来，形成这样的一个clip对象，再通过set_clip_path剪裁，就得到了不漏一岛的白化图。
所以，多边形的多少不是限制set_clip_path的原因，而在于巧用codes。否则，多个绘图对象也会造成内存的消耗和运行速率的减慢。

登录/注册后可看大图

那么回过头来看，如果想要同时保留山东和辽宁的数据（这回，我们把两省的大大小小所有岛屿全包含进来），其他地区不要，该怎么做呢？
还是拿bou2_4p这张地图举例。
先打开meteoinfo（也可以在python里直接运行查看，但为了更方便），查看bou2_4p的属性，如图：
QQ截图20160311173326.png

可以发现，山东对应的ADCODE99是370000，辽宁对应的是210000。
那我们就需要对maskout.py进行修改，直接跳到“这块是已经修改的地方”

#coding=utf-8
###################################################################################################################################
#####This module enables you to maskout the unneccessary data outside the interest region on a matplotlib-plotted output instance
####################in an effecient way,You can use this script for free ########################################################
###############################by 平流层的萝卜####################################################################################
#####USAGE: INPUT include 'originfig':the matplotlib instance##
# 'ax': the Axes instance
# 'shapefile': the shape file used for generating a basemap A
# 'region':the name of a region of on the basemap A,outside the region the data is to be maskout
# OUTPUT is 'clip' :the the masked-out or clipped matplotlib instance.
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def shp2clip(originfig,ax,shpfile,region):
sf = shapefile.Reader(shpfile)
vertices = [] ######这块是已经修改的地方
codes = [] ######这块是已经修改的地方
for shape_rec in sf.shapeRecords():
# if shape_rec.record[3] == region: ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。
if shape_rec.record[4] in region: ######这块是已经修改的地方
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)
for contour in originfig.collections:
contour.set_clip_path(clip)
return clip

复制代码

而主程序就变成了：（直接跳到“这块就是已经修改的地方”）

#coding=utf-8
import maskout
# from matplotlib.path import Path
# from matplotlib.patches import PathPatch
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.font_manager import FontProperties
import numpy
import netCDF4 as nc
# import matplotlib
# import shapefile
font = FontProperties(fname=r"C:\Windows\Fonts\simkai.ttf")
# font1 = FontProperties(fname=u"C:\\Windows\\Fonts\\华文行楷.ttf")
fig = plt.figure(figsize=(17,8.85),dpi=80)# figsize=(16,9.6)
ax=fig.add_axes([.2,.2,.6,.68]) #[*left*, *bottom*, *width*,*height*]
def makedegreelabel(degreelist):
labels=[str(x)+u'°E' for x in degreelist]
return labels
ncdata=nc.Dataset('pres.mon.ltm.nc')
data=ncdata.variables['pres'][0,:,:]
lat=ncdata.variables['lat'][:]
lon=ncdata.variables['lon'][:]
lon1,lon2=lon[0],lon[-1]
lat1,lat2=lat[-1],lat[0]
nx=data.shape[1];ny=data.shape[0]
m = Basemap(llcrnrlon=lon1,llcrnrlat=lat1,urcrnrlon=lon2,urcrnrlat=lat2,projection = 'cyl')
xx,yy=m.makegrid(nx,ny)
yy=yy[::-1]
# m.readshapefile('country1','ss',color='gray')
m.readshapefile('bou2_4l','whatevername',color='k',linewidth=1)
minval,maxval=int(numpy.amin(data)),int(numpy.amax(data))+1
cs = ax.contourf(xx,yy,data,range(500,1040),cmap= plt.cm.get_cmap('rainbow') )
clip=maskout.shp2clip(cs,ax,'bou2_4p',[210000,370000]) #########这块就是已经修改的地方
plt.show()
exit()

复制代码

那么得到的图像就是这样的：

放大看，是这样：

再放大看，是这样：

所有岛屿大大小小都是有数据，有颜色的。
如果看不清，再看一张大连周边：

所以，哪怕岛屿的数量再多，岛屿的面积再小，也一个都不放过，真正做到寸土必争。如果需要三个或更多的省份的白化，原理完全一样，只要在主程序的clip=maskout.shp2clip(cs,ax,'bou2_4p',[210000,370000])中，把[210000,370000]里边的数字替换成其他的就好了，其他市、县、区域、国家等等，也都是一样的。尽量做到白化的简洁高效。
至于maskout.py，是否可以不用修改，一劳永逸地使用。我想，这个好像不容易做到，因为还是那句话，shp的记录各不相同，怎么选取标识符，需要因地制宜。

本次更新完毕，以后继续。

登录/注册后可看大图

2016.3.23 更新

改进了maskout.py文件，增加了一个输入参数，使maskout适用于多种投影方式，比如lambert,stereographic等。（对于只使用等经纬度投影的朋友，旧的maskout.py足够使用。）
先上图：
等经纬度投影（上），兰伯托等角投影(Lambert Conformal)（中），球面投影(Stereographic)（下）

兰伯托等角投影(Lambert Conformal)

球面投影(Stereographic)

maskout.py里，修改的地方如下图红框里。
QQ截图20160323202827.png

主程序调用变成了这样，增加了一个参数m，就是basemap的一个实例，instance.
QQ截图20160323202631.png

新的maskou.py附件在这里：

maskout.py (2.25 KB, 下载次数: 1201)

2016-3-23 20:30 上传

点击文件名下载附件，

感谢@lovechang1314的建议。

以后有时间继续跟进。

登录/注册后可看大图

北京时间2019.3.22凌晨，开贴三年后再次更新
@晋陵小生在此基础上对等值线上的label进行了进一步的优化。用我最初的脚本是没有办法去掉区域外等值线上的label的，这会限制maskout函数对contour功能的使用效果，如下：

登录/注册后可看大图

经过@晋陵小生的优化，效果如下：

登录/注册后可看大图

在此，分享该帖的地址：http://bbs.06climate.com/forum.p ... 30&page=1#pid997832。感谢@晋陵小生的优化。

登录/注册后可看大图

2020.4.23，统一回复一个常问的问题：‘local variable 'clip' referenced before assignment’
我翻了一些回复，有人在使用过程中，出现这样的提示错误：local variable 'clip' referenced before assignment。当时由于工作等原因，没有及时回复大家的这个疑问。现在统一回复如下：出现这样的错误，是因为在maskout.py里，if shape_rec.record[n] in region:(n是一个数字) 判断在循环内都是否，最终导致，maskout.py并没有找到符合条件的、用于剪裁的多边形，即maskout.py里后续的代码中的clip没有被赋值，就直接被contour.set_clip_path(clip)引用。
导致这个错误的深层原因，就是if shape_rec.record[n] in region这个判断语句没有使用好。举个例子，用python打印某条shape_rec.record，结果是['CH', 'CHN', 'China', 'China', 1281008318.0, 9600000.0, 3605750.0, 'Renminbi Yuan', 'CNY', 'N', '2']，那么明显可以看到第三列和第四列是中国（其实通过打印更多的shape_rec.record，可知第四列才是中国，第三列是辖区，比如对于UK这种，第四列是united kingdom，第三列可以是british columbia之类的各种海外岛屿）。所以正确的做法是，选择第四列作为标识符，即if shape_rec.record[n] in region里的n应该是3，对应主程序里的clip=maskout.shp2clip(cf,ax,m,'country1',['China'])；如果随便瞎选一个，让n=7，对应的是货币，并非国别，所以无论怎么循环，这一列也没有=‘China’，所以就导致了clip最终为空，出现local variable 'clip' referenced before assignment这个错误提示。

因此，在使用maskout.py的时候，必需事先了解自己的shp文件的规则信息，有哪列信息可以作为判断的标识，来提取感兴趣区，这是前提工作。

lovechang1314 · 发表于 2016-3-22 17:02:52

本帖最后由 lovechang1314 于 2016-3-22 17:06 编辑

简单测试了以下，可以用。
有个问题，projection='cyl' 应该是不需要转换坐标的吧？
如果用'stere'这个时候画图的时候需要转换坐标，这个时候就没办法clip，可否把坐标转换的代码加入到你的脚本之中？
比如
x= np.arange(70,140.1, 0.5)
y= np.arange(20,70.1,0.5)
lon,lat=np.meshgrid(x,y)
inlon, inlat = m(lon, lat)
那么在maskout.py中读取shp文件的点时就进行转换在append，这样应该就可以了吧？
如图：

C:\Users\Administrator\Desktop\basemap.png

平流层的萝卜 · 发表于 2016-3-23 10:18:21

lovechang1314 发表于 2016-3-22 17:02
简单测试了以下，可以用。
有个问题，projection='cyl' 应该是不需要转换坐标的吧？
如果用'stere'这个时 ...

非常好的建议。当时没有考虑到投影的问题，我再测试一下，改一改。

wanghp · 发表于 2016-5-4 21:42:41

关于只画海上的学习，参照楼主的教程，选用高精度的海洋数据（10m_ocean.shp）

非对称 · 发表于 2017-5-30 08:31:46

平流层的萝卜发表于 2017-5-29 19:08
哦？愿闻其详啊

@staticmethod
def MergePaths(paths):
# paths为path对象列表数组
path = Path.make_compound_path(*paths)
return path

紫依 · 发表于 2020-8-31 16:28:58

UnicodeDecodeError: 'utf-8' codec can't decode byte 0xba in position 0: invalid start byte
请问总是报这样的错是什么原因？

qq492947833 · 发表于 2020-10-19 12:37:58

平流层的萝卜发表于 2020-10-18 09:44
clip没有被赋值，请参见我帖子内容的结尾部分

我已经找到原因了楼主，在您最初的maskout里面用的命令是“if shape_rec.record[7] == region:”，这样就只能mask一个省，如果把命令改成“if shape_rec.record[7] in region:”就可以使用list，mask多个省了

kayoshina · 发表于 2016-3-5 11:20:59

赞，感谢分享！

平流层的萝卜 · 发表于 2018-6-14 13:55:32

夭夭坂田淘发表于 2018-6-13 20:11
还得再麻烦一下楼主。。。我把maskout.py放在和我的脚本同一个文件夹里，但是import maskout后，仍提示“No ...

额，不好用么？那就放在你的python安装目录下的\Lib\site-packages下，再试试

xiaodixue · 发表于 2022-10-10 15:58:12

下载都是php文件，打不开怎么办，下载了好几遍

平流层的萝卜 · 发表于 2020-4-14 19:03:10

偷影子的人发表于 2020-4-14 14:46
楼主，我想画在兰伯特投影的图上，但是我改了程序之后，maskout按你说的改了，然后主程序里改了clip那里， ...

数据的投影方式这一块，你是否指 xx,yy=m(xx,yy) 这里？确保这个要写上去。

偷影子的人 · 发表于 2020-4-14 14:46:27

楼主，我想画在兰伯特投影的图上，但是我改了程序之后，maskout按你说的改了，然后主程序里改了clip那里，还有 m=Basemap(projection = 'laea', lat_0 = 33, lon_0 = 102., lat_ts = 33, width = 7000000, height = 5000000, resolution = 'l'),
但我这样写了之后，画出图来是空白的。
是什么原因呢？我看网上有人说“数据的投影方式也要改”，我对这个不是很明白，能问一下您具体应该怎么做吗？

zhouling · 发表于 2016-3-4 20:49:22

赞一个，不错

{:5_213:}

非对称 · 发表于 2016-3-4 20:51:21

支持原创，十分有用

贫道敬孔 · 发表于 2016-3-4 21:43:25

{:5_213:}虽然不用Python，但是图好漂亮

曦尘 · 发表于 2016-3-4 23:17:24

好棒，很强大👍

云知道 · 发表于 2016-3-5 00:07:48

腻害！！！

fengkp · 发表于 2016-3-5 08:39:41

你的帅气惊动了玉帝哥哥！！

星雨 · 发表于 2016-3-5 10:05:40

{:5_213:}好东西啊

牧场物语之主 · 发表于 2016-3-5 10:16:40

厉害，厉害！！

		自动登录	找回密码
密码			立即注册

[经验总结] Python完美白化

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

评分

本帖被以下淘专辑推荐:

相关帖子