python - 使用纬度经度坐标从卫星图像获取最近的像素值

Question

我有一个卫星图像文件。加载到 dask 数组中。我想获取感兴趣的纬度、经度的像素值(最近的)。

卫星图像采用 GEOS 投影。我有 2D numpy 数组形式的经度和纬度信息。

Satellite Image file

我已将其加载到 dask 数据数组中

from satpy import Scene
import matplotlib as plt
import os

cwd = os.getcwd()

fn = os.path.join(cwd, 'EUMETSAT_data/1Jan21/MSG1-SEVI-MSG15-0100-NA-20210101185741.815000000Z-20210101185757-1479430.nat')

files = [fn]

scn = Scene(filenames=files, reader='seviri_l1b_native')
scn.load(["VIS006"])
da = scn['VIS006']

这就是 dask 数组的样子:

我在 satpy 的帮助下从区域属性中读取了 lon lats:

lon, lat = scn['VIS006'].attrs['area'].get_lonlats()
print(lon.shape)
print(lat.shape)

(1179, 808)
(1179, 808)

我每个得到一个 2d numpy 数组，用于坐标的经度和纬度，但我不能使用它们进行切片或选择。

获取最近的纬度、像素信息的最佳实践/方法是什么？如何将数据投影到经纬度坐标上，然后将其用于索引以获得像素值。

最后，我想获得感兴趣的经纬度的像素值(最接近)。

提前致谢!!!

Answer 1

您正在使用的 AreaDefinition 对象 (.attrs['area']) 有几种获取不同坐标信息的方法。

area = scn['VIS006'].attrs['area']
col_idx, row_idx = area.get_xy_from_lonlat(lons, lats)

scn['VIS006'].values[row_idx, col_idx]

请注意，行和列已翻转。 get_xy_from_lonlat 方法应该适用于数组或标量。

如果您感兴趣的话，还有其他方法可以获取每个像素的 X/Y 坐标。