Schreiben eines numpy-Arrays in eine Rasterdatei

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?

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:

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.

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

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

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)

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)