Suite

Transformer le shapefile en masque et calculer la moyenne


J'utilise Python pour traiter d'énormes tableaux stockés sous forme de fichiers NetCDF. Je voudrais calculer la moyenne d'une zone définie par un fichier de formes. Je viens d'installer GDAL mais s'il y a d'autres outils que je devrais utiliser, merci de me le faire savoir.

Ma question principale est donc de savoir comment puis-je transformer mon fichier de formes en masque?

Ensuite, utiliser mon masque pour calculer la moyenne NumPy de mon tableau dans la zone du masque ?

shp="country.shp" array=mynumpyarray mask=useGDALtochangeshapefiletomask Meanofarea=N.mean(array,limits=mask)

J'ai pu convertir mon fichier de formes en raster à l'aide du lien gdal_RasterizeLayer fourni, puis d'un tableau masqué à l'aide

maskarray=mask_ds.GetRasterBand(1).ReadAsArray()

Maintenant, il ne s'agit plus que de faire correspondre mon étendue NetCDF avec mon raster/maskarray. Voici mes étendues comme demandé.

mon étendue shapefile/raster : x_min, x_max, y_min, y_max 140.962408758 149.974994992 -39.1366533667 -33.9813898583 mon étendue de fichier netcdf : min longitude : 139,8 longitude max : 150,0 min latitude -39,2 latitude max : -33,6 taille LAT 106 taille LON 193

Vous recherchez la fonction gdal.RasterizeLayer.

Vous pouvez ensuite utiliser ReadAsArray pour transformer le polygone rastérisé en un tableau numpy.

En fonction de l'étendue de votre fichier NetCDF et des lignes/colonnes, le code suivant devrait vous générer un masque numpy 0-1 qui correspond exactement au NetCDF.

shapefile=r'quel que soit le chemin de votre fichier de formes' xmin,ymin,xmax,ymax=[139.8,-39.2,150.0,-33.6] #Vos étendues comme indiqué ci-dessus ncols,nrows=[193,106] #Vos lignes/cols comme indiqué ci-dessus maskvalue = 1 xres=(xmax-xmin)/float(ncols) yres=(ymax-ymin)/float(nrows) geotransform=(xmin,xres,0,ymax,0, -yres) src_ds = ogr.Open(shapefile) ) src_lyr=src_ds.GetLayer() dst_ds = gdal.GetDriverByName('MEM').Create(", ncols, nrows, 1 ,gdal.GDT_Byte) dst_rb = dst_ds.GetRasterBand(1) dst_rb.Fill(0) #initialiser raster(0) avec des zéros dst_rb.SetNoDataValue(0) dst_ds.SetGeoTransform(geotransform) err = gdal.RasterizeLayer(dst_ds, [1], src_lyr, burn_values=[maskvalue]) dst_ds.FlushCache(1) mask_arr=dst_Baster.A.G )

La solution la plus simple consiste peut-être à appeler gdal_rasterize à partir de votre code Python (utilisezsous-processus.appelou alorssubprocess.check_call), soit en utilisant votre fichier NetCDF comme destination, soit en créant un fichier image séparé (GeoTIFF est toujours un bon pari) et en le chargeant dans un tableau numpy. Il est peut-être possible d'utiliser le format In Memory Raster de GDAL, mais je ne suis pas sûr de la prise en charge de Python et de la manière dont il jouerait avec les tableaux numpy.

Il y a une surcharge avec l'appel d'un sous-processus et le chargement d'une image, mais cela sera probablement amorti par le temps qu'il faudra pour faire la rastérisation réelle.

Si vous voulez lancer le vôtre, vous pouvez écrire un simple moteur de rendu de lignes de balayage, mais je soupçonne que ce serait plus lourd et peut-être plus lent s'il était écrit en Python pur.


Maintenant en 2017, j'utiliserasterio.mask.masketfionapour pixelliser le fichier de forme, puis l'appliquer au jeu de données netCDF à l'aidexarray.


Maintenant en 2019, il y a https://github.com/corteva/rioxarray pour faciliter l'intégration derasterioetxarray. Le découpage d'un ensemble de données (par exemple à partir de NetCDF) avec une géométrie est répertorié à titre d'exemple ici : https://corteva.github.io/rioxarray/stable/examples/clip_geom.html


Voir la vidéo: Convertir des polygones en lignes (Octobre 2021).