我正在试图找到距离气象站最近的网格点，给出一个纬度和经度。当我使用df=df.sel(latitude=Lat.to_xarray(), longitude=Lon.to_xarray(), method='nearest')找到最近的网格点时，返回的网格点充满了NAN值。正因为如此，我希望找到第二个最近的网格点，希望它包含数据。我不确定如何使用上面代码的修改版本来实现这一点，所以我尝试删除作为最近值返回的原始网格点(经度=42.36056，经度=-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’]Not Found in Axis”。如何删除此网格点？或者，有没有其他方法可以找到第二个最近的网格点？下面是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从数据中提取它们(并完全跳过x数组的最近邻查找)。

完整示例

设置

我将设置一个快速示例数据集：

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})

还有一个在随机位置有“Stations”的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]

重新索引xarrayDataSet以符合您的点列表

最后，您可以使用这些最近的有效桩号经度/经度从数据中提取点：

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

numpy 如何删除xarray中的网格点？

1条答案

TL；DR

完整示例

相关问题

热门标签

最新问答