2

I have a satellite image file. Loaded into dask array. I want to get pixel value (nearest) of a latitude, longitude of interest.

Satellite image is in GEOS projection. I have longitude and latitude information as 2D numpy arrays.

Satellite Image file

I have loaded it into a dask data array

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']

This is what the dask array looks like: enter image description here enter image description here

I read lon lats from the area attribute with the help of satpy:

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

(1179, 808)
(1179, 808)

I get a 2d numpy array each, for longitude and latitude that are coordinates but I can not use them for slicing or selecting.

What is the best practice/method to get nearest lat long, pixel information? How do I project the data onto lat long coordinates that I can then use for indexing to arrive at the pixel value.

At the end, I want to get pixel value (nearest) of lat long of interest.

Thanks in advance!!!

Rex
  • 117
  • 8
  • 1
    You have the proj parameters of the GEOS projection used by the satellite image. I would do the opposite: convert you lat lon points into that geos projection. From there a simple rounding will give you the nearest pixel(s). – Serge Ballesta Mar 02 '21 at 07:35

3 Answers3

3

The AreaDefinition object you are using (.attrs['area']) has a few methods for getting different coordinate information.

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

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

Note that row and column are flipped. The get_xy_from_lonlat method should work for arrays or scalars.

There are other methods for getting X/Y coordinates of each pixel if that is what you're interesting in.

djhoese
  • 3,567
  • 1
  • 27
  • 45
1

You can find the location with following:

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
msi_gerva
  • 2,021
  • 3
  • 22
  • 28
  • 1
    Hmm... Proj4 does a good job at handling coordinate systems and cares for the earth not being a true sphere. Why do jou suggest to compute by hand? – Serge Ballesta Mar 02 '21 at 07:31
  • Because it's simple math and arrays of lon/lat are already present. One can always take more precise/complicated methods, but I thought the question is about slicing and using the existing arrays. – msi_gerva Mar 02 '21 at 07:38
  • I can understand your point, and your are answering the exact question OP asked (I did not DV...). Yet I think that it is not the best way to handle to actual problem, the reason why I commented here and on the question. – Serge Ballesta Mar 02 '21 at 07:45
  • Thankyou guys for the comments. Thanks Serge for inputs. So the best way is to project my lat long of interest onto the GEOS projection coordinate reference system? Im looking at ways to do that now. – Rex Mar 02 '21 at 07:48
1

@serge ballesta - thanks for the direction

Answering my own question.

Project the latitude and longitude (platecaree projection) onto the GEOS projection CRS. Find x and y. Use this x and y and nearest select method of xarray to get pixel value from dask array.

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')
Rex
  • 117
  • 8