Schreiben eines numpy-Arrays in eine Rasterdatei

30

Ich bin neu bei GIS.

Ich habe einen Code, der Infrarotbilder des Mars in thermische Trägheitskarten umwandelt, die dann als 2D-Numpy-Arrays gespeichert werden. Ich habe diese Karten als hdf5-Dateien gespeichert, möchte sie aber wirklich als Rasterbilder speichern, damit ich sie in QGIS verarbeiten kann. Ich habe mehrere Suchanfragen durchgeführt, um herauszufinden, wie das geht, aber ohne Glück. Ich habe versucht, den Anweisungen im Tutorial unter http://www.gis.usu.edu/~chrisg/python/ zu folgen, aber die Dateien, die ich mit seinem Beispielcode erzeuge, werden beim Importieren in QGIS als graue Kästchen geöffnet. Ich habe das Gefühl, wenn jemand einem vereinfachten Beispiel für das, was ich tun möchte, die einfachste Vorgehensweise vorschlagen könnte, könnte ich vielleicht Fortschritte erzielen. Ich habe QGIS und GDAL, ich würde mich sehr freuen, andere Frameworks zu installieren, die jeder empfehlen kann. Ich benutze Mac OS 10.7.

Wenn ich zum Beispiel eine Reihe von thermischen Trägheiten habe, die wie folgt aussehen:

TI = ( (0.1, 0.2, 0.3, 0.4),
       (0.2, 0.3, 0.4, 0.5),
       (0.3, 0.4, 0.5, 0.6),
       (0.4, 0.5, 0.6, 0.7) )

Und für jedes Pixel habe ich den Breiten- und Längengrad:

lat = ( (10.0, 10.0, 10.0, 10.0),
        ( 9.5,  9.5,  9.5,  9.5),
        ( 9.0,  9.0,  9.0,  9.0),
        ( 8.5,  8.5,  8.5,  8.5) )
lon = ( (20.0, 20.5, 21.0, 21.5),
        (20.0, 20.5, 21.0, 21.5),
        (20.0, 20.5, 21.0, 21.5),
        (20.0, 20.5, 21.0, 21.5) ) 

Welches Verfahren wird empfohlen, um diese Daten in eine Rasterdatei zu konvertieren, die ich in QGIS öffnen kann?

EddyTheB
quelle
Auf welche Folie im Tutorial beziehen Sie sich?
RK

Antworten:

23

Eine mögliche Lösung für Ihr Problem: Konvertieren Sie es in ein ASCII-Raster, dessen Dokumentation Sie hier finden . Dies sollte mit Python relativ einfach zu bewerkstelligen sein.

Wenn Sie die obigen Beispieldaten verwenden, erhalten Sie in einer .asc-Datei Folgendes:

ncols 4
nrows 4
xllcorner 20
yllcorner 8.5
cellsize 0.5
nodata_value -9999
0.1 0.2 0.3 0.4
0.2 0.3 0.4 0.5
0.3 0.4 0.5 0.6
0.4 0.5 0.6 0.7

Dies fügt QGIS und ArcGIS erfolgreich hinzu und sieht in ArcGIS folgendermaßen aus: Raster-Version von oben

Nachtrag: Während Sie es wie beschrieben zu QGIS hinzufügen können, bleibt QGIS 1.8.0 hängen, wenn Sie versuchen, die Eigenschaften dafür aufzurufen (um es zu stilisieren). Ich werde das als Fehler melden. Wenn Ihnen das auch passiert, gibt es viele andere kostenlose GIS.

GIS-Jonathan
quelle
Das ist fantastisch, danke. Und ich stelle mir vor, dass ich mein Array, nachdem ich es als ASCII-Datei geschrieben habe, mit einer vorab geschriebenen Konvertierungsfunktion in ein Binärformat konvertieren könnte.
EddyTheB
Zu Ihrer Information, ich hatte keine Probleme mit QGIS, ich habe auch die Version 1.8.0.
EddyTheB
31

Unten ist ein Beispiel, das ich für einen Workshop geschrieben habe, in dem die Python-Module numpy und gdal verwendet werden. Es liest Daten aus einer .tif-Datei in ein Numpy-Array, führt eine Neuklassifizierung der Werte im Array durch und schreibt sie dann wieder in ein .tif-Array.

Nach Ihrer Erklärung klingt es so, als ob es Ihnen gelungen ist, eine gültige Datei zu erstellen, Sie müssen sie jedoch nur in QGIS symbolisieren. Wenn ich mich richtig erinnere, werden beim erstmaligen Hinzufügen eines Rasters häufig alle Farben angezeigt, wenn Sie keine bereits vorhandene Farbzuordnung haben.

import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *


# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
if inDs is None:
  print 'Could not open image file'
  sys.exit(1)

# read in the crop data and get info about it
band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize

cropData = band1.ReadAsArray(0,0,cols,rows)

listAg = [1,5,6,22,23,24,41,42,28,37]
listNotAg = [111,195,141,181,121,122,190,62]

# create the output image
driver = inDs.GetDriver()
#print driver
outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
if outDs is None:
    print 'Could not create reclass_40.tif'
    sys.exit(1)

outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.int16)


for i in range(0, rows):
    for j in range(0, cols):

    if cropData[i,j] in listAg:
        outData[i,j] = 100
    elif cropData[i,j] in listNotAg:
        outData[i,j] = -100
    else:
        outData[i,j] = 0


# write the data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(-99)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

del outData
DavidF
quelle
1
+1 für das Spülen - habe meinen Kopf gegen die Wand geschlagen, um herauszufinden, wie ich das Ding 'retten' kann!
Badgley
Ich musste hinzufügen outDs = None, um es zu retten
JaakL
23

Ich bin endlich auf diese Lösung gestoßen, die ich aus dieser Diskussion gewonnen habe ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Ich mag es, weil ich direkt von einem Numpy-Array zu einer TIF-Rasterdatei wechseln kann. Ich wäre sehr dankbar für Kommentare, die die Lösung verbessern könnten. Ich werde es hier posten, falls jemand nach einer ähnlichen Antwort sucht.

import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

array = np.array(( (0.1, 0.2, 0.3, 0.4),
                   (0.2, 0.3, 0.4, 0.5),
                   (0.3, 0.4, 0.5, 0.6),
                   (0.4, 0.5, 0.6, 0.7),
                   (0.5, 0.6, 0.7, 0.8) ))
# My image array      
lat = np.array(( (10.0, 10.0, 10.0, 10.0),
                 ( 9.5,  9.5,  9.5,  9.5),
                 ( 9.0,  9.0,  9.0,  9.0),
                 ( 8.5,  8.5,  8.5,  8.5),
                 ( 8.0,  8.0,  8.0,  8.0) ))
lon = np.array(( (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5) ))
# For each pixel I know it's latitude and longitude.
# As you'll see below you only really need the coordinates of
# one corner, and the resolution of the file.

xmin,ymin,xmax,ymax = [lon.min(),lat.min(),lon.max(),lat.max()]
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)   
# That's (top left x, w-e pixel resolution, rotation (0 if North is up), 
#         top left y, rotation (0 if North is up), n-s pixel resolution)
# I don't know why rotation is in twice???

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.
                                             # Anyone know how to specify the 
                                             # IAU2000:49900 Mars encoding?
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system 
                                                   # to the file
output_raster.GetRasterBand(1).WriteArray(array)   # Writes my array to the raster

output_raster.FlushCache()
EddyTheB
quelle
3
Die "Drehung erfolgt zweimal", um die Auswirkung eines gedrehten Bits von y auf x und des gedrehten Bits von x auf y zu berücksichtigen. Siehe lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html, in dem versucht wird, die Wechselbeziehungen zwischen den "Rotations" -Parametern zu erklären.
Dave X
Dieser Beitrag ist wirklich nützlich, danke. In meinem Fall wird jedoch eine TIF-Datei angezeigt, die vollständig schwarz ist, wenn ich sie als Bild außerhalb von ArcGIS öffne. Mein räumlicher Bezug ist das British National Grid (EPSG = 27700), und die Einheiten sind Meter.
FaCoffee
Ich habe hier eine Frage gestellt: gis.stackexchange.com/questions/232301/…
FaCoffee
Haben Sie herausgefunden, wie die Mars-Codierung für IAU2000: 49900 festgelegt wird?
K.-Michael Aye
4

Es gibt auch eine schöne Lösung im offiziellen GDAL / OGR-Kochbuch für Python.

Dieses Rezept erstellt ein Raster aus einem Array

import gdal, ogr, os, osr
import numpy as np


def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):

    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    reversed_arr = array[::-1] # reverse array so the tif looks like the array
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster


if __name__ == "__main__":
    rasterOrigin = (-123.25745,45.43013)
    pixelWidth = 10
    pixelHeight = 10
    newRasterfn = 'test.tif'
    array = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])


    main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)
Adam Erickson
quelle
Dieses Rezept ist gut, aber es gibt ein Problem mit der endgültigen TIFF-Datei. Die Lat-Lon-Werte der Pixel sind nicht korrekt.
Shubham_geo
Möglicherweise treten
Adam Erickson,
Eine Sache, auf die ich gestoßen bin, ist, dass die von Ihnen erwähnte Art und Weise das Array mit Leichtigkeit in ein Raster verwandelt. Aber wir müssen dieses Raster mit Hilfe von gdal_translate links oben und rechts unten georeferenzieren. Eine Möglichkeit, dies zu tun, sind die beiden folgenden Schritte: 1) Ermitteln Sie zuerst die oberen linken und unteren rechten Lat-Lon-Werte über gdalinfo. um es mit den Lat-Lon-Koordinaten oben links und unten rechts zu georeferenzieren.
Shubham_geo
0

Eine Alternative zu dem in den anderen Antworten vorgeschlagenen Ansatz ist die Verwendung des rasterioPakets. Ich hatte Probleme beim Generieren dieser mit gdalund fand diese Site nützlich.

Angenommen, Sie haben eine andere tif-Datei ( other_file.tif) und ein numpy-Array ( numpy_array) mit der gleichen Auflösung und Ausdehnung wie diese Datei. Dies ist der Ansatz, der für mich funktioniert hat:

import rasterio as rio    

with rio.open('other_file.tif') as src:
    ras_data = src.read()
    ras_meta = src.profile

# make any necessary changes to raster properties, e.g.:
ras_meta['dtype'] = "int32"
ras_meta['nodata'] = -99

with rio.open('outname.tif', 'w', **ras_meta) as dst:
    dst.write(numpy_array, 1)
Tim Williams
quelle