Extraire des données de NetCDF en utilisant un polygone de fichier de formes en Python [duplicate]

J'ai besoin de sous-définir les données d'un NetCDF à l'aide d'un fichier de formes spécifique. Les données sont la température de la surface de la mer et la couleur de l'océan à une résolution de 1/4 degré. J'ai 4 polygones décrivant les États-Unis. Grand écosystème marin du plateau continental nord-est et ses sous-composants que je dois utiliser pour extraire les données. Je travaille avec des fichiers composites mensuels de 1982 à 2014, donc cette routine d'extraction de données doit être automatisée. Les fichiers sont déjà classés dans la grille de la zone de travail approximative de [35, 45, -80, -60].

Auparavant, nous convertissions les fichiers de données HDF5 en rasters dans R et les traitions de cette façon, mais cette méthode est vraiment inefficace et je suis sûr qu'il existe une meilleure solution en Python en utilisant les fichiers NetCDF actuels.

Jusqu'à présent, j'ai utilisé GDAL et Fiona pour lire les fichiers de formes et NetCDF4 pour charger les fichiers de données. Je ne sais pas comment procéder pour sous-définir les données. J'ai trouvé ça:

GDAL pour Python: extraire des sous-domaines du fichier NetCDF?

Mais je n'ai pas la moindre idée de la façon de sous-définir un fichier NetCDF en utilisant autre chose qu'une simple boîte englobante, ce que ces polygones ne sont certainement pas.

Le point dans les routines polygonales prendrait probablement une éternité pour fonctionner, mais je pourrais peut-être sous-définir les données en utilisant un cadre de délimitation plus petit qui est tourné pour s'adapter à ces formes comme ceci comme point de départ initial, puis faire une recherche point par point:

Sous-ensemble d'un fichier netCDF curviligne (sortie du modèle ROMS) à l'aide d'un cadre de délimitation lon / lat.

Des idées?

EDIT 1:

Je viens de découvrir le package OpenClimateGIS qui semble pouvoir parfaitement convenir ... J'essaierai avec cela pour voir si je peux le faire fonctionner: http://ncpp.github.io/ocgis/examples. html # advanced-subsetting

Cela pourrait être adaptable à vos besoins.

Si cela ne vous dérange pas d'appeler la ligne de commande à partir de python, vous pouvez faire quelque chose comme gdalwarp -cutline clip.shp -cl clip -crop_to_cutline input_raster output_raster_clipped.tif. -cwhereet -csqlpourrait être des options gdalwarp plus appropriées pour sélectionner l'un des quatre polygones à découper.

Jetez un œil à ceci: /programming/34585582/how-to-mask-the-specific-array-data-based-on-the-shapefile

Ce que vous voulez garder à l'esprit, c'est qu'une fois que vous avez chargé votre NetCDF, vous travaillez avec un tableau NumPy.

Qu'essayez-vous de produire? Statistiques récapitulatives basées sur les zones de polygones?

Quoi qu'il en soit, voici ce que je ferais:

1. Chargez votre fichier de formes et mettez vos zones dans un format pris en charge (viser le processus de masque matplotlib dans le lien ci-dessus sonne bien)
2. Chargez votre fichier NetCDF et récupérez les données dans un seul tableau numpy X, Y, T
3. Masquez ce tableau en utilisant les polygones (un à la fois?)
4. Exportez vos statistiques récapitulatives.

Vous pouvez utiliser rioxarray. Voici un exemple: 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)