Carregar totalmente a varredura em uma matriz numpy?

26

Eu tenho tentado verificar meus filtros no DEM raster para reconhecimento de padrões e isso sempre resulta em falta das últimas linhas e colunas (como 20) . Eu tentei com a biblioteca PIL, carregamento de imagem. Então com entorpecido. A saída é a mesma.

Eu pensei que algo estava errado com meus loops, ao verificar os valores na matriz (apenas escolhendo pixels com Identificação no ArcCatalog), percebi que os valores de pixel não eram carregados em uma matriz.

Então, basta abrir, colocar no array e salvar a imagem do array:

a=numpy.array(Image.open(inraster)) #raster is .tif Float32, size 561x253
newIm=Image.new(Im.mode, Im.size)
Image.fromarray(a).save(outraster)

Resulta no corte das últimas linhas e colunas. Desculpe, não posso postar a imagem

Alguém poderia ajudar a entender o porquê? E aconselhar alguma solução?

EDITAR:

Então, consegui carregar pequenos rasters em um array numpy com a ajuda de caras, mas ao ter uma imagem maior, começo a receber erros. Suponho que seja sobre os limites da matriz numpy e, portanto, a matriz é automaticamente remodelada ou algo assim ... Então ex:

Traceback (most recent call last):
  File "<pyshell#36>", line 1, in <module>
    ima=numpy.array(inDs.GetRasterBand(1).ReadAsArray())
  File "C:\Python25\lib\site-packages\osgeo\gdal.py", line 835, in ReadAsArray
    buf_xsize, buf_ysize, buf_obj )
  File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 140, in BandReadAsArray
    ar = numpy.reshape(ar, [buf_ysize,buf_xsize])
  File "C:\Python25\lib\site-packages\numpy\core\fromnumeric.py", line 108, in reshape
    return reshape(newshape, order=order)
ValueError: total size of new array must be unchanged

O ponto é que eu não quero ler bloco por bloco, pois preciso filtrar, várias vezes com filtros diferentes, tamanhos diferentes. Existe alguma solução alternativa ou devo aprender a radar por blocos: O

najuste
fonte

Respostas:

42

se você tiver ligações python-gdal:

import numpy as np
from osgeo import gdal
ds = gdal.Open("mypic.tif")
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())

E pronto:

myarray.shape
(2610,4583)
myarray.size
11961630
myarray
array([[        nan,         nan,         nan, ...,  0.38068664,
     0.37952521,  0.14506227],
   [        nan,         nan,         nan, ...,  0.39791253,
            nan,         nan],
   [        nan,         nan,         nan, ...,         nan,
            nan,         nan],
   ..., 
   [ 0.33243281,  0.33221543,  0.33273876, ...,         nan,
            nan,         nan],
   [ 0.33308044,  0.3337177 ,  0.33416209, ...,         nan,
            nan,         nan],
   [ 0.09213851,  0.09242494,  0.09267616, ...,         nan,
            nan,         nan]], dtype=float32)
nickves
fonte
Sim, com o gdal, acho que não tive problemas, mas estou tentando usar como menos bibliotecas ... E o numpy parecia tão popular por isso 'enquanto pesquisava'. Alguma idéia, de fato, por que numpy / PIL para de carregar ???
najuste
Eu não sei. O PIL deve ser robusto o suficiente para ser enviado com python. Mas imho geotiff são mais que imagens - eles carregam muitos metadados, por exemplo - e o PIL não é (novamente imho) a ferramenta certa.
nickves 10/09/12
Às vezes, odeio esses requisitos de cotação e barra, ao abrir dados. Mas e quanto a escrever um array numpy de volta ao Raster? Ele funciona com a biblioteca PIL, mas usando outputRaster.GetRasterBand (1) .WriteArray (myarray) produz raster inválido ..
najuste
não esqueça de liberar os dados para o disco, com outBand.FlushCache (). Você pode encontrar alguns tutoriais aqui: gis.usu.edu/~chrisg/python/2009
nickves
11
Marque " lists.osgeo.org/pipermail/gdal-dev/2010-January/023309.html " - parece que você acabou ou se foi.
nickves 23/09/12
21

Você pode usar o rasterio para fazer interface com matrizes NumPy. Para ler uma varredura em uma matriz:

import rasterio

with rasterio.open('/path/to/raster.tif', 'r') as ds:
    arr = ds.read()  # read all raster values

print(arr.shape)  # this is a 3D numpy array, with dimensions [band, row, col]

Isso lerá tudo em uma matriz numpy 3D arr, com dimensões [band, row, col].


Aqui está um exemplo avançado para ler, editar um pixel e salve-o novamente na varredura:

with rasterio.open('/path/to/raster.tif', 'r+') as ds:
    arr = ds.read()  # read all raster values
    arr[0, 10, 20] = 3  # change a pixel value on band 1, row 11, column 21
    ds.write(arr)

A varredura será gravada e fechada no final da instrução "with" .

Mike T
fonte
Por que não podemos ver todos os valores quando escrevo print (arr). Separa valores com isso ..., ...,?
Mustafa Uçar
@ MustafaUçar é assim que o NumPy imprime matrizes, que você pode modificar . Ou corte uma janela da matriz para imprimir, entre muitos outros truques do Numpy.
Mike T
Uma pergunta geral. Se eu quiser gerar uma única matriz com várias cenas, sendo quatro dimensões, como (cena, altura, largura, faixas), como devo modificar esse trecho?
Ricardo Barros Lourenço
@ RicardoBarrosLourenço Acho que sua quarta dimensão (cena?) Está armazenada em cada arquivo. Primeiro eu preenchia uma matriz numpy 4D vazia, depois passava por cada arquivo (cena) e inseria a parte 3D de cada uma. Pode ser necessário arr.transpose((1, 2, 0))obter (altura, largura, faixas) de cada arquivo.
Mike T
@ MikeT essa população seria como np.append()?
Ricardo Barros Lourenço
3

É verdade que estou lendo uma imagem png antiga, mas isso funciona usando o scipy ( imsaveembora use PIL):

>>> import scipy
>>> import numpy
>>> img = scipy.misc.imread("/home/chad/logo.png")
>>> img.shape
(81, 90, 4)
>>> array = numpy.array(img)
>>> len(array)
81
>>> scipy.misc.imsave('/home/chad/logo.png', array)

Meu png resultante também tem 81 x 90 pixels.

Chad Cooper
fonte
Obrigado, mas estou tentando usar como menos bibliotecas .. E por enquanto posso fazê-lo com gdal + numpy ... (espero sem PIL).
najuste
11
@najuste Em que SO estão? Mac e a maioria dos sabores Linux vêm com scipye numpy.
Chad Cooper
Aparentemente ... Estou no Windows, várias versões do Win. : /
najuste 11/09/12
2

Minha solução usando gdal se parece com isso. Eu acho que é muito reutilizável.

import gdal
import osgeo.gdalnumeric as gdn

def img_to_array(input_file, dim_ordering="channels_last", dtype='float32'):
    file  = gdal.Open(input_file)
    bands = [file.GetRasterBand(i) for i in range(1, file.RasterCount + 1)]
    arr = np.array([gdn.BandReadAsArray(band) for band in bands]).astype(dtype)
    if dim_ordering=="channels_last":
        arr = np.transpose(arr, [1, 2, 0])  # Reorders dimensions, so that channels are last
    return arr
MonsterMax
fonte
0

estou usando uma imagem hiperespectral com 158 bandas. Eu quero calcular raster. mas eu entendo

import gdal # Import GDAL library bindings
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
import pylab as plt
import numpy as np
import xlrd
# The file that we shall be using
# Needs to be on current directory
filename = ('C:/Users/KIFF/Desktop/These/data/Hyperion/10th_bandmathref')
outFile = ('C:/Users/KIFF/Desktop/These/data/Hyperion/Math')
XLS=('C:/Users/KIFF/Desktop/These/data/Coef/bcoef.xlsx')
wb = xlrd.open_workbook(XLS)
sheet = wb.sheet_by_index(0)
sheet.cell_value(0, 0)


g = gdal.Open(filename, GA_ReadOnly)

# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
    print ("Problem opening file %s!" % filename)
else:
    print ("File %s opened fine" % filename )

#band_array = g.ReadAsArray()
#print(band_array)
print ("[ RASTER BAND COUNT ]: ", g.RasterCount)

for band in range( g.RasterCount ):
    print (band)
    band += 1
    outFile = ('C:/Users/KIFF/Desktop/These/data/Results/Temp/Math_1_sur_value'+str(band)+'.tiff')
    #print ("[ GETTING BAND ]: ", band )
    srcband = g.GetRasterBand(band)
    if srcband is None:
        continue
    data1 = BandReadAsArray(srcband).astype(np.float)
    print(data1)
   # for i in range(3,sheet.nrows):
    b=sheet.cell_value(band+2,1)
    #print(b)
    dataOut = (1/data1)
    driver = gdal.GetDriverByName("ENVI")
    dsOut = driver.Create(outFile, g.RasterXSize, g.RasterYSize, 1)
    CopyDatasetInfo(g,dsOut)
    bandOut=dsOut.GetRasterBand(1)
    BandWriteArray(bandOut, dataOut)

para o print(data1)eu tenho apenas alguns "1", mas os valores reais são alguns carros alegóricos

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. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
2

Valor de pixel 0,139200

Plz ajuda a encontrar o erro

Kais Tounsi
fonte