从 NETCDF 文件中提取数据的有效方法
Efficient way to extract data from NETCDF files
我有一些坐标(大约 20000 个),我需要从多个 NetCDF 文件中提取数据,每个文件大约有 30000 个时间步长(未来气候情景)。使用解决方案 效率不高,原因是在每个 i,j 将“dsloc”转换为“dataframe”所花费的时间(查看下面的代码)。
** 示例 NetCDF 文件可以从 here 下载 **
import pandas as pd
import xarray as xr
import time
#Generate some coordinates
coords_data = [{'lat': 68.04, 'lon': 15.20, 'stid':1},
{'lat':67.96, 'lon': 14.95, 'stid': 2}]
crd= pd.DataFrame(coords_data)
lat = crd["lat"]
lon = crd["lon"]
stid=crd["stid"]
NC = xr.open_dataset(nc_file)
point_list = zip(lat,lon,stid)
start_time = time.time()
for i,j,id in point_list:
print(i,j)
dsloc = NC.sel(lat=i,lon=j,method='nearest')
print("--- %s seconds ---" % (time.time() - start_time))
DT=dsloc.to_dataframe()
DT.insert(loc=0,column="station",value=id)
DT.reset_index(inplace=True)
temp=temp.append(DT,sort=True)
print("--- %s seconds ---" % (time.time() - start_time))
结果是:
68.04 15.2
--- 0.005853414535522461 seconds ---
--- 9.02660846710205 seconds ---
67.96 14.95
--- 9.028568267822266 seconds ---
--- 16.429600715637207 seconds ---
这意味着每个 i,j 需要大约 9 秒的时间来处理。给定大量时间步长的坐标和 netcdf 文件,我想知道是否有一种 pythonic 方式可以优化代码。
我也可以使用 CDO 和 NCO 运算符,但我也发现了类似的问题。
我有一个可能的解决方案。这个想法是首先将 xarray 数据数组转换为 pandas,然后根据 lat/lon 条件获取 pandas 数据帧的子集。
# convert xarray data to a pandas dataframe
def xr_to_df(data):
data = data.to_dataframe()
data.reset_index(inplace=True)
return data
# convert your xarray data to a pandas dataframe
full_df = xr_to_df(full_xarray)
# create a 2 columns pandas dataframe containing your target coordinates
points = pd.DataFrame({'lat':target_lat, 'lon':target_lon})
# get the values at your target points only via merging on the left
subset = pd.merge(points,full_df)
对于您的数据量,我不确定这有多快。但至少,这避免了循环。我想它应该更快?
我注意到你的点是随机分布的(不在网格中心)。要解决这个问题,您可以先编写自己的代码将它们重新网格化到 netcdf 分辨率上,使用 np.argmin(abs(lat - lat_netcdf))
之类的东西找到最近的纬度和经度。
这是 xarray advanced indexing 使用 DataArray 索引的完美用例。
# Make the index on your coordinates DataFrame the station ID,
# then convert to a dataset.
# This results in a Dataset with two DataArrays, lat and lon, each
# of which are indexed by a single dimension, stid
crd_ix = crd.set_index('stid').to_xarray()
# now, select using the arrays, and the data will be re-oriented to have
# the data only for the desired pixels, indexed by 'stid'. The
# non-indexing coordinates lat and lon will be indexed by (stid) as well.
NC.sel(lon=crd_ix.lon, lat=crd_ix.lat, method='nearest')
数据中的其他维度将被忽略,因此如果您的原始数据具有维度 (lat, lon, z, time)
,您的新数据将具有维度 (stid, z, time)
。
我有一些坐标(大约 20000 个),我需要从多个 NetCDF 文件中提取数据,每个文件大约有 30000 个时间步长(未来气候情景)。使用解决方案
import pandas as pd
import xarray as xr
import time
#Generate some coordinates
coords_data = [{'lat': 68.04, 'lon': 15.20, 'stid':1},
{'lat':67.96, 'lon': 14.95, 'stid': 2}]
crd= pd.DataFrame(coords_data)
lat = crd["lat"]
lon = crd["lon"]
stid=crd["stid"]
NC = xr.open_dataset(nc_file)
point_list = zip(lat,lon,stid)
start_time = time.time()
for i,j,id in point_list:
print(i,j)
dsloc = NC.sel(lat=i,lon=j,method='nearest')
print("--- %s seconds ---" % (time.time() - start_time))
DT=dsloc.to_dataframe()
DT.insert(loc=0,column="station",value=id)
DT.reset_index(inplace=True)
temp=temp.append(DT,sort=True)
print("--- %s seconds ---" % (time.time() - start_time))
结果是:
68.04 15.2
--- 0.005853414535522461 seconds ---
--- 9.02660846710205 seconds ---
67.96 14.95
--- 9.028568267822266 seconds ---
--- 16.429600715637207 seconds ---
这意味着每个 i,j 需要大约 9 秒的时间来处理。给定大量时间步长的坐标和 netcdf 文件,我想知道是否有一种 pythonic 方式可以优化代码。 我也可以使用 CDO 和 NCO 运算符,但我也发现了类似的问题。
我有一个可能的解决方案。这个想法是首先将 xarray 数据数组转换为 pandas,然后根据 lat/lon 条件获取 pandas 数据帧的子集。
# convert xarray data to a pandas dataframe
def xr_to_df(data):
data = data.to_dataframe()
data.reset_index(inplace=True)
return data
# convert your xarray data to a pandas dataframe
full_df = xr_to_df(full_xarray)
# create a 2 columns pandas dataframe containing your target coordinates
points = pd.DataFrame({'lat':target_lat, 'lon':target_lon})
# get the values at your target points only via merging on the left
subset = pd.merge(points,full_df)
对于您的数据量,我不确定这有多快。但至少,这避免了循环。我想它应该更快?
我注意到你的点是随机分布的(不在网格中心)。要解决这个问题,您可以先编写自己的代码将它们重新网格化到 netcdf 分辨率上,使用 np.argmin(abs(lat - lat_netcdf))
之类的东西找到最近的纬度和经度。
这是 xarray advanced indexing 使用 DataArray 索引的完美用例。
# Make the index on your coordinates DataFrame the station ID,
# then convert to a dataset.
# This results in a Dataset with two DataArrays, lat and lon, each
# of which are indexed by a single dimension, stid
crd_ix = crd.set_index('stid').to_xarray()
# now, select using the arrays, and the data will be re-oriented to have
# the data only for the desired pixels, indexed by 'stid'. The
# non-indexing coordinates lat and lon will be indexed by (stid) as well.
NC.sel(lon=crd_ix.lon, lat=crd_ix.lat, method='nearest')
数据中的其他维度将被忽略,因此如果您的原始数据具有维度 (lat, lon, z, time)
,您的新数据将具有维度 (stid, z, time)
。