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

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来获取我的扫描线/ground_pixel索引:

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

这给了我(据称)距离罗马最近的地面像素的扫描线/ground_pixel坐标:

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

我将回答我自己的问题，因为我相信我想出了一个不错的解决方案，我的 blog post 上对此进行了更详细的讨论。关于这个话题。

地理距离

首先，将地球表面上两点之间的距离简单地定义为两个纬度/经度对之间的欧几里德距离可能会导致结果不准确，具体取决于两点之间的距离。因此最好将坐标转换为 ECEF首先坐标并在变换后的坐标上构建 KD 树。假设行星表面上的点(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-Tree 输入数据需要为 (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))

在这样做时，我们需要解开原始数据集形状上的索引，以获得扫描线/ground_pixel索引:

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

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

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

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

da[index]

请注意，我们还可以一次查询多个点，并且通过如上所述构建索引器，我们仍然可以通过一次调用对数据集建立索引:

da[index]

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

进一步阅读

• 在纬度/经度元组上使用欧几里得范数对于较小的距离来说可能足够准确(它使地球近似为平坦的，它在小范围内工作)。有关地理距离的更多详细信息 here .
• 使用 KD 树查找最近邻居并不是解决此问题的唯一方法，请参阅这篇非常全面的 article .
• KD-Tree 直接在 xarray 中的实现 is in the pipeline .
• 我的blog post就这个主题而言。