从NetCDF中提取局部数据

发表于 2019-07-03 分类于闭门造车阅读次数： Disqus：

前言

前面我们已经简单介绍了NetCDF文件以及如何用Python处理NetCDF数据。在前面的示例中，我们可以看到NetCDF经常是全球的数据，那么如果我们只希望针对其中一部分，就像剪纸一样提取特定区域或者形状的数据，该如何操作呢？这里我们就简单介绍从NetCDF中裁取特定区域的子集所使用的两个Python包：salem和regionmask。

说起剪纸，我第一反应就是动画片《狐狸打猎人》，小时候看故事书的时候并没有觉得什么，看动画片的时候就觉得有点慌了，现在回想起来都有一股寒意……

fox_hunter

使用salem提取NetCDF子集

salem实际上功能较为简单，相当于对xarray库增加了选取子集、画图等方面的功能，使用起来也很容易上手。

1	%matplotlib inline

import matplotlib.pyplot as plt
import matplotlib as mpl
import xarray
import pandas as pd
import numpy as np
import salem
import cartopy.crs as ccrs
import regionmask
import cartopy.feature as cfeature
from salem.utils import get_demo_file

mpl.rcParams['figure.figsize']=(16,9)
mpl.style.use('ggplot')

准备数据

我们准备的NetCDF文件是由NOAA提供的再分析资料中的2018年近地面层（sig995）的空气温度，我们可以用salem，实际上是用xarray打开，数据的具体信息如下：

1 2	data = salem.open_xr_dataset(r'air.sig995.2018.nc') data

<xarray.Dataset>
Dimensions:  (lat: 73, lon: 144, time: 1460)
Coordinates:
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) datetime64[ns] 2018-01-01 ... 2018-12-31T18:00:00
Data variables:
    air      (time, lat, lon) float32 ...
Attributes:
    Conventions:    COARDS
    title:          4x daily NMC reanalysis (2014)
    history:        created 2017/12 by Hoop (netCDF2.3)
    description:    Data is from NMC initialized reanalysis\n(4x/day).  These...
    platform:       Model
    References:     http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reana...
    dataset_title:  NCEP-NCAR Reanalysis 1
    pyproj_srs:     +units=m +proj=latlong +datum=WGS84

虽然再分析资料是全球的网格数据（ $2.5\times2.5$ ），但是我们只想看中国的怎么办？所以我们需要用到中国的shp文件来提取，shp文件可以从Natural Earth下载。

1
2
3

# GeoPandas' GeoDataFrame
df_shp = salem.read_shapefile(r'D:\OneDrive\PyProjects\ZGis\NEV\10m_cultural\ne_10m_admin_0_countries.shp')
china_shp = df_shp.loc[df_shp['name'].isin(['China'])]

提取数据

salem中与提取数据相关的两个方法分别是.salem.subset()和.salem.roi()，其中前者提取的是shp文件所在的一个矩形，而后者则是限定在感兴趣区域（region of interest, roi）边界周围，其余的数据都会被剔除（masked）。

1 2	subset = data.salem.subset(shape=china_shp, margin=1) roi = subset.salem.roi(shape=china_shp, all_touched=True)

作图

如前所述，salem为Dataset和DataArray添加了作图功能：

f, axes = plt.subplots(2,1,figsize=(20,16))
roi.air.isel(time=1).salem.quick_map(ax=axes[0])
subset.air.isel(time=1).salem.quick_map(ax=axes[1])
plt.show()

output_16_0

使用regionmask提取NetCDF子集

salem的开发者在技术文档中提到，regionmask提供与salem类似的功能，而且会更加通用一些。的确如此，regionmask可以直接读取shapely生成的Polygon，所以功能也就更加底层（但是似乎regionmask不能够直接读取shapefile)。

准备数据

这里我们使用CMIP5模式在RCP8.5情景下的降水数据作为示例，具体的含义我也是不求甚解，有兴趣的可以查一下相关资料，数据的相关信息如下：

1 2	data = xarray.open_dataset(r'data\PATTERN_pr_MONS_ACCESS1-3_rcp85.nc') data

