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

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

请注意，纬度/经度坐标标识中心像素，像素边界由 xarray 自动推断。

现在，假设我想识别距离罗马最近的地面像素。

到目前为止，我想到的最好方法是在堆叠的扁平纬度/经度数组上使用 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) 的输入。

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 树查找最近邻居并不是解决此问题的唯一方法，请参阅这篇非常全面的文章。
• KD-Tree 直接在 xarray 中的实现正在酝酿之中。
• 我关于这个主题的博客文章。