Gravando numpy array em arquivo raster

30

Eu sou novo no GIS.

Eu tenho algum código que converte imagens infravermelhas de Marte em mapas de inércia térmica, que são armazenados como matrizes numpy 2D. Eu tenho salvado esses mapas como arquivos hdf5, mas eu realmente gostaria de salvá-los como imagens rasterizadas, para poder processá-los no QGIS. Eu passei por várias pesquisas para descobrir como fazer isso, mas sem sorte. Eu tentei seguir as instruções no tutorial em http://www.gis.usu.edu/~chrisg/python/, mas os arquivos que eu produzo usando o código de exemplo são abertos como caixas cinza simples quando os importo para o QGIS. Sinto que se alguém pudesse sugerir o procedimento mais simples possível para um exemplo simplificado do que eu gostaria de fazer, talvez eu pudesse fazer algum progresso. Eu tenho QGIS e GDAL, ficaria muito feliz em instalar outras estruturas que qualquer pessoa poderia recomendar. Eu uso o Mac OS 10.7.

Portanto, se, por exemplo, eu tenho uma matriz numpy de inércia térmica que se parece com:

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

E para cada pixel eu tenho a latitude e 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) ) 

Qual procedimento as pessoas recomendariam para converter esses dados em um arquivo raster que eu possa abrir no QGIS?

EddyTheB
fonte
A que slide do tutorial você está se referindo?
RK

Respostas:

23

Uma solução possível para o seu problema: converta-o em uma varredura ASCII, cuja documentação está aqui . Isso deve ser bastante fácil de fazer com python.

Portanto, com seus dados de exemplo acima, você terminaria com o seguinte em um arquivo .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

Isso adiciona com sucesso ao QGIS e ao ArcGIS e, estilizado no ArcGIS, fica assim: versão raster acima

Adendo: Embora você possa adicioná-lo ao QGIS conforme observado, se você tentar acessar as propriedades dele (para estilizá-lo), o QGIS 1.8.0 trava. Estou prestes a relatar isso como um bug. Se isso acontecer com você também, existem muitos outros SIG gratuitos por aí.

GIS-Jonathan
fonte
Isso é fantástico, obrigado. E imagino que, tendo escrito minha matriz como um arquivo ascii, eu possa convertê-la em um formato binário usando uma função de conversão pré-escrita.
EddyTheB 22/10/12
Para sua informação, não tive o problema pendente do QGIS, também tenho a versão 1.8.0.
EddyTheB 22/10/12
31

Abaixo está um exemplo que escrevi para um workshop que utiliza os módulos numpy e gdal Python. Ele lê os dados de um arquivo .tif em uma matriz numpy, faz uma reclassificação dos valores na matriz e os grava de volta em um .tif.

Pela sua explicação, parece que você conseguiu escrever um arquivo válido, mas você só precisa simbolizá-lo no QGIS. Se bem me lembro, quando você adiciona um raster pela primeira vez, ele geralmente mostra todas as cores se você não possui um mapa de cores 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
fonte
1
+1 por rubor - estava batendo minha cabeça contra a parede, tentando descobrir como 'salvar' a coisa!
Badgley
Eu tive que adicionar outDs = Nonepara obtê-lo salvo
JaakL
23

Finalmente encontrei esta solução, obtida com essa discussão ( http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html ). Gosto porque posso ir direto de uma matriz numpy para um arquivo raster tif. Ficaria muito grato por comentários que poderiam melhorar a solução. Vou postá-lo aqui, caso alguém procure uma resposta semelhante.

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
fonte
3
A "rotação ocorre duas vezes" para explicar o efeito de um bit girado de y em x e o bit girado de x em y. Veja lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html, que tenta explicar as inter-relações entre os parâmetros de "rotação".
Dave X
Este post é realmente útil, obrigado. No meu caso, no entanto, estou recebendo um arquivo tif completamente preto quando o abro como uma imagem fora do ArcGIS. Minha referência espacial é a British National Grid (EPSG = 27700), e as unidades são metros.
FaCoffee
Eu postei uma pergunta aqui: gis.stackexchange.com/questions/232301/…
FaCoffee
Você descobriu como definir a codificação IAU2000: 49900 Mars?
K.-Michael Aye
4

Há também uma boa solução no livro de receitas oficial do GDAL / OGR para Python.

Esta receita cria uma varredura a partir de uma 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
fonte
Esta receita é boa, mas há um problema com o arquivo tiff final. Os valores lat-lon de pixels não são adequados.
Shubham_geo
Você pode estar vendo incompatibilidades estranhas entre o ESRI WKT e o OGC WKT: gis.stackexchange.com/questions/129764/…
Adam Erickson
Uma coisa que encontrei é que a maneira mencionada por você definitivamente mudará a matriz para uma varredura com facilidade. Mas precisamos georreferenciar essa varredura com a parte superior esquerda e a parte inferior direita coordenadas usando gdal_translate. Uma maneira de fazer isso é seguir os dois passos: 1) Primeiro encontre os valores de lat-lon superior esquerdo e inferior direito através do gdalinfo. para georreferenciá-lo com as coordenadas lat-lon superior esquerda e inferior direita.
Shubham_geo 18/08
0

Uma alternativa à abordagem sugerida nas outras respostas é usar o rasteriopacote. Ocorreu um problema ao gerá-los gdale achei este site útil.

Supondo que você tenha outro arquivo tif ( other_file.tif) e uma matriz numpy ( numpy_array) que tenha a mesma resolução e extensão que este arquivo, esta é a abordagem que funcionou para mim:

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
fonte