看《武林外传》的时候，最喜欢同福客栈的屋顶，每一个难过惆怅的人，都可以坐在上面，喝酒发呆，仰望星空，并且总会被人发现，然后就会有人过来陪着。屋檐下面是江湖，屋檐上只有柔软的人心。每个人的世界都会下雨，我没办法在你内心撑伞，但是在这一方小天地上，我愿意陪你....
——— 摘自网易云评论
前几天有网友私信问我，能否分享一篇绘图的白化方法（在gis行业我们习惯叫掩模裁剪或者栅格裁剪），这个需求应该是相当常见的，因为我们的感兴趣区（研究区）总是与数据源的空间范围不统一，我以前习惯于将数据转成栅格影像，然后用arcgis进行裁剪，或者写一个arcpy脚本进行批量裁剪，这一期作者再分享一种如何利用python matplotlib包进行区域掩模（白化）的方法。
image
下面我们以本地文件夹下的PERSIANN-CDR全球降水数据集（0.25°x0.25°，daily）为例，演示如何利用python在制图的时候进行感兴趣区外的白化处理。
image
首先我们展示一下PERSIANN-CDR降水数据的全球分布：
image
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : GlobalMapping.py
"""
Created by s.k zeng in Chengdu on 2020/7/12 22:23
"""
import numpy as np
import netCDF4 as nc
import cartopy.crs as ccrs
import cartopy.mpl.ticker as mticker
import matplotlib.pyplot as plt
def read_perssiancdr():
filename = r"E:\Scripts\Test\PERSIANN-CDR\PERSIANN-CDR_v01r01_20150615_c20151028.nc"
fid = nc.Dataset(filename)
print(fid.variables)
lon = fid.variables['lon'][:]
lat = fid.variables['lat'][:]
rain = fid.variables['precipitation'][:]
rain = np.squeeze(rain)
return lon, lat, rain
def main():
lon, lat, rain = read_perssiancdr()
masked_rain = np.ma.masked_where(rain <= 0, rain) # 这里是将降水量为0的区域也进行白化处理
proj = ccrs.PlateCarree(central_longitude=0)
fig = plt.figure(figsize=(9, 6), dpi=300)
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(color='k')
ax.set_global()
clevs = [0, 2.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25, 30]
cs = ax.contourf(lon, lat, np.transpose(masked_rain), clevs,
transform=proj, extend='max', cmap='coolwarm') # RdBu_r nipy_spectral
# masked out
cbar = fig.colorbar(cs, orientation='vertical', pad=0.02, aspect=20, shrink=0.6)
cbar.set_label('mm/day')
ax.set_xticks([-180, -120, -60, 0, 60, 120, 180], crs=proj)
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=proj)
lon_formatter = mticker.LongitudeFormatter(zero_direction_label=True)
lat_formatter = mticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
plt.show()
if __name__ == '__main__':
main()
接下来， 我们对该全球降水数据进行中国区域及其周边区域显示：
image
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : AreaMapping.py
"""
Created by s.k zeng in Chengdu on 2020/7/12 22:42
"""
import numpy as np
import netCDF4 as nc
import utils.maskclip as maskclip
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as mticker
import matplotlib.pyplot as plt
def read_perssiancdr():
filename = r"E:\Scripts\Test\PERSIANN-CDR\PERSIANN-CDR_v01r01_20150615_c20151028.nc"
fid = nc.Dataset(filename)
print(fid.variables)
lon = fid.variables['lon'][:]
lat = fid.variables['lat'][:]
rain = fid.variables['precipitation'][:]
rain = np.squeeze(rain)
return lon, lat, rain
def main():
lon, lat, rain = read_perssiancdr()
extent = [70, 140, 0, 60]
masked_rain = np.ma.masked_where(rain <= 0, rain)
proj = ccrs.PlateCarree(central_longitude=0)
fig = plt.figure(figsize=(9, 6), dpi=300)
ax = fig.add_subplot(1, 1, 1, projection=proj)
clevs = [0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0, 25.0, 30.0]
cs = ax.contourf(lon, lat, np.transpose(masked_rain), clevs,
transform=proj, extend='max', cmap='RdBu_r') # RdBu_r nipy_spectral
ax.add_feature(cfeature.BORDERS, linewidth=1.5, linestyle='-') # 添加国家边界线
ax.coastlines(linewidth=1.5) # 添加海岸线
ax.set_extent(extent, proj) # 设置显示范围
cbar = fig.colorbar(cs, orientation='vertical', pad=0.02, aspect=20)
cbar.set_label('mm/day')
ax.set_xticks([70, 80, 90, 100, 110, 120, 130, 140], crs=proj)
ax.set_yticks([0, 10, 20, 30, 40, 50, 60], crs=proj)
lon_formatter = mticker.LongitudeFormatter(zero_direction_label=True)
lat_formatter = mticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
plt.show()
if __name__ == '__main__':
main()
最后（今天的关键），就是给定一个shape file文件，如何将shp感兴趣区内的降水值保留，而区域外的值进行白化（mask out）处理
首先我们需要给定一个shape file文件，这里作者以中国区域shape file文件为例，进行中国区域的降水量掩模。
image
这里我们需要指定这个shape file文件中 某一个要素feature（即对于要裁剪图形）的属性和属性值（这里我们以BOU1_4M_ID = 507）为例来唯一确定该裁剪要素。
**作者编写了一个裁剪函数如下： **
说明：该裁剪函数需要输入 绘图的Figure, 正在绘图的坐标系ax, 裁剪文件shape file, 要素的属性值fieldVals。
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : maskclip.py
"""
Created by s.k zeng in Chengdu on 2020/7/12 16:55
"""
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def shp2clip(originfig, ax, shpfile, fieldVals):
"""
This method enables you to maskout the unneccessary data
outside the interest region
:param ax: the Axes instance
:param shpfile: the shape file used for clip
:param fieldVals: thi features attributes value list in shape file,
outside the region the data is to be masked out
:return:
"""
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if shape_rec.record[3] in fieldVals: # 注意这里需要指定你的字段的索引号，我的是第3个字段
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
然后编写主函数调用上面的maskclip代码，并进行绘图：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : MaskAreaMapping.py
"""
Created by s.k zeng in Chengdu on 2020/7/9 20:08
"""
import numpy as np
import netCDF4 as nc
import utils.maskclip as maskclip
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as mticker
import matplotlib.pyplot as plt
def read_perssiancdr():
"""
数据读取
:param filename:
:return:
"""
filename = r"E:\Scripts\Test\PERSIANN-CDR\PERSIANN-CDR_v01r01_20150615_c20151028.nc"
fid = nc.Dataset(filename)
print(fid.variables)
lon = fid.variables['lon'][:]
lat = fid.variables['lat'][:]
rain = fid.variables['precipitation'][:]
rain = np.squeeze(rain)
return lon, lat, rain
def main():
"""
绘图
:return:
"""
lon, lat, rain = read_perssiancdr()
extent = [70, 140, 0, 60]
masked_rain = np.ma.masked_where(rain <= 0, rain)
proj = ccrs.PlateCarree(central_longitude=0)
shpfn = r'E:\GisData\SHP\China_shp\bou2_4m\bou2_4l.shp'
reader = shpreader.Reader(shpfn)
statesFeat = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')
fig = plt.figure(figsize=(9, 6), dpi=300)
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_feature(statesFeat, linewidth=1.5, edgecolor='k')
ax.add_feature(cfeature.OCEAN)
ax.coastlines() # 添加海岸线
ax.set_extent(extent, proj)
clevs = [0, 1, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 15.0]
cs = ax.contourf(lon, lat, np.transpose(masked_rain), clevs,
transform=proj, extend='max', cmap='RdYlGn_r') # RdBu_r nipy_spectral
# masked out
maskclip.shp2clip(cs, ax, r'E:\GisData\SHP\China_shp\bou2_4m\bou2_4p.shp', [507])
cbar = fig.colorbar(cs, orientation='vertical', pad=0.02, aspect=20)
cbar.set_label('PERSIANN-CDR(mm)')
ax.set_xticks([70, 80, 90, 100, 110, 120, 130, 140], crs=proj)
ax.set_yticks([0, 10, 20, 30, 40, 50, 60], crs=proj)
lon_formatter = mticker.LongitudeFormatter(zero_direction_label=True)
lat_formatter = mticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
plt.show()
if __name__ == '__main__':
main()
最后的结果显示：
image
- 我们还可以选择shape file矢量文件的其他featrue要素的字段和属性值来进行裁剪（比如我们只需要显示四川和江苏的降水量分布）：
image
- 我们发现中国省级界线shape file文件的ADCODE99字段中四川省和江苏省的属性值为510000和320000，可以唯一确定四川和重庆的边界，因此：
image
修改后的效果图：
image
感谢支持，本人能力水平有限，如文章中存在错误和高见欢迎留言！
