基于cartopy绘制contour并对中国地区进行白化（包含南海）
1. 写在前面
利用cartopy画填色图已经掌握，这一篇主要记录了在填色的基础上叠加白化。主要参考了气象家园的两篇帖子，并进行了整理绘图。帖子1 /forum.php?mod=viewthread&tid=100563&highlight=%B0%D7%BB%AF 主要总结了论坛里的几种白化方法，并提供了maskout.py文件，本人主要采用了第四种方法，测试数据也从那获取。 帖子2 /forum.php?mod=viewthread&tid=100551&highlight=cartopy 主要获取了各种shp文件，感谢大神们！
2. 效果图
先看效果图，（a）©两图给出了完整的中国地图和完整的南海九段线，其中（a）添加了海岸线；（b）(d)两图中南海被放到了右下角，这也是常用的方法。
3. 导入库
主要用的库有读取netcdf的库，绘图的cartopy和matplotlib。最重要的是maskout.py这个库，可以从帖子1中第四种方法中获取。
import osimport maskoutfrom netCDF4 import Datasetimport numpy as npimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterfrom cartopy.io.shapereader import Readerimport matplotlib.pyplot as pltimport matplotlib.ticker as mticker
4. 读取数据
测试数据在帖子1中获取，shp文件在帖子2中获取，不再重复。感谢两位楼主大佬。
## 读数据f = Dataset(r'-06-01-1-5.nc')lat = f.variables['latitude'][:] lon = f.variables['longitude'][:]t = f.variables['t'][2,1,:,:] # 气温###########################################################白化中国地图，加南海九段线，加海岸线SHP = r'C:\Users\qiuyu\.local\share\cartopy\shapefiles\natural_earth\cultural\china_shp'
5. 定义绘图函数
这一部分的思路就是，先绘制标准地图和标准的填色图，然后利用maskout.shp2clip进行白化，把中国地图以外的区域的填色给取消。可以选择添加或者不添加海岸线。
def make_map(ax,box,lon,lat,var,proj,title,if_coast, if_nanhai):projection = ccrs.PlateCarree()# 加国界ax.add_geometries(Reader(os.path.join(SHP, 'cnmap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)# 加海岸线if if_coast:ax.add_geometries(Reader(os.path.join(SHP, 'coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)#标注坐标轴ax.set_extent([box[0],box[1],box[2],box[3]])ax.set_xticks(np.linspace(box[0], box[1],5), crs=projection) ax.set_yticks(np.linspace(box[2], box[3],5), crs=projection)#zero_direction_label=True 有度的标识，False则去掉'''lon_formatter = LongitudeFormatter(zero_direction_label=True) lat_formatter = LatitudeFormatter()ax.xaxis.set_major_formatter(lon_formatter)ax.yaxis.set_major_formatter(lat_formatter)#添加网格线ax.gridlines(linestyle='--',alpha=0.4) # plot cf = ax.contourf(lon,lat,var,cmap = mpl.cm.RdBu_r,transform=ccrs.PlateCarree()) plt.colorbar(cf,ax=ax, extend='both',orientation='vertical') maskout.shp2clip(cf,ax,shpfile=os.path.join(SHP, 'country1.shp'),region='China',proj= proj)ax.set_title(title)return ax
6. 绘图
定义经纬度区间，建立画布和轴，调用make_map函数进行画图。此时，画出的（b）（d）两图没有包含南海小地图。
box1=[70,140,0,50] #经纬度范围box2=[70,140,15,50] #经纬度范围proj=ccrs.PlateCarree()fig = plt.figure(figsize=(15,7))ax1 = fig.add_subplot(221,projection = ccrs.PlateCarree())ax2 = fig.add_subplot(222,projection = ccrs.PlateCarree())ax3 = fig.add_subplot(223,projection = ccrs.PlateCarree())ax4 = fig.add_subplot(224,projection = ccrs.PlateCarree())make_map(ax1,box1,lon,lat,t,proj,title='(a) With coastline',if_coast=True, if_nanhai=True)make_map(ax2,box2,lon,lat,t,proj,title='(b) With coastline + Nanhai',if_coast=True, if_nanhai=False)make_map(ax3,box1,lon,lat,t,proj,title='(c) Without coastline',if_coast=False, if_nanhai=False)make_map(ax4,box2,lon,lat,t,proj,title='(d) Without coastline + Nanhai',if_coast=False, if_nanhai=True)
7. 添加南海小地图
添加南海小地图的思路就是（1）先用fig.add_axes生成新的ax_nanhai插入到figure中，插入的位置由pos确定；（2）然后再重复画大地图的方法，用ax_nanhai绘制，选定新的区域（box_nanhai）。
#----------添加南海小地图------------------def add_nanhai (ax,pos,if_coast):#--------------右下角添加南海地图------------------------------------------box_nanhai=[103,125,2,25]ax_nanhai = fig.add_axes(pos,projection = ccrs.PlateCarree())# 加国界ax_nanhai.add_geometries(Reader(os.path.join(SHP, 'cnmap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)if if_coast:ax_nanhai.add_geometries(Reader(os.path.join(SHP, 'coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)ax_nanhai.set_extent([box_nanhai[0],box_nanhai[1],box_nanhai[2],box_nanhai[3]])pos1 = [0.757, 0.54, 0.1, 0.1] #南海小地图在figure中的位置和大小pos2 = [0.757, 0.124, 0.1, 0.1]add_nanhai(ax2,pos1,if_coast=True)add_nanhai(ax4,pos2,if_coast=False)plt.savefig('map.png')
8. 完整代码
import osimport maskoutfrom netCDF4 import Datasetimport numpy as npimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterfrom cartopy.io.shapereader import Readerimport matplotlib.pyplot as pltimport matplotlib.ticker as mticker############################################################ 读数据f = Dataset(r'-06-01-1-5.nc')lat = f.variables['latitude'][:] lon = f.variables['longitude'][:]t = f.variables['t'][2,1,:,:] # 气温###########################################################白化中国地图，加南海九段线，加海岸线SHP = r'C:\Users\qiuyu\.local\share\cartopy\shapefiles\natural_earth\cultural\china_shp'def make_map(ax,box,lon,lat,var,proj,title,if_coast, if_nanhai):projection = ccrs.PlateCarree()# 加国界ax.add_geometries(Reader(os.path.join(SHP, 'cnmap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)# 加海岸线if if_coast:ax.add_geometries(Reader(os.path.join(SHP, 'coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)#标注坐标轴ax.set_extent([box[0],box[1],box[2],box[3]])ax.set_xticks(np.linspace(box[0], box[1],5), crs=projection) ax.set_yticks(np.linspace(box[2], box[3],5), crs=projection)#zero_direction_label=True 有度的标识，False则去掉'''lon_formatter = LongitudeFormatter(zero_direction_label=True) lat_formatter = LatitudeFormatter()ax.xaxis.set_major_formatter(lon_formatter)ax.yaxis.set_major_formatter(lat_formatter)#添加网格线ax.gridlines(linestyle='--',alpha=0.4) # plot cf = ax.contourf(lon,lat,var,cmap = mpl.cm.RdBu_r,transform=ccrs.PlateCarree()) plt.colorbar(cf,ax=ax, extend='both',orientation='vertical') maskout.shp2clip(cf,ax,shpfile=os.path.join(SHP, 'country1.shp'),region='China',proj= proj)ax.set_title(title)return ax# make plot box1=[70,140,0,50]box2=[70,140,15,50]proj=ccrs.PlateCarree()fig = plt.figure(figsize=(15,7))ax1 = fig.add_subplot(221,projection = ccrs.PlateCarree())ax2 = fig.add_subplot(222,projection = ccrs.PlateCarree())ax3 = fig.add_subplot(223,projection = ccrs.PlateCarree())ax4 = fig.add_subplot(224,projection = ccrs.PlateCarree())make_map(ax1,box1,lon,lat,t,proj,title='(a) With coastline',if_coast=True, if_nanhai=True)make_map(ax2,box2,lon,lat,t,proj,title='(b) With coastline + Nanhai',if_coast=True, if_nanhai=False)make_map(ax3,box1,lon,lat,t,proj,title='(c) Without coastline',if_coast=False, if_nanhai=False)make_map(ax4,box2,lon,lat,t,proj,title='(d) Without coastline + Nanhai',if_coast=False, if_nanhai=True)#----------添加南海小地图------------------def add_nanhai (ax,pos,if_coast):#--------------右下角添加南海地图------------------------------------------box_nanhai=[103,125,2,25]ax_nanhai = fig.add_axes(pos,projection = ccrs.PlateCarree())# 加国界ax_nanhai.add_geometries(Reader(os.path.join(SHP, 'cnmap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)if if_coast:ax_nanhai.add_geometries(Reader(os.path.join(SHP, 'coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.7)ax_nanhai.set_extent([box_nanhai[0],box_nanhai[1],box_nanhai[2],box_nanhai[3]])pos1 = [0.757, 0.54, 0.1, 0.1]pos2 = [0.757, 0.124, 0.1, 0.1]add_nanhai(ax2,pos1,if_coast=True)add_nanhai(ax4,pos2,if_coast=False)plt.savefig('map.png')
