¿Cómo se extrae el valor de un ráster en función de la latitud y la longitud?

Question

Tengo varios archivos netcdf que puedo leer utilizando la biblioteca netCDF4 en python. Los datos son cuadriculados en 4320 filas y 2160 columnas.

También tengo una lista de lat/lon (en grados decimales) y estoy tratando de encontrar el valor de la cuadrícula en la que se encuentran las coordenadas.

Actualmente, estoy convirtiendo las coordenadas en límites de cuadrícula mediante una simple transformación lineal y consultando el archivo netcdf utilizando estos nuevos puntos. ¿Es válido este enfoque?

# prelim
from netCDF4 import Dataset
import numpy as np

# coordinate data
lons = [40.822651, 40.69685, 40.750038, 40.798931, 40.74612]
lats = [19.83832, 20.007561, 19.97426, 19.86334, 19.843889]

# find boundaries of gridcells corresponding to lat/lon
lons_grid = np.ceil(12*(lons + 180))
lats_grid = np.ceil(12*(lats + 90))

# find value in gridcell for given lat/lon
fnc = Dataset(filename.nc, 'r')
lat = fnc.variables['latitude'][:]
lon = fnc.variables['longitude'][:]
variable = fnc.variables['variable'][:]
print variable[lons_grid, lats_grid]

Si no es así, ¿hay alguna forma de intersecar datos de puntos con un raster en python sin usar arcpy? (Mis archivos netcdf son muy grandes y no se pueden leer con arcpy).

Answer 1

Para obtener las coordenadas de las celdas adecuadas a partir de la latitud y la longitud, es necesario conocer la cobertura del archivo netCDF (por ejemplo, la coordenada de la celda superior izquierda y el tamaño de las celdas en las direcciones x e y). La solución de @Luke a continuación funciona para un raster srtaight x/y sin rotación de celdas. Si no tienes eso, podrías mirar el affine para transformar la latitud y la longitud en coordenadas de celda. Dado GeoTransformación GDAL en la forma:

1. Celda superior izquierda x coordenadas
2. x tamaño de la célula en unidades proyectadas
3. x rotación direccional (normalmente 0)
4. Celda superior izquierda y coordenadas
5. y rotación direccional (normalmente 0)
6. Negativo y tamaño de la célula en unidades proyectadas

Así, por ejemplo, el -237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0 se traduciría como una celda superior izquierda de la trama con coordenadas -237481.5, 237536.4 y un 425.0 cuadrado de la unidad sin rotación.

Utilizando el affine se puede transformar esto en un objeto de transformación así:

from affine import Affine
aff = Affine.from_gdal(-237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)

Que puede transformar las coordenadas de la celda en coordenadas proyectadas:

xs = np.array([0, 1, 2])
ys = np.array([3, 4, 5])

x_coords, y_coords = aff * (xs, ys)

O (en su caso) de coordenadas a coordenadas de celdas:

xs, ys = ~aff * (np.array(lons), np.array(lats))

Estos valores son flotantes, por lo que tendrás que transformarlos en enteros para obtener las coordenadas de las celdas que puedes utilizar.

xs = np.round(xs).astype(np.int)
ys = np.round(ys).astype(np.int)

A continuación, puede utilizarlos como índices de su matriz netCDF4 ( utilizar la última versión de netCDF4 - 1.2.1 ya que esto significa que no es necesario ordenar o eliminar los duplicados).

variable = fnc.variables['variable'][xs, ys] # be careful of your index order

Esto devolverá una matriz cuadrada debido a la indexación ligeramente diferente de netCDF, pero puede obtener los valores reales que busca como el diagnoal:

values = variable.diagonal()

Answer 2

Si conoces el tamaño de los píxeles/celdas y la longitud/longitud mínima de la cuadrícula, es un cálculo sencillo.

px = (lon -  minlon) / xcellsize
py = (lat -  minlat) / ycellsize

También puede utilizar min/max lon/lat y el número de cols/rows:

xcellsize =  (maxlon -  minlon) / ncols
ycellsize =  (maxlat -  minlat) / nrows
#then use previous calculation to get pixel coordinates