使用 shapefile 屏蔽 NetCDF 并计算 shapefile 中所有多边形的平均值和异常值

Posted 2023-02-18

【中文标题】使用 shapefile 屏蔽 NetCDF 并计算 shapefile 中所有多边形的平均值和异常值

【英文标题】：mask NetCDF using shapefile and calculate average and anomaly for all polygons within the shapefile

【发布时间】：2022-01-11 09:32:24

【问题描述】：

有几个教程（example 1、example 2、example 3）关于使用 shapefile 屏蔽 NetCDF 和计算平均度量。但是，我对那些关于屏蔽 NetCDF 和提取平均值等度量的工作流程感到困惑，并且这些教程不包括提取异常（例如，2019 年的温度与基线平均温度之间的差异）。

我在这里举个例子。我下载了 2000 年到 2019 年的月度温度 (download temperature file) 和美国州级 shapefile (download shapefile)。我想根据 2000 年至 2019 年的月平均温度以及 2019 年相对于 2000 年至 2010 年的基准温度的温度异常来获得州级平均温度。具体而言，最终数据框如下所示：

state	avg_temp	anom_temp2019
AL	xx	xx
AR	xx	xx
...	...	...
WY	xx	xx

# Load libraries
%matplotlib inline

import regionmask
import numpy as np
import xarray as xr
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

# Read shapefile
us = gpd.read_file('./shp/state_cus.shp')

# Read gridded data
ds = xr.open_mfdataset('./temp/monthly_mean_t2m_*.nc')
......

非常感谢您提供的可以完成上述任务的明确工作流程的帮助。非常感谢。

【参考方案1】：

这可以使用区域掩码来实现。我不使用您的文件，而是使用美国各州的 xarray 教程数据和 naturalearth 数据。

import numpy as np
import regionmask
import xarray as xr

# load polygons of US states
us_states_50 = regionmask.defined_regions.natural_earth.us_states_50

# load an example dataset
air = xr.tutorial.load_dataset("air_temperature")

# turn into monthly time resolution
air = air.resample(time="M").mean()

# create a mask
mask3D = us_states_50.mask_3D(air)

# latitude weights
wgt = np.cos(np.deg2rad(air.lat))

# calculate regional averages
reg_ave = air.weighted(mask3D * wgt).mean(("lat", "lon"))

# calculate the average temperature (over 2013-2014)
avg_temp = reg_ave.sel(time=slice("2013", "2014")).mean("time")

# calculate the anomaly (w.r.t. 2013-2014)
reg_ave_anom = reg_ave - avg_temp

# select a single timestep (January 2013)
reg_ave_anom_ts = reg_ave_anom.sel(time="2013-01")

# remove the time dimension
reg_ave_anom_ts = reg_ave_anom_ts.squeeze(drop=True)

# convert to a pandas dataframe so it's in tabular form
df = reg_ave_anom_ts.air.to_dataframe()

# set the state codes as index
df = df.set_index("abbrevs")

# remove other columns
df = df.drop(columns="names")

您可以在regionmask docs (Working with geopandas) 上找到如何使用自己的 shapefile 的信息。

免责声明：我是regionmask的主要作者。

【讨论】：

感谢您的帮助。我对权重不是很清楚，为什么要使用权重，为什么要使用基于纬度的权重，我们可以使用基于经度的权重，如果可以，基于纬度和基于经度的权重有什么区别？

这取决于您使用的网格。当您向极点移动时，网格单元会变小。在常规的 lat-lon 网格上，它们使用 lat 的余弦进行缩放，对于非常规网格则不同。另一方面，经度在任何地方都是相同的长度。这有帮助吗？

是的，这是有道理的。我有另一个问题。 us_states_50 有 51 个区域（或者，也许我可以称它为状态），但是，在屏蔽之后，mask3D 只有 43 个区域（缺少区域 6、7、8、11、20、30、39 和 49），为什么会有缺少几个区域？

不包含任何点的区域将被删除（例如夏威夷）。您可以通过设置mask3D = us_states_50.mask_3D(air, drop=False) 来保留这些区域。在df 中，它们将显示为NaN 值。

谢谢。与此相关的另一个问题。对我来说，一个 NC 文件由网格单元（由 lat-lon 对表示）组成，每个单元存储几个值（例如，温度、降水等）。这意味着每个区域都应该有一些细胞落入其中，对吗？为什么有些区域不包含细胞（或点）？也许区域太小了？