Escribir matriz numpy en archivo ráster

30

Soy nuevo en SIG.

Tengo un código que convierte las imágenes infrarrojas de Marte en mapas de inercia térmica, que luego se almacenan como matrices numpy 2D. He guardado estos mapas como archivos hdf5, pero realmente me gustaría guardarlos como imágenes ráster para poder procesarlos en QGIS. He realizado varias búsquedas para encontrar cómo hacerlo, pero sin suerte. Intenté seguir las instrucciones del tutorial en http://www.gis.usu.edu/~chrisg/python/ pero los archivos que produzco usando su código de ejemplo se abren como cuadros grises cuando los importo a QGIS. Siento que si alguien pudiera sugerir el procedimiento más simple posible a un ejemplo simplificado de lo que me gustaría hacer, entonces podría progresar un poco. Tengo QGIS y GDAL, me encantaría instalar otros marcos que cualquiera podría recomendar. Yo uso Mac OS 10.7.

Entonces, si, por ejemplo, tengo una gran variedad de inercia térmica que se parece a:

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

Y para cada píxel tengo la latitud y longitud:

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

¿Qué procedimiento recomendaría la gente para convertir estos datos en un archivo ráster que pueda abrir en QGIS?

EddyTheB
fuente
¿A qué diapositiva del tutorial te refieres?
RK

Respuestas:

23

Una posible solución a su problema: conviértalo en un ráster ASCII, cuya documentación está aquí . Esto debería ser bastante fácil de hacer con Python.

Entonces, con sus datos de ejemplo anteriores, terminaría con lo siguiente en un archivo .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

Esto se agrega con éxito a QGIS y ArcGIS, y estilizado en ArcGIS se ve así: versión ráster de arriba

Anexo: Si bien puede agregarlo a QGIS como se indicó, si intenta acceder a las propiedades (para estilizarlo), QGIS 1.8.0 se bloquea. Estoy a punto de reportar eso como un error. Si esto te sucede a ti también, entonces hay muchos otros SIG gratuitos disponibles.

GIS-Jonathan
fuente
Eso es fantástico, gracias. E imagino que habiendo escrito mi matriz como un archivo ascii podría convertirlo a un formato binario usando una función de conversión preescrita.
EddyTheB
FYI, no tuve el problema pendiente con QGIS, también tengo la versión 1.8.0.
EddyTheB
31

A continuación hay un ejemplo que escribí para un taller que utiliza los módulos Python numpy y gdal. Lee datos de un archivo .tif en una matriz numpy, hace una reclasificación de los valores en la matriz y luego los vuelve a escribir en un .tif.

Según su explicación, parece que podría haber logrado escribir un archivo válido, pero solo necesita simbolizarlo en QGIS. Si no recuerdo mal, cuando agregas un ráster por primera vez, a menudo muestra todos los colores si no tienes un mapa de colores preexistente.

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
fuente
1
¡+1 por sonrojarse - estaba golpeando mi cabeza contra la pared tratando de descubrir cómo 'salvar' la cosa!
badgley
Tuve que agregar outDs = Nonepara guardarlo
JaakL
23

Finalmente encontré esta solución, que obtuve de esta discusión ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Me gusta porque puedo pasar directamente de una matriz numpy a un archivo ráster tif. Estaría muy agradecido por los comentarios que podrían mejorar la solución. Lo publicaré aquí en caso de que alguien más busque una respuesta similar.

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
fuente
3
La "rotación es dos veces" para tener en cuenta el efecto de un bit girado de y sobre x y el bit girado de x sobre y. Ver lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html que intenta explicar las interrelaciones entre los parámetros de "rotación".
Dave X
Esta publicación es realmente útil, gracias. En mi caso, sin embargo, obtengo un archivo tif que está completamente negro cuando lo abro como una imagen fuera de ArcGIS. Mi referencia espacial es la Red Nacional Británica (EPSG = 27700), y las unidades son metros.
FaCoffee
He publicado una pregunta aquí: gis.stackexchange.com/questions/232301/…
FaCoffee
¿Descubrió cómo configurar la codificación IAU2000: 49900 Mars?
K.-Michael Aye
4

También hay una buena solución en el libro de cocina oficial de GDAL / OGR para Python.

Esta receta crea un ráster a partir de una matriz.

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
fuente
Esta receta es buena pero hay un problema con el archivo tiff final. Los valores lat-lon de píxeles no son adecuados.
Shubham_geo
Es posible que vea incompatibilidades extrañas entre ESRI WKT y OGC WKT: gis.stackexchange.com/questions/129764/…
Adam Erickson
Una cosa que encontré es que la forma mencionada por usted definitivamente cambiará la matriz a un ráster con facilidad. Pero necesitamos georreferenciar este ráster con la parte superior izquierda e inferior derecha coordinadas usando gdal_translate. Una forma de hacerlo es siguiendo los dos pasos: 1) Primero encuentre los valores lat-lon superiores izquierdo e inferior derecho a través de gdalinfo 2) Luego, a través de gdal_translate, utilice el geotiff (generado con el enfoque mencionado anteriormente de convertir la matriz en ráster) para georreferenciarlo con las coordenadas lat-lon superior izquierda e inferior derecha.
Shubham_geo
0

Una alternativa al enfoque sugerido en las otras respuestas es usar el rasteriopaquete. Tuve problemas para generarlos usando gdaly encontré que este sitio era útil.

Suponiendo que tiene otro archivo tif ( other_file.tif) y una matriz numpy ( numpy_array) que tiene la misma resolución y extensión que este archivo, este es el enfoque que funcionó para mí:

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
fuente