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

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的帮助下从area属性读取lon lats：

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

(1179, 808)
(1179, 808)

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

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

最后，我想得到感兴趣的lat long的像素值（最近的）

提前感谢

Answer 1

您可以通过以下方式找到位置：

import numpy as np
px,py = (23.0,55.0) # some location to take out values:

dist = np.sqrt(np.cos(lat*np.pi/180.0)*(lon-px)**2+(lat-py)**2); # this is the distance matrix from point (px,py)
kkout = np.squeeze(np.where(np.abs(dist)==np.nanmin(dist))); # find location where distance is minimum
print(kkout) # you will see the row and column, where to take out data