运用Xarray读取nc文件并进行数组运算的小bug:
问题描述
我这里主要是针对一个长时间序列的三维数组针对后面两个维度进行加减运算,在这个过程中发现利用xarray进行运算时会报错。
import sys
import warnings
warnings.filterwarnings("ignore")
import os
import numpy as np
import geopandas as gpd
import xarray as xr
import salem
shp_path = "/Data/LIPING/cordex/dataset/masktp/TPBoundary_HF/TPBoundary_HF_LP.shp"
tp_shp = gpd.read_file(shp_path)
hgt = xr.open_dataset('CN05-hgt.nc').topo
hgt_ma = hgt.salem.roi(shape=tp_shp)
ds = xr.open_dataset(f'masktp_cn05_tas_daily_1970-2005.nc')
tas = ds.tas
tas_hc = tas - 0.0065 * hgt_ma
解决方案1:
经过尝试我发现是因为hgt存储的经纬度的属性为:latitude和longitude而,tas中存储的经纬度的属性为:lat和lon,因此我对topo这个变量重新进行属性赋值之后发现可以成功进行运算了。
import sys
import warnings
warnings.filterwarnings("ignore")
import os
import numpy as np
import geopandas as gpd
import xarray as xr
import salem
shp_path = "/Data/LIPING/cordex/dataset/masktp/TPBoundary_HF/TPBoundary_HF_LP.shp"
tp_shp = gpd.read_file(shp_path)
lon = np.arange(69.75, 140.5, 0.25)
lat = np.arange(14.75, 55.5, 0.25)
hgt = xr.open_dataset('CN05-hgt.nc').topo
hgt = xr.DataArray(hgt,coords=[lat,lon], dims=['lat','lon'])
#如果要让后续的代码能够顺利的进行运算,必须使得两个DataArray的后两个维度的属性名称和大小一致
hgt_ma = hgt.salem.roi(shape=tp_shp)
ds = xr.open_dataset(f'masktp_cn05_tas_daily_1970-2005.nc')
tas = ds.tas
tas_hc = tas - 0.0065 * hgt_ma
解决方案2:
将所有的数组转换城array运用numpy进行数组运算。
import sys
import warnings
warnings.filterwarnings("ignore")
import os
import numpy as np
import geopandas as gpd
import xarray as xr
import salem
shp_path = "/Data/LIPING/cordex/dataset/masktp/TPBoundary_HF/TPBoundary_HF_LP.shp"
tp_shp = gpd.read_file(shp_path)
hgt = xr.open_dataset('CN05-hgt.nc').topo
hgt_ma = hgt.salem.roi(shape=tp_shp)
ds = xr.open_dataset(f'masktp_cn05_tas_daily_1970-2005.nc')
tas = ds.tas
tas = np.array(tas)
hgt_ma = np.array(hgt_ma)
tas_hc = np.zeros_like(tas)
for i in range(tas.shape[0]):
tas_hc[i,:,:] = tas[i,:,:] - 0.0065 * hgt_ma