关于python salem掩膜处理的问题

在利用salem对nc文件的shp进行掩膜处理时，程序为：t = blh.salem.roi(shape=shp)，
然后报错OSError: exception: access violation reading 0x00000000000000A0
请问如何改正

结合报错程序的位置，参考下面文章描述排查。
https://www.codeleading.com/article/10833012124/

导入需要的库

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()

建议：
目前只支持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'}。