Python气象数据可视化学习笔记5——基于cartopy绘制contour并对中国地区进行白化（包含南海）

基于cartopy绘制contour并对中国地区进行白化（包含南海）

1. 写在前面

利用cartopy画填色图已经掌握，这一篇主要记录了在填色的基础上叠加白化。主要参考了气象家园的两篇帖子，并进行了整理绘图。帖子1 http://bbs.06climate.com/forum.php?mod=viewthread&tid=100563&highlight=%B0%D7%BB%AF 主要总结了论坛里的几种白化方法，并提供了maskout.py文件，本人主要采用了第四种方法，测试数据也从那获取。 帖子2 http://bbs.06climate.com/forum.php?mod=viewthread&tid=100551&highlight=cartopy 主要获取了各种shp文件，感谢大神们！

2. 效果图

先看效果图，（a）©两图给出了完整的中国地图和完整的南海九段线，其中（a）添加了海岸线；（b）(d)两图中南海被放到了右下角，这也是常用的方法。

3. 导入库

主要用的库有读取netcdf的库，绘图的cartopy和matplotlib。最重要的是maskout.py这个库，可以从帖子1中第四种方法中获取。

import os
import maskout
from netCDF4 import Dataset
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

4. 读取数据

测试数据在帖子1中获取，shp文件在帖子2中获取，不再重复。感谢两位楼主大佬。

## 读数据
f = Dataset(r'2019-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 os
import maskout
from netCDF4 import Dataset
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

##########################################################
## 读数据
f = Dataset(r'2019-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')