Écriture d'un tableau numpy dans un fichier raster

30

Je suis nouveau au SIG.

J'ai un code qui convertit les images infrarouges de Mars en cartes d'inertie thermique, qui sont ensuite stockées sous forme de tableaux numpy 2D. J'ai enregistré ces cartes en tant que fichiers hdf5 mais j'aimerais vraiment les enregistrer en tant qu'images raster afin de pouvoir les traiter dans QGIS. J'ai parcouru plusieurs recherches pour trouver comment le faire, mais sans succès. J'ai essayé de suivre les instructions du didacticiel à http://www.gis.usu.edu/~chrisg/python/ mais les fichiers que je produis en utilisant son exemple de code s'ouvrent sous forme de cases grises lorsque je les importe dans QGIS. J'ai l'impression que si quelqu'un pouvait suggérer la procédure la plus simple possible à un exemple simplifié de ce que j'aimerais faire, je pourrais peut-être progresser. J'ai QGIS et GDAL, je serais très heureux d'installer d'autres frameworks que n'importe qui pourrait recommander. J'utilise Mac OS 10.7.

Donc, par exemple, si j'ai un tableau numpy d'inertie thermique qui ressemble à:

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

Et pour chaque pixel j'ai la latitude et la longitude:

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

Quelle procédure les gens recommandent-ils pour convertir ces données en un fichier raster que je peux ouvrir dans QGIS?

EddyTheB
la source
À quelle diapositive du didacticiel faites-vous référence?
RK

Réponses:

23

Une solution possible à votre problème: convertissez-le en un raster ASCII, dont la documentation est ici . Cela devrait être assez facile à faire avec python.

Donc, avec vos données d'exemple ci-dessus, vous vous retrouveriez avec ce qui suit dans un fichier .asc:

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

Cela s'ajoute avec succès à QGIS et ArcGIS, et stylisé dans ArcGIS, il ressemble à ceci: version raster ci-dessus

Addendum: Bien que vous puissiez l'ajouter à QGIS comme indiqué, si vous essayez d'aller dans ses propriétés (pour le styliser), QGIS 1.8.0 se bloque. Je suis sur le point de signaler cela comme un bug. Si cela vous arrive aussi, il existe de nombreux autres SIG gratuits.

SIG-Jonathan
la source
C'est fantastique, merci. Et j'imagine qu'après avoir écrit mon tableau sous forme de fichier ascii, je pourrais le convertir en un format binaire en utilisant une fonction de conversion pré-écrite.
EddyTheB
Pour info, je n'ai pas eu le problème de blocage avec QGIS, j'ai aussi la version 1.8.0.
EddyTheB
31

Voici un exemple que j'ai écrit pour un atelier qui utilise les modules Python numpy et gdal. Il lit les données d'un fichier .tif dans un tableau numpy, reclasse les valeurs du tableau, puis les réécrit dans un fichier .tif.

D'après votre explication, il semble que vous ayez réussi à écrire un fichier valide, mais il vous suffit de le symboliser dans QGIS. Si je me souviens bien, lorsque vous ajoutez un raster pour la première fois, il affiche souvent toutes une couleur si vous n'avez pas de carte de couleurs préexistante.

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
la source
1
+1 pour le rinçage - me cognais la tête contre le mur en essayant de comprendre comment «sauver» la chose!
badgley
J'ai dû ajouter outDs = Nonepour le sauver
JaakL
23

J'ai finalement trouvé cette solution, que j'ai tirée de cette discussion ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Je l'aime parce que je peux passer directement d'un tableau numpy à un fichier raster tif. Je serais très reconnaissant pour les commentaires qui pourraient améliorer la solution. Je le posterai ici au cas où quelqu'un d'autre chercherait une réponse similaire.

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
la source
3
La "rotation est en deux fois" pour tenir compte de l'effet d'un bit de rotation de y sur x et du bit de rotation de x sur y. Voir lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html qui essaie d'expliquer les interrelations entre les paramètres de "rotation".
Dave X
Ce message est vraiment utile, merci. Dans mon cas, cependant, j'obtiens un fichier tif qui est complètement noir lorsque je l'ouvre en tant qu'image en dehors d'ArcGIS. Ma référence spatiale est le British National Grid (EPSG = 27700), et les unités sont en mètres.
FaCoffee
J'ai posté une question ici: gis.stackexchange.com/questions/232301/…
FaCoffee
Avez-vous découvert comment définir l'encodage IAU2000: 49900 Mars?
K.-Michael Aye
4

Il y a aussi une belle solution dans le livre de recettes officiel GDAL / OGR pour Python.

Cette recette crée un raster à partir d'un tableau

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
la source
Cette recette est bonne mais il y a un problème avec le fichier tiff final. Les valeurs lat-lon des pixels ne sont pas correctes.
Shubham_geo
Vous pourriez voir des incompatibilités étranges entre ESRI WKT et OGC WKT: gis.stackexchange.com/questions/129764/…
Adam Erickson
Une chose que j'ai rencontrée est que la façon dont vous avez mentionné changera certainement le tableau en raster avec facilité. Mais nous devons géoréférencer ce raster avec le coin supérieur gauche et le coin inférieur droit coordonnés à l'aide de gdal_translate. Une façon de le faire est de suivre les deux étapes: 1) Trouvez d'abord les valeurs lat-lon en haut à gauche et en bas à droite via gdalinfo 2) Ensuite, via gdal_translate, utilisez le géotiff (généré avec l'approche mentionnée ci-dessus de conversion de tableau en raster) pour le géoréférencer avec les coordonnées lat-lon en haut à gauche et en bas à droite.
Shubham_geo
0

Une alternative à l'approche suggérée dans les autres réponses est d'utiliser le rasteriopackage. J'ai eu des problèmes pour les générer gdalet j'ai trouvé ce site utile.

En supposant que vous avez un autre fichier tif ( other_file.tif) et un tableau numpy ( numpy_array) qui ont la même résolution et l'étendue que ce fichier, c'est l'approche qui a fonctionné pour moi:

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