<xarray.Dataset>
Dimensions:      (lat: 145, lon: 192, month: 12)
Coordinates:
  * lat          (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 87.5 88.75 90.0
  * lon          (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1
Dimensions without coordinates: month
Data variables:
    pattern      (month, lat, lon) float32 ...
    CI95         (month, lat, lon) float32 ...
    r2           (month, lat, lon) float32 ...
    error        (month, lat, lon) float32 ...
    climatology  (month, lat, lon) float32 ...
Attributes:
    title:              LSR pattern scaling
    creation_date:      Mon Feb 13 16:35:45 EST 2017
    methods:            PRlocal regressed onto TGAV via least squared regression
    TGAV:               Globally averaged & weighted by latitude monthly TAS
    time_increment:     Months
    source_model:       ACCESS1-3
    original_variable:  pr
    realm:              Amon
    NCL_function1:      wgt_areaave (weighted area average)
    NCL_function2:      regCoef (LSR regression)
    Lat_weights:        sqrt(cos(rad*lat))
    editor:             CLynch
    editor_locale:      JGCRI
    forcing:            rcp85
    months:             start month = Jan; end month = Dec
    _FillValue:         1e+20
    units:              Celsius
    cell_methods:       time: mean
    cell_measures:      area: areacella
    history:            2012-03-14T03:46:55Z altered by CMOR: Treated scalar ...
    missing_value:      1e+20
    associated_files:   baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocatio...

我们可以看一下对6月份降水的模式模拟结果：

proj = ccrs.Robinson(central_longitude=100)
ax = plt.axes(projection=proj)
im = data.pattern.isel(month=5).plot.pcolormesh(ax=ax,
                                                transform=ccrs.PlateCarree(),
                                                cmap='coolwarm',
                                                add_colorbar=False)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
cbar = plt.colorbar(im,fraction=0.024, pad=0.04)
cbar.set_label('pattern[mm/month per TGAV in $^o C$]')
plt.show()

output_23_0

在选择区域的时候，regionmask大概可以分为四种方法：

利用NaturalEarth提供的数据，提取国家或者陆地区域；
regionmask本身定义了在科学文献中常用的地理分区，包括Giorgi分区和SREX分区
使用numpy或者xarray中的区域；
开发者使用shapely等自己创建感兴趣区域。

示例介绍

Scientific Regions

看到regionmask的说明文档以后，我对其中的Scientific Regions比较感兴趣，就简单了解了一下，并且使用SREX定义的分区：SREX分区出自联合国政府间气候变化专门委员会（IPCC）关于极端事件风险管理和环境变化的报告，将全球主要的人类居住地分为了26个地理区域，具体分区见下：

1
2
3

from regionmask.defined_regions import giorgi, natural_earth, srex
regionmask.defined_regions.srex.plot()
plt.tight_layout()

output_28_0

从上面我们可以看到，掩码的顺序是按从阿拉斯加开始，按照从北到南、自西向东的顺序编码，那么我们可以从中提取亚非拉和大洋洲的区域作为ROI：

mask = srex.mask(data, wrap_lon=True)
ax = plt.axes(projection=proj)
ptn_regs = data.where(mask > 10)
im = ptn_regs.sel(month=5).pattern.plot.pcolormesh(ax=ax,
                                                   transform=ccrs.PlateCarree(),
                                                   cmap='coolwarm',
                                                   add_colorbar=False)
srex.plot(ax=ax, add_ocean=False, coastlines=False, label='abbrev')
ax.coastlines()
cbar = plt.colorbar(im,fraction=0.024, pad=0.04)
cbar.set_label('pattern[mm/month per TGAV in $^o C$]')
plt.tight_layout()

output_30_0

基于Natural Earth

mask = natural_earth.countries_50[['Australia', 'Argentina',
                                   'RUS', 'BR', 'CA', 'CN',
                                   'Taiwan', 'US']].mask(data, wrap_lon=True)

ptn_con = data.where(mask>=0)

pro = ccrs.Robinson(central_longitude=180)
ax = plt.axes(projection=proj)
ptn_con.sel(month=5).pattern.plot.pcolormesh(ax=ax,
                                             transform=ccrs.PlateCarree(),
                                             cmap='coolwarm',
                                             add_colorbar=False)

natural_earth.countries_110.plot(ax=ax, add_ocean=False, resolution='50m',
                                 proj=proj, add_label=False)
cbar = plt.colorbar(im,fraction=0.024, pad=0.04)
cbar.set_label('pattern[mm/month per TGAV in $^o C$]')
plt.show()

output_32_0

自定义提取区域

个人觉得能够自定义提取区域是最重要的功能，在regionmask中与自定义提取区域的方法是regionmask.Regions_cls()和regionmask.Region_cls()，两者的区别从名字上可以瞧出一些端倪：前者接受列表作为参数，一次可以定义多个区域，而后者是提取单个区域。

接下来，我们按照CORDEX (COordinated Regional Downscaling Experiment) 计划中定义的"domain"作为示例。

CORDEX一共定义了14个domains，一个domain是覆盖一个区域的多边形，由一系列参数定义（例如旋转坐标系的旋转极点等、非旋转坐标系下的8个顶点和1个中心点等），感兴趣的可以查阅相关资料。

from regionmask import Regions_cls

outline = dict()
outline[1] = ((51.59, 50.50), (116.70, 51.90), (181.50, 50.31), (165.94, 25.56),
               (156.08, -0.24), (116.51, 6.90), (76.91, -0.10), (67.11, 25.72))

midpoint = dict()
midpoint[1] = (116.57, 34.40)

short_name = {1: 'eas'}
name = {1: 'East Asia'}

numbers = range(1, 2)
source = ('https://www.cordex.org/domains/')
cordex = Regions_cls('CORDEX', numbers, name,
                     short_name, outline,source=source)

ax = plt.axes(projection=ccrs.Orthographic(central_longitude=120, central_latitude=25))
ax.coastlines()
ax.set_global()
cordex[['eas']].plot(ax=ax, label='name')
plt.show()

output_35_0

小结

在NetCDF入门的基础上，本文介绍了如何使用salem和regionmask从NetCDF文件中按照区域提取数据。经过这些日子的尝试，Python下关于地理分析的库非常多，许多在功能上会有所重复，但是又各有特色，但是最为麻烦的是，有时候不同包的依赖可能会存在冲突，所以在使用的时候要谨慎选择、对症下药。如果有一天，能有一个更加通用和强大的包就好了，然而我感觉这个领域有些小众，短时间内是看不到希望了。

参考文献

fox_hunter_2