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

Question

我有卫星图像文件。加载到 dask 阵列中。我想获得感兴趣的纬度、经度的像素值（最近）。

卫星图像在 GEOS 投影中。我有经度和纬度信息作为 2D numpy 数组。

我已经将它加载到 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 数组，经度和纬度是坐标，但我不能使用它们进行切片或选择。

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

最后，我想获得感兴趣的经纬度的像素值（最近）。

提前致谢！！！

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

Answer 2

@serge ballesta - 感谢指导

回答我自己的问题。

将纬度和经度（platecaree投影）投影到GEOS投影CRS上。找到 x 和 y。使用 x 和 y 以及 xarray 的最近 select 方法从 dask 数组中获取像素值。

import cartopy.crs as ccrs

data_crs = ccrs.Geostationary(central_longitude=41.5, satellite_height=35785831, false_easting=0, false_northing=0, globe=None, sweep_axis='y')

lon = 77.541677 # longitude of interest
lat = 8.079148 # latitude of interst

# lon lat system in 
x, y = data_crs.transform_point(lon, lat, src_crs=ccrs.PlateCarree())

dn = ds.sel(x=x,y=y, method='nearest')

Answer 3

您正在使用的 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 坐标。

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

Get nearest pixel value from satellite image using latitude longitude coordinates

python

gis

netcdf

python-xarray

satellite-image