Comment obtenir GDAL pour créer des statistiques pour GTiff en Python

14

Je crée régulièrement mes propres rasters GeoTIFF avec GDAL en Python, par exemple:

from osgeo import gdal
from numpy import random
data = random.uniform(0, 10, (300, 200))
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create('MyRaster.tif', 200, 300)
band = ds.GetRasterBand(1)
band.WriteArray(data)
ds = band = None # save, close

cependant, lorsque le résultat est affiché avec ArcCatalog / ArcGIS, il apparaît en noir ou en gris, car il n'a pas de statistiques. Ceci est résolu soit en cliquant avec le bouton droit sur le raster et en choisissant "Calculer les statistiques ..." dans ArcCatalog (il existe plusieurs autres façons de le faire), soit en utilisant gdalinfo dans une invite de commande:

gdalinfo -stats MyRaster.tif

va générer MyRaster.tif.aux.xml, qui est utilisé par ArcGIS pour dimensionner correctement le raster. Le fichier PAM (Persistent Auxiliary Metadata) contient les statistiques, notamment les valeurs minimales et maximales:

<PAMDataset>
  <PAMRasterBand band="1">
    <Metadata>
      <MDI key="STATISTICS_MINIMUM">0</MDI>
      <MDI key="STATISTICS_MAXIMUM">10</MDI>
      <MDI key="STATISTICS_MEAN">5.0189833333333</MDI>
      <MDI key="STATISTICS_STDDEV">2.9131294111984</MDI>
    </Metadata>
  </PAMRasterBand>
</PAMDataset>

Ma question: existe-t-il un moyen intégré pour que GDAL crée un fichier de statistiques (autre que l'utilisation de la gdalinfo -statscommande)? Ou dois-je écrire le mien?

python gdal geotiff-tiff statistics Mike T
la source

Réponses:

13

Vous pouvez utiliser la méthode GetStatistics pour obtenir les statistiques.

par exemple.

stats =   ds.GetRasterBand(1).GetStatistics(0,1)

il reviendra (Min, Max, Mean, StdDev)

donc le xml peut être lu:

<PAMDataset>
  <PAMRasterBand band="1">
    <Metadata>
      <MDI key="STATISTICS_MINIMUM">stats[0]</MDI>
      <MDI key="STATISTICS_MAXIMUM">stats[1]</MDI>
      <MDI key="STATISTICS_MEAN">stats[2]</MDI>
      <MDI key="STATISTICS_STDDEV">stats[3]</MDI>
    </Metadata>
  </PAMRasterBand>
</PAMDataset>

Je ne connais aucun moyen pythonique de créer / manipuler un fichier xml, mais étant donné la nature simpliste du xml qui l'accompagne, il devrait être assez simple de le créer avec des opérations d'E / S sur fichier

nickves
la source
5
Il s'avère que band.GetStatistics(0,1)cela calculera réellement les statistiques et les ajoutera aux métadonnées GeoTIFF dans le fichier unique. Aucun autre fichier requis. Cependant, à partir des tests avec les produits Esri, il ne fonctionne qu'avec ArcGIS 10.0 et versions ultérieures, et non ArcGIS 9.3 ou antérieur.
Mike T
La fonction est décrite sur la page GDAL . Sur cette base, les deux arguments passés à la fonction sont bApproxOK (si les statistiques TRUE peuvent être calculées sur la base de vues d'ensemble ou d'un sous-ensemble de toutes les tuiles) et bForce (si les statistiques FALSE ne seront renvoyées que si cela peut être fait sans rescanner l'image) .
3

Si les statistiques sont déjà calculées et incluses dans le fichier en interne, gdalinfo -statsne créez pas de fichier de statistiques PAM supplémentaire (.aux.xml) pour utiliser GDAL 2.1.0. Mais il est très facile d'implémenter le .xml pour vous-même. Voici quelques modules Python intégrés expliqués pour faire ce genre de choses. Pour moi, j'ai utilisé l' API XML ElementTree avec le code ci-dessous:

import xml.etree.cElementTree as ET

stats = file.GetRasterBand(band).GetStatistics(0,1)

pamDataset = ET.Element("PAMDataset")
pamRasterband = ET.SubElement(pamDataset, "PAMRasterBand", band="1")
metadata = ET.SubElement(pamRasterband, "Metadata")
ET.SubElement(metadata, "MDI", key = "STATISTICS_MAXIMUM").text = str(stats[1])
ET.SubElement(metadata, "MDI", key = "STATISTICS_MEAN").text = str(stats[2])
ET.SubElement(metadata, "MDI", key = "STATISTICS_MINIMUM").text = str(stats[0])
ET.SubElement(metadata, "MDI", key = "STATISTICS_STDDEV").text = str(stats[3])

tree = ET.ElementTree(pamDataset)
tree.write(destFilePath + ".aux.xml")

Le résultat ressemble à:

<PAMDataset>
    <PAMRasterBand band="1">
        <Metadata>
            <MDI key="STATISTICS_MINIMUM">-40.65</MDI>
            <MDI key="STATISTICS_MEAN">10.2929293137</MDI>
            <MDI key="STATISTICS_MAXIMUM">45.050012207</MDI>
            <MDI key="STATISTICS_STDDEV">17.4892321447</MDI>
        </Metadata>
    </PAMRasterBand>
</PAMDataset>
Manuel
la source