Jak zmusić GDAL do tworzenia statystyk dla GTiff w Pythonie

14

Regularnie tworzę własne rastry GeoTIFF z GDAL w Pythonie, np .:

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

Jednak gdy wynik jest wyświetlany za pomocą ArcCatalog / ArcGIS, wygląda na czarny lub szary, ponieważ nie ma statystyk. Można to rozwiązać albo klikając prawym przyciskiem myszy raster i wybierając „Oblicz statystyki ...” w ArcCatalog (istnieje kilka innych sposobów, aby to zrobić), lub używając gdalinfo w wierszu polecenia:

gdalinfo -stats MyRaster.tif

wygeneruje MyRaster.tif.aux.xml, który jest używany przez ArcGIS do prawidłowego skalowania rastra. Plik PAM (Persistent Auxiliary Metadata) zawiera statystyki, w szczególności wartości minimalne i maksymalne:

<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>

Moje pytanie: czy istnieje wbudowany sposób zmuszenia GDAL do utworzenia pliku statystyk (innego niż użycie gdalinfo -statspolecenia)? Czy też muszę pisać własne?

— Mike T.
źródło

13

Aby uzyskać statystyki, możesz użyć metody GetStatistics.

na przykład.

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

zwróci (Min, Max, Mean, StdDev)

więc xml można odczytać:

<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>

Nie znam żadnego pythonicznego sposobu tworzenia / manipulowania plikiem xml, ale biorąc pod uwagę uproszczony charakter towarzyszącego mu pliku xml, tworzenie go za pomocą operacji we / wy pliku powinno być dość proste

— nickves
źródło

5

Okazuje się, że band.GetStatistics(0,1)faktycznie obliczy statystyki i doda je do metadanych GeoTIFF w jednym pliku. Żadne inne pliki nie są wymagane. Jednak z testów produktów Esri działa tylko z ArcGIS 10.0 i nowszymi, a nie ArcGIS 9.3 lub wcześniejszymi.

— Mike T

Funkcja jest opisana na stronie GDAL . Na tej podstawie dwa argumenty przekazane do funkcji to bApproxOK (jeśli statystyki PRAWDA mogą być obliczone na podstawie przeglądów lub podzbioru wszystkich kafelków) i bForce (jeśli statystyki FAŁSZ zostaną zwrócone tylko wtedy, gdy można to zrobić bez ponownego skanowania obrazu) .

3

Jeśli statystyki są już obliczone i uwzględnione w pliku wewnętrznie, gdalinfo -statsnie należy tworzyć dodatkowego pliku statystyk PAM (.aux.xml) do korzystania z GDAL 2.1.0. Ale bardzo łatwo jest zaimplementować plik .xml dla własnego. Oto kilka wbudowanych modułów Pythona, które wyjaśniają, jak to zrobić. Dla siebie użyłem API ElementTree XML z poniższym kodem:

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")

Wynik wygląda następująco:

<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
źródło