我正在尝试在给定纬度和经度的情况下找到距离气象站最近的网格点。当我使用 df=df.sel(latitude=Lat.to_xarray(), longitude=Lon.to_xarray(), method='nearest') 找到最近的网格点时返回的网格点
充满了Nan的值(value)观。因此，我想找到第二个最近的网格点，希望它包含数据。我不确定如何使用上面代码的修改版本来做到这一点，所以我尝试删除作为最近的返回的原始网格点(lat = 42.36056，lon = -71.01056)，然后重新运行上面的行。我试图通过这样做来删除这一点
import os
from netCDF4 import Dataset as netcdf_dataset
import numpy as np
import xarray as xr
import pandas as pd
#open gridded data
NUM_DAYS=20
df=xr.open_mfdataset('/glacier1/mmartin/data/ERA5_LandOnly_???????.nc', chunks={'time':24*NUM_DAYS, 'latitude':271, 'longitude':601})
#drop grid point
df=df.drop_sel(latitude=['42.36056'],longitude=['-71.01056'])
但是当我这样做时，我收到以下错误:KeyError:“在轴中找不到['42.36056']”。如何删除这个网格点？或者有其他方法可以找到第二个最近的网格点吗？这是print(df)看起来像。
<xarray.Dataset>
Dimensions: (latitude: 271, longitude: 601, time: 25933)
Coordinates:
* time (time) datetime64[ns] 1951-01-01 1951-01-02 ... 2021-12-31
* longitude (longitude) float32 -125.0 -124.9 -124.8 ... -65.2 -65.1 -65.0
* latitude (latitude) float32 50.0 49.9 49.8 49.7 ... 23.3 23.2 23.1 23.0
Data variables:
t2m (time, latitude, longitude) float32 dask.array<chunksize=(1, 271, 601), meta=np.ndarray>
该数据集不是原始数据。这是我找到每日最高温度后的结果。原始数据集如下所示:
<xarray.Dataset>
Dimensions: (latitude: 271, longitude: 601, time: 613632)
Coordinates:
* longitude (longitude) float32 -125.0 -124.9 -124.8 ... -65.2 -65.1 -65.0
* latitude (latitude) float32 50.0 49.9 49.8 49.7 ... 23.3 23.2 23.1 23.0
* time (time) datetime64[ns] 1951-01-01 ... 2021-12-31T23:00:00
Data variables:
t2m (time, latitude, longitude) float32 dask.array<chunksize=(480, 271, 601), meta=np.ndarray>
Attributes:
Conventions: CF-1.6
history: 2022-10-03 03:29:52 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...
如果效果更好的话，我可以在每日最大计算之前删除该点。
最佳答案
TL;DR
您不能从数组中间删除任意点。该数组是一个[超]立方体，因此不可能从立方体中间“删除”一个点。相反，如果您尝试提取最近的非空邻居，则需要设置自定义插值器以方便提取数据。幸运的是，这并没有那么糟糕。
首先，找到要包含在插值中的有效点集。确保堆叠数据，以便您可以删除任何具有 NaN 的纬度/经度组合。然后，使用scipy.spatial.KDTree构建可重用的最近邻插值引擎，并找到要从数组中提取的最近的非空点。一旦您知道要为每个站/点提取哪个数据像素，您就可以使用 .sel 从数据中提取它们(并完全跳过 xarray 的最近邻查找)。
完整示例
设置
我将设置一个快速示例数据集:
import pandas as pd, numpy as np, xarray as xr, scipy.spatial
lons = np.arange(-109.75, -99.9, 0.5)
lats = np.arange(23.25, 28.01, 0.5)
time = pd.date_range('2020-01-01', freq='D', periods=100)
land_mask = xr.DataArray(
np.random.random(size=(10, 20)) > 0.3,
dims=['lat', 'lon'],
coords=[lats, lons],
)
da = xr.DataArray(
np.random.random(size=(10, 20, 100)),
dims=['lat', 'lon', 'time'],
coords=[lats, lons, time],
).where(land_mask)
ds = xr.Dataset({"t2m": da})
还有一个在随机位置具有“站”的 DataFrame:
stations = pd.DataFrame({
'station_id': np.arange(100000, 1000000, 10000),
'latitude': np.random.random(size=90) * 5 + 23,
'longitude': np.random.random(size=90) * 10 - 110,
}).set_index("station_id")
现在的数据看起来很像您的数据，具有固定的 NaN 空间模式和(大)时间维度:
In [3]: ds
Out[3]:
<xarray.Dataset>
Dimensions: (lat: 10, lon: 20, time: 100)
Coordinates:
* lat (lat) float64 23.25 23.75 24.25 24.75 ... 26.25 26.75 27.25 27.75
* lon (lon) float64 -109.8 -109.2 -108.8 -108.2 ... -101.2 -100.8 -100.2
* time (time) datetime64[ns] 2020-01-01 2020-01-02 ... 2020-04-09
Data variables:
t2m (lat, lon, time) float64 nan nan nan nan ... 0.9747 0.3858 0.9034
然后您想要为点列表选择与最近的非 NaN 点对应的数据:
In [5]: stations
Out[5]:
latitude longitude
station_id
100000 23.547167 -100.674304
110000 23.641703 -108.543307
120000 23.704048 -104.567338
130000 24.858875 -107.999671
140000 24.357413 -102.789371
... ... ...
950000 23.879972 -109.887476
960000 25.718888 -107.929292
970000 25.223900 -101.083424
980000 26.847443 -108.199510
990000 24.248193 -103.473922
[90 rows x 2 columns]
设置插值引擎
第一步是找到非 NaN 的点集并将它们堆叠起来，以便您可以从中提取一组有效的 x 和 y 点:
In [7]: non_null_points = ds.t2m.notnull().all(dim='time').stack(point=('lat', 'lon'))
...: non_null_points = non_null_points.where(non_null_points, drop=True)
...: non_null_points
Out[7]:
<xarray.DataArray 't2m' (point: 132)>
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
Coordinates:
* point (point) MultiIndex
- lat (point) float64 23.25 23.25 23.25 23.25 ... 27.75 27.75 27.75 27.75
- lon (point) float64 -109.2 -108.8 -107.2 ... -101.8 -100.8 -100.2
In [8]: valid_x = non_null_points.lon.values
...: valid_y = non_null_points.lat.values
现在，您可以使用scipy.spatial.KDTree构建可重用的最近邻插值引擎:
In [9]: tree = scipy.spatial.KDTree(np.stack([valid_x, valid_y]).T)
使用您的积分查询最近(有效)邻居
您现在可以使用您的站点纬度和经度查询此信息，并将相应的最近有效点分配回您的站点 DataFrame:
In [10]: dist, ind = tree.query(stations[["longitude", "latitude"]].values)
In [11]: stations["nearest_x"] = valid_x[ind]
...: stations["nearest_y"] = valid_y[ind]
In [12]: stations
Out[12]:
latitude longitude nearest_x nearest_y
station_id
100000 23.547167 -100.674304 -100.75 23.75
110000 23.641703 -108.543307 -108.75 23.75
120000 23.704048 -104.567338 -104.75 23.75
130000 24.858875 -107.999671 -108.25 24.75
140000 24.357413 -102.789371 -103.25 24.25
... ... ... ... ...
950000 23.879972 -109.887476 -109.75 23.75
960000 25.718888 -107.929292 -107.75 25.75
970000 25.223900 -101.083424 -101.25 25.25
980000 26.847443 -108.199510 -108.25 26.75
990000 24.248193 -103.473922 -103.25 24.25
[90 rows x 4 columns]
重新索引 xarray 数据集以符合您的点列表
最后，您可以使用这些最近的有效站点纬度/经度从数据中提取点:
In [13]: reindexed = ds.sel(lat=stations.nearest_y.to_xarray(), lon=stations.nearest_x.to_xarray())
In [14]: reindexed
Out[14]:
<xarray.Dataset>
Dimensions: (station_id: 90, time: 100)
Coordinates:
lat (station_id) float64 23.75 23.75 23.75 ... 25.25 26.75 24.25
lon (station_id) float64 -100.8 -108.8 -104.8 ... -108.2 -103.2
* time (time) datetime64[ns] 2020-01-01 2020-01-02 ... 2020-04-09
* station_id (station_id) int64 100000 110000 120000 ... 970000 980000 990000
Data variables:
t2m (station_id, time) float64 0.9344 0.6062 ... 0.6152 0.8736
请注意，重新索引的数据没有任何 NaN:
In [15]: reindexed.isnull().any()
Out[15]:
<xarray.Dataset>
Dimensions: ()
Data variables:
t2m bool False
关于python - 如何删除xarray中的网格点？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/74197804/