寻找给定坐标在不规则网格上最近的地面像素

Question

我使用的是不规则二维网格组织的卫星数据，其维度为扫描线（沿轨道方向）和地面像素（横向跨越方向）。每个地面像素的纬度和经度信息存储在辅助坐标变量中。

给定一个（纬度，经度）点，我想要确定我的数据集中最接近的地面像素。

让我们构建一个10x10的玩具数据集：

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
%matplotlib inline

lon, lat = np.meshgrid(np.linspace(-20, 20, 10), 
                       np.linspace(30, 60, 10))
lon += lat/10
lat += lon/10
da = xr.DataArray(data = np.random.normal(0,1,100).reshape(10,10), 
                  dims=['scanline', 'ground_pixel'],
                  coords = {'lat': (('scanline', 'ground_pixel'), lat),
                            'lon': (('scanline', 'ground_pixel'), lon)})

ax = plt.subplot(projection=ccrs.PlateCarree())
da.plot.pcolormesh('lon', 'lat', ax=ax, cmap=plt.cm.get_cmap('Blues'), 
                   infer_intervals=True);
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.tight_layout()

请注意，lat/lon坐标标识中心像素，并且像素边界由xarray自动推断。

现在，假设我想要确定离罗马最近的地面像素。

到目前为止，我想到的最好方法是在堆叠的平铺的lat/lon数组上使用scipy的kdtree：

from scipy import spatial
pixel_center_points = np.stack((da.lat.values.flatten(), da.lon.values.flatten()), axis=-1)
tree = spatial.KDTree(pixel_center_points)

rome = (41.9028, 12.4964)
distance, index = tree.query(rome)
print(index)
# 36

然后我需要使用`unravel_index`函数来获取我的扫描线/地面像素索引：

pixel_coords = np.unravel_index(index, da.shape)
print(pixel_coords)
# (3, 6)

这段话的意思是：“这里给出了扫描线/地面像素坐标，即与罗马最近（理论上）的地面像素。”

ax = plt.subplot(projection=ccrs.PlateCarree())
da.plot.pcolormesh('lon', 'lat', ax=ax, cmap=plt.cm.get_cmap('Blues'), 
                   infer_intervals=True);
ax.scatter(da.lon[pixel_coords], da.lat[pixel_coords], 
           marker='x', color='r', transform=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.tight_layout()

我相信一定有更优雅的方法来解决这个问题。特别是，我想摆脱展平/解开步骤（我尝试在二维数组上构建kdtree的所有尝试都失败了），尽可能利用我的xarray数据集的变量（例如添加新的center_pixel维度，并将其用作输入到KDTree）。

- stm4tt

Answer 1

我将回答自己的问题，因为我相信我想出了一个不错的解决方案，该方案在我关于这个主题的博客文章中有更详细的讨论。

地理距离

首先，将地球表面两点之间的距离定义为两个纬度/经度对之间的欧几里得距离可能会导致不准确的结果，具体取决于两个点之间的距离。因此最好先将坐标转换为ECEF坐标，并在转换后的坐标上构建KD-Tree。假设点位于行星表面（h=0），则进行坐标转换如下：

def transform_coordinates(coords):
    """ Transform coordinates from geodetic to cartesian

    Keyword arguments:
    coords - a set of lan/lon coordinates (e.g. a tuple or 
             an array of tuples)
    """
    # WGS 84 reference coordinate system parameters
    A = 6378.137 # major axis [km]   
    E2 = 6.69437999014e-3 # eccentricity squared    

    coords = np.asarray(coords).astype(np.float)

    # is coords a tuple? Convert it to an one-element array of tuples
    if coords.ndim == 1:
        coords = np.array([coords])

    # convert to radiants
    lat_rad = np.radians(coords[:,0])
    lon_rad = np.radians(coords[:,1]) 

    # convert to cartesian coordinates
    r_n = A / (np.sqrt(1 - E2 * (np.sin(lat_rad) ** 2)))
    x = r_n * np.cos(lat_rad) * np.cos(lon_rad)
    y = r_n * np.cos(lat_rad) * np.sin(lon_rad)
    z = r_n * (1 - E2) * np.sin(lat_rad)

    return np.column_stack((x, y, z))

构建KD树

我们可以通过转换数据集的坐标来构建KD树，需要注意将2D网格平铺为一维的纬度/经度元组序列。这是因为KD树的输入数据需要是（N，K），其中N是点的数量，K是维度（在我们的情况下，K=2，因为我们假设没有高度分量）。

# reshape and stack coordinates
coords = np.column_stack((da.lat.values.ravel(),
                          da.lon.values.ravel()))

# construct KD-tree
ground_pixel_tree = spatial.cKDTree(transform_coordinates(coords))

查询树和索引xarray数据集

现在查询树的方法就是将点的经纬度坐标转换为ECEF坐标并传递给树的query方法：

rome = (41.9028, 12.4964)
index = ground_pixel_tree.query(transform_coordinates(rome))

在这样做的过程中，我们需要解开原始数据集形状上的索引，以获取扫描线/地面像素索引：

index = np.unravel_index(index, self.shape)

我们现在可以使用这两个组件来索引原始的xarray数据集，但我们也可以构建两个索引器来与xarray 点对点索引 功能一起使用：

index = xr.DataArray(index[0], dims='pixel'), \
        xr.DataArray(index[1], dims='pixel')

获取最接近的像素现在既简单又优雅：

da[index]

请注意，我们也可以一次查询多个点，并通过以上方式构建索引器，仍然可以使用单个调用对数据集进行索引。

da[index]

然后，它将返回数据集的子集，其中包含最接近我们查询点的地面像素。

进一步阅读

• 对于较小的距离，使用经纬度元组上的欧几里得范数可能足够准确（将地球近似为平面，这在小尺度上有效）。有关地理距离的更多详细信息请参见此处。
• 使用KD-Tree查找最近邻居并不是解决此问题的唯一方法，请参见这篇非常全面的文章。
• 直接将KD-Tree实现到xarray中正在进行中。
• 我关于此主题的博客文章。

- stm4tt

---

正如您所指出的那样，这个功能已经在xarray的范围内，只等着一个有兴趣的用户/开发者来实现它。看起来您已经在努力实现它了。如果您想在包中看到此功能，我鼓励您向xarray项目提交PR。 - jhamman