Extraiga datos de NetCDF utilizando un polígono de archivo de forma en Python [duplicado]

Necesito subconjuntar datos de un NetCDF usando un archivo de forma 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 grados. Tengo 4 polígonos que describen los Estados Unidos. El ecosistema marino grande de la plataforma continental del noreste y sus subcomponentes que necesito usar para extraer los datos. Estoy trabajando con archivos compuestos mensuales desde 1982 hasta 2014, por lo que esta rutina de extracción de datos debe ser automatizada. Los archivos ya están subconjuntos a la cuadrícula de área de trabajo aproximada de [35, 45, -80, -60].

Anteriormente, estábamos convirtiendo archivos de datos HDF5 en rásteres en R y procesándolos de esta manera, pero este método es realmente ineficiente y estoy seguro de que hay una mejor solución en Python usando los archivos NetCDF actuales.

Hasta ahora he estado usando GDAL y Fiona para leer en los archivos de forma y NetCDF4 para cargar los archivos de datos. No estoy seguro de cómo subconjugar los datos. Encontré esto:

GDAL para Python: ¿extraer subdominios del archivo NetCDF?

Pero no tengo la menor idea acerca de cómo subconjugar un archivo NetCDF utilizando algo más que un simple cuadro delimitador, que estos polígonos ciertamente no son.

El punto en las rutinas de polígonos probablemente tomaría una eternidad para funcionar, pero tal vez podría subconjugar los datos usando un cuadro delimitador más pequeño que se gira para ajustarse a estas formas como este como un punto de inicio inicial y luego hacer una búsqueda de punto en poli:

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

¿Algunas ideas?

EDITAR 1:

Acabo de encontrar el paquete OpenClimateGIS que parece que puede encajar perfectamente ... Voy a probar esto para ver si puedo hacerlo funcionar: http://ncpp.github.io/ocgis/examples. html # subconjunto avanzado

Esto podría ser adaptable a sus necesidades.

Si no le importa llamar a la línea de comando desde python, podría hacer algo como esto gdalwarp -cutline clip.shp -cl clip -crop_to_cutline input_raster output_raster_clipped.tif. -cwherey -csqlpodrían ser opciones de gdalwarp más apropiadas para seleccionar uno de los cuatro polígonos para el recorte.

Eche un vistazo a esto: /programming/34585582/how-to-mask-the-specific-array-data-based-on-the-shapefile

Lo que debe tener en cuenta es que una vez que haya cargado su NetCDF, estará trabajando con una matriz NumPy.

¿Qué estás tratando de generar? Resumen de estadísticas basadas en las áreas de polígono?

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

1. Cargue su archivo de formas y obtenga sus áreas en un formato compatible (con el objetivo del proceso de máscara matplotlib en el enlace anterior suena bien)
2. Cargue su archivo NetCDF y obtenga los datos en una sola matriz numpy X, Y, T
3. Enmascara esa matriz usando los polígonos (¿uno a la vez?)
4. Exporta tus estadísticas de resumen.

Puedes usar 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)