Suite

Extraction de la proportion en pixels des classes d'un raster plus fin dans un raster plus grossier

Extraction de la proportion en pixels des classes d'un raster plus fin dans un raster plus grossier


J'ai deux rasters, l'un a une résolution de 1 km. L'autre a une résolution spatiale de 30 m. Ce dernier pixel est une classification. Je souhaite superposer le raster de 1 km sur le raster de classification et extraire pour chaque pixel de 1 km la distribution des classes correspondant à la classification de 30 m.

La classification 30 m comprend 7 classes donc je voudrais que ce processus produise 7 nouvelles couches à 1 km. Un pour chaque classe. A chaque pixel de 1 km sera désormais associé un % correspondant de présence de chaque classe présente dans le classement 30 m. Si l'on faisait un tableau où chaque ligne est un pixel, on aurait un tableau comme celui-ci :

Une difficulté supplémentaire est que les rasters sont assez gros. J'aimerais vraiment y parvenir en utilisant python.

ÉDITER:

J'ai ajouté des sous-ensembles de mes deux ensembles de données à ce lien (à des résolutions de 30 m et 1000 m respectivement), au cas où quelqu'un voudrait proposer une preuve de concept; serait génial.


Étape 1

Créez des rasters binaires pour chacune des classes uniques. Il peut s'agir de rasters à 1 bande pour chaque classe ou d'un seul raster avec une bande pour chaque classe (par exemple, GeoTIFF). Si vous utilisez GTiff, vous pouvez utiliser l'option de créationNBITS=1pour économiser de l'espace. Vous pouvez également envisager des rasters à deux bits pour stocker une logique à trois valeurs où le troisième (par exemple 2) est NODATA, ce qui nécessiteraitNBITS=2.

Pour chaque bande/classe, les pixels seront 0 là où cette classe n'existe pas, ou 1 là où cette classe existe.

import numpy as np from osgeo import gdal gdal.UseExceptions() # Récupère les données du raster avec les classifications ds = gdal.Open('cropped_30m.tif') band = ds.GetRasterBand(1) class_ar = band.ReadAsArray() gt = ds .GetGeoTransform() pj = ds.GetProjection() ds = band = None # close # Définit les valeurs raster pour chaque classe, à relier à chaque bande class_ids = (np.arange(31) + 1).tolist() # Make un nouveau bit rasters drv = gdal.GetDriverByName('GTiff') ds = drv.Create('bit_raster.tif', class_ar.shape[1], class_ar.shape[0], len(class_ids), gdal.GDT_Byte, [ 'NBITS=1', 'COMPRESS=LZW', 'INTERLEAVE=BAND']) ds.SetGeoTransform(gt) ds.SetProjection(pj) pour bidx dans range(ds.RasterCount) : bande = ds.GetRasterBand(bidx + 1) ) # créer une sélection booléenne = (class_ar == class_ids[bidx]) band.WriteArray(selection.astype('B')) ds = band = None # enregistrer, fermer

Par exemple. pour Rouge = Bande 10, Vert = Bande 13 et Bleu = Bande 27, on peut visualiser :

Le noir indique qu'il s'agit d'une autre classe/bande.

Étape 2

Utilisez gdalwarp pour obtenir le moyenne des valeurs 0 et 1 à partir d'une grille d'échantillonnage plus grossière. Cette méthode nécessite GDAL >= 1.10.0. Un moyen très simple de le faire est de spécifier une résolution cible sur la ligne de commande, par ex. 1km :

$ gdalwarp -ot Float32 -tr 1000 1000 -r moyenne bit_raster.tif moyenne_1km.tif

Le résultat aura autant de bandes que de classes, et chaque bande décrira la fraction de cette classe pour chaque pixel, entre 0 et 1. Si vous préférez le pourcentage, multipliez-le par 100 (par exemple, voir gdal_calc.py).

Il existe également une interface Python vers gdalwarp (par exemple, cette réponse), qui peut utiliser un modèle raster pour la forme et où NODATA doit être utilisé. Par example:

# Ouvrir le raster de l'étape 1 src_ds = gdal.Open('bit_raster.tif') # Ouvrir un modèle ou copier un tableau, pour les dimensions et le masque NODATA cpy_ds = gdal.Open('cropped_1000m.tif') band = cpy_ds.GetRasterBand(1 ) cpy_mask = (band.ReadAsArray() == band.GetNoDataValue()) # Raster de résultat, avec la même résolution et la même position que la copie raster dst_ds = drv.Create('average2_1km.tif', cpy_ds.RasterXSize, cpy_ds.RasterYSize, len(class_ids), gdal.GDT_Float32, ['INTERLEAVE=BAND']) dst_ds.SetGeoTransform(cpy_ds.GetGeoTransform()) dst_ds.SetProjection(cpy_ds.GetProjection()) # Faites la même chose que gdalwarp -r average ; cela peut prendre un certain temps pour terminer gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Average) # Convertir toutes les fractions en pourcentage et appliquer le même masque NODATA à partir du raster de copie NODATA = -9999 pour bidx dans la plage (dst_ds .RasterCount): band = dst_ds.GetRasterBand(bidx + 1) ar = band.ReadAsArray() * 100.0 ar[cpy_mask] = NODATA band.WriteArray(ar) band.SetNoDataValue(NODATA) # Enregistrer et fermer tous les rasters src_ds = cpy_ds = dst_ds = bande = Aucun

Par exemple. les trois mêmes bandes (classes) peuvent être visualisées en RVB :

Ou tirez-le dans QGIS et utilisez l'outil Identifier les entités pour inspecter la fraction pour chaque bande/classification à un seul endroit :

Ces nombres totalisent 100%, ce qui est un bon contrôle.


Voir la vidéo: Calcul statistique dans un raster, surface et pourcentage. statistic in raster, area, pourcent