Extraer datos de NetCDF utilizando un polígono shapefile en Python

Question

Necesito hacer un subconjunto de datos de un NetCDF utilizando un shapefile específico. Los datos son la temperatura de la superficie del mar y el color del océano con una resolución de 1/4 de grado. Tengo 4 polígonos que describen la Tengo 4 polígonos que describen el gran ecosistema marino de la plataforma continental del noreste de los Estados Unidos y sus subcomponentes que necesito utilizar para extraer los datos. Estoy trabajando con archivos compuestos mensuales de 1982-2014, por lo que esta rutina de extracción de datos debe ser automatizada. Los archivos ya están subconjuntados en la cuadrícula aproximada del área de trabajo de [35, 45, -80, -60].

Anteriormente, convertíamos los archivos de datos HDF5 en rásters en R y los procesábamos de esta manera, pero este método es realmente ineficiente y estoy seguro de que hay una solución mejor en Python utilizando los archivos NetCDF actuales.

Hasta ahora he utilizado GDAL y Fiona para leer los archivos shape y NetCDF4 para cargar los archivos de datos. No estoy seguro de cómo hacer el subconjunto de datos. He encontrado esto:

GDAL for Python: ¿extraer subdominios del archivo NetCDF?

Pero no tengo ni la más remota idea de cómo subdividir un archivo NetCDF utilizando algo que no sea un simple cuadro delimitador, lo que ciertamente no son estos polígonos.

Las rutinas de punto en polígono probablemente tardarían una eternidad en funcionar, pero tal vez podría subconjuntar los datos utilizando un cuadro delimitador más pequeño que se gira para ajustarse a estas formas como un punto de partida inicial y luego hacer una búsqueda de punto en polígono:

Subconjunto de un archivo netCDF curvilíneo (salida del modelo ROMS) utilizando un cuadro delimitador lon/lat.

¿Alguna idea?

EDITAR 1:

Acabo de encontrar el paquete OpenClimateGIS que parece que puede encajar perfectamente... Voy a tener un ir con esto para ver si puedo conseguir que funcione: http://ncpp.github.io/ocgis/examples.html#advanced-subsetting

Answer 1

Este puede ser adaptable a sus necesidades.

Si no te importa llamar a la línea de comandos desde python, podrías hacer algo como gdalwarp -cutline clip.shp -cl clip -crop_to_cutline input_raster output_raster_clipped.tif . -cwhere y -csql podría ser más apropiado opciones de gdalwarp para seleccionar uno de los cuatro polígonos para el recorte.

Answer 2

Mira esto: https://stackoverflow.com/questions/34585582/how-to-mask-the-specific-array-data-based-on-the-shapefile

Lo que debes tener en cuenta es que una vez que has cargado tu NetCDF, estás trabajando con un array de NumPy.

¿Qué quiere producir? ¿Estados resumidos basados en las áreas de los polígonos?

De todos modos, esto es lo que yo haría:

1. Cargue su shapefile y obtenga sus áreas en un formato compatible (apuntando al proceso de máscara matplotlib en el enlace anterior suena bien)
2. Cargue su archivo NetCDF y obtenga los datos en una única matriz numpy X, Y, T
3. Enmascarar esa matriz utilizando los polígonos (¿de uno en uno?)
4. Exporta tus estadísticas resumidas.

Answer 3

Puedes utilizar rioxarray. Aquí hay un ejemplo: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import rioxarray
import geopandas

geodf = geopandas.read_file(...)
xds = rioxarray.open_rasterio(...)
clipped = xds.rio.clip(geodf.geometry.apply(mapping), geodf.crs)