Geographic Information Systems Asked by Gpetr on January 24, 2021
I’m trying to retrieve image coordinates (X, Y) from geographic coordinates in NetCDF file (Sentinel 3 data) using Numpy. The problem is that the final pixel is too close to the desired pixel but not the closest to it (about one or two pixel far away from the closest in x, y or both x,y). I’m using Numpy’s norm()
, min()
and where()
function. The code is presented below:
from netCDF4 import Dataset
import numpy as np
# Read the NetCDF file with the geographic coordinates information
nc_coord = Dataset('coords.nc', 'r')
lat = nc_coord.variables['latitude'][:]
# Python buffer object pointing to the start of the array’s data
lat = lat.data
lon = nc_coord.variables['longitude'][:]
lon = lon.data
# define your point:
target_lon, target_lat = 23.779675, 35.33786
# lat and lon are x,y matrices, we will add another dimension z in order to stack them
lat = lat[:, :, np.newaxis]
lon = lon[:, :, np.newaxis]
grid = np.concatenate([lat, lon], axis=2)
vector = np.array([target_lat, target_lon]).reshape(1, 1, -1)
subtraction = vector - grid
# vector subtraction
dist = np.linalg.norm(subtraction, axis=2)
# query where in the 2D matrix is the closest point to your points of interest (?)
result = np.where(dist == dist.min())
# this is your lat/lon geo point translated to a x, y position in the image matrix
target_x_y = result[0][0], result[1][0]
Can you figure out what is happening?
From your question, you are using a swath dataset (i.e. not using a regular grid). If you were using a regular grid, you'd just calculate the pixel number from the e.g. corner coordinates. I think
geolocation = Dataset(fname)
points = np.array([geolocation.variables['lon'][:],
geolocation.variables['lat'][:]
])
# Exact point
point = np.array([3.60829699, 44.41841389])
dist = np.linalg.norm(point[:, None, None] - points, axis=0)
i, j = np.where(dist == dist.min())
# Should evaluate to True
print(np.allclose(point, points[:, i, j].squeeze()))
print(dist.min())
# Not Exact point, actually closest to neighbouring pixel
point = np.array([3.605, 44.421])
dist = np.linalg.norm(point[:, None, None] - points, axis=0)
i, j = np.where(dist == dist.min())
# Should evaluate to False, as test point not in arrays
print(np.allclose(point, points[:, i, j].squeeze()))
print(dist.min())
Answered by Jose on January 24, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP