使用 shapefile 屏蔽 NetCDF 并计算 shapefile 中所有多边形的平均值和异常值
mask NetCDF using shapefile and calculate average and anomaly for all polygons within the shapefile
有几个教程 (example 1, example 2, example 3) 关于使用 shapefile 屏蔽 NetCDF 和计算平均度量。但是,我对那些关于掩蔽 NetCDF 和提取平均值等措施的工作流程感到困惑,并且那些教程不包括提取异常(例如,2019 年的温度与基线平均温度之间的差异)。
我这里打个比方。我已经下载了每月的温度 (download temperature file) from 2000 to 2019 and the state-level US 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')
......
非常感谢您的帮助,提供了可以完成上述任务的明确工作流程。非常感谢。
这可以使用 regionmask 来实现。我不使用你的文件,而是 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 的主要作者。
有几个教程 (example 1, example 2, example 3) 关于使用 shapefile 屏蔽 NetCDF 和计算平均度量。但是,我对那些关于掩蔽 NetCDF 和提取平均值等措施的工作流程感到困惑,并且那些教程不包括提取异常(例如,2019 年的温度与基线平均温度之间的差异)。
我这里打个比方。我已经下载了每月的温度 (download temperature file) from 2000 to 2019 and the state-level US 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')
......
非常感谢您的帮助,提供了可以完成上述任务的明确工作流程。非常感谢。
这可以使用 regionmask 来实现。我不使用你的文件,而是 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 的主要作者。