python绘图 | salem一招解决所有可视化中的掩膜(Mask)问题

引言

对于空间数据，我们感兴趣的往往是其中的某一部分，对于不需要的部分需要做一些掩膜(Mask)。比如只关注海洋的数值变化，那么陆地上的数值对我其实是一种干扰，就要想办法掩盖掉。又比如我有全国的数据变量，但是只想研究其中某几个省份，那也需要对非相关省份进行掩盖。

salem一招解决掩膜

Mask的方法其实很多，但我比较懒，只想一招鲜吃遍天，掌握一个在任何场景都适用的方法，恰巧salem就满足这个需求，只需要一个shp文件就可以搞定所有问题。 salem是xarray的扩展包，集成了一些地球科学数据处理的小工具，其中.roi函数可以根据shp文件提取感兴趣的区域。

代码

全部文件及脚本后台回复salem

导入需要的库

import xarray as xr
import numpy as np
from xarray.backends import NetCDF4DataStore
import salem
from datetime import datetime
from siphon.catalog import TDSCatalog
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.io.shapereader import Reader, natural_earth
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import geopandas

获取数据

best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
                    'Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
ncss = best_gfs.datasets[0].subset()
query = ncss.query()
query.lonlat_box(north=50, south=0, east=150, west=90).time(datetime.utcnow())
query.variables('Temperature_surface')
query.accept('netcdf4')
nc = ncss.get_data(query)
data = xr.open_dataset(NetCDF4DataStore(nc))
temp = data['Temperature_surface'].isel(time=0)

定义地图函数

def create_map():
    # --创建画图空间
    proj = ccrs.PlateCarree()  # 创建坐标系
    fig = plt.figure(figsize=(6, 8), dpi=400)  # 创建页面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
    # --设置地图属性
    provinces = cfeat.ShapelyFeature(
        Reader('./cn_shp/Province_9/Province_9.shp').geometries(),
        proj, edgecolor='k',
        facecolor='none'
    )
    ax.add_feature(provinces, linewidth=0.6, zorder=2)
    ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=10)
    ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=10)
    ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=10)
     # --设置网格属性
    gl = ax.gridlines(
        crs = ccrs.PlateCarree(),
        draw_labels = False,
        linewidth = 0.9,
        color = 'k',
        alpha = 0.5,
        linestyle = '--'
    )
    gl.xlabels_top = gl.ylabels_right = gl.ylabels_left = gl.ylabels_bottom = False  # 关闭经纬度标签
    # --设置刻度
    ax.set_xticks(np.arange(90, 145 + 5, 5))
    ax.set_yticks(np.arange(0, 50 + 5, 5))
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.xaxis.set_minor_locator(plt.MultipleLocator(1))
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.yaxis.set_minor_locator(plt.MultipleLocator(1))
    ax.tick_params(axis='both', labelsize=5, direction='out')
    # -- 设置范围
    ax.set_extent([90, 140, 0, 50], crs=ccrs.PlateCarree())
    return ax

设置colorbar

cbar_kwargs = {
    'orientation': 'horizontal',
    'label': 'Potential',
    'shrink': 0.8,
}

设置level

levels = np.arange(270, 310, 1)

画图

temp.plot.contourf(
    ax=create_map(), 
    cmap='Spectral_r', 
    levels=levels, 
    cbar_kwargs=cbar_kwargs, 
    transform=ccrs.PlateCarree(), 
    extend='both'
)

读取陆地shp，并使用salem.roi来提取感兴趣的区域。

shp_path = './ne_10m_land_scale_rank/'
shp = geopandas.read_file(shp_path + 'ne_10m_land_scale_rank.shp')
t = temp.salem.roi(shape=shp)
t.plot.contourf(
    ax=create_map(), 
    cmap='Spectral_r', 
    levels=levels, 
    cbar_kwargs=cbar_kwargs, 
    transform=ccrs.PlateCarree()
)

读取海洋shp，并使用salem.roi来提取感兴趣的区域。

shp_path = './ne_10m_ocean_scale_rank/'
shp = geopandas.read_file(shp_path + 'ne_10m_ocean_scale_rank.shp')
t = temp.salem.roi(shape=shp)
t.plot.contourf(
    ax=create_map(), 
    cmap='Spectral_r', 
    levels=levels, 
    cbar_kwargs=cbar_kwargs, 
    transform=ccrs.PlateCarree()
)

读取中国各省shp，并使用salem.roi来提取感兴趣的区域。

shp_path = './China_shp_lqy/'
shp = geopandas.read_file(shp_path + 'province.shp')
shp.crs = {'init': 'epsg:4326'}
criterion = (shp['省']=='广东省') | (shp['省']=='湖南省') | (shp['省']=='福建省') | (shp['省']=='江西省')
shp = shp[criterion]
t = temp.salem.roi(shape=shp)
t.plot.contourf(
    ax=create_map(), 
    cmap='Spectral_r', 
    levels=levels, 
    cbar_kwargs=cbar_kwargs, 
    transform=ccrs.PlateCarree()

Tips

目前只支持pip安装，建议pip install salem==0.2.4.
第一次import salem的时候会自动下载salem-sample-data的zip包，但最近一直下不下来，需要把之前下好的zip包存到~/.salem_cache/路径下。(zip包获取后台回复“salem”)
用salem.roi提取区域的时候，维度名称要用标准的lat，lon，如果出现lat_0，lon_0是无法识别的。
有些shp没有投影信息，需要自行添加，如：shp.crs = {'init': 'epsg:4326'}。