¿Cargar completamente el ráster en una matriz numpy?

26

He estado tratando de verificar mis filtros en el ráster DEM para el reconocimiento de patrones y siempre resulta en la falta de las últimas filas y columnas (como ... 20) . He intentado con la biblioteca PIL, carga de imágenes. Luego con numpy. La salida es la misma.

Pensé, algo estaba mal con mis bucles, al verificar los valores en la matriz (simplemente seleccionando píxeles con identificación en ArcCatalog) me di cuenta de que los valores de píxeles no se cargaban en una matriz.

Entonces, simplemente abriendo, colocando en la matriz y guardando la imagen de la matriz:

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 en cortar las últimas filas y columnas. Lo sentimos, no puedo publicar la imagen

¿Alguien podría ayudar a entender por qué? ¿Y aconsejar alguna solución?

EDITAR:

Entonces, logré cargar pequeños rásteres en una matriz numpy con la ayuda de muchachos, pero cuando tengo una imagen más grande empiezo a obtener errores. Supongo que se trata de los límites de la matriz numpy, por lo que la matriz se reforma automáticamente o algo así ... Entonces, por ejemplo:

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

El punto es que no quiero leer bloque por bloque, ya que necesito filtrar, varias veces con diferentes filtros, diferentes tamaños ... ¿Hay alguna solución o debo aprender rading por bloques: O

najuste
fuente

Respuestas:

42

si tiene enlaces python-gdal:

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

Y tu estas listo:

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
fuente
Sí, con gdal, supongo que no tuve problemas, pero estoy tratando de usar como menos bibliotecas ... Y numpy parecía tan popular para eso 'mientras buscaba en Google'. ¿Alguna idea, de hecho, por qué numpy / PIL deja de cargarse?
Najuste
No lo sé. PIL debe ser lo suficientemente robusto para que se envíe con python. Pero imho geotiff son más que imágenes, por ejemplo, contienen muchos metadatos, y PIL no es (de nuevo en mi humilde opinión) la herramienta adecuada.
nickves
Solo que a veces odio esos requisitos de comillas y barras al abrir los datos ... Pero, ¿qué hay de escribir una matriz numpy en Raster? Se trabaja con la biblioteca PIL, pero utilizando outputRaster.GetRasterBand (1) .WriteArray (myarray) produce trama no válida ..
najuste
no olvide vaciar los datos al disco, con outBand.FlushCache (). Puede encontrar algunos tutoriales aquí: gis.usu.edu/~chrisg/python/2009
nickves
1
Marque " lists.osgeo.org/pipermail/gdal-dev/2010-January/023309.html ", parece que se ha quedado sin memoria RAM.
nickves
21

Puede usar rasterio para interactuar con matrices NumPy. Para leer un ráster en una 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]

Esto leerá todo en una matriz numpy 3D arr, con dimensiones [band, row, col].


Aquí hay un ejemplo avanzado para leer, editar un píxel y luego guardarlo de nuevo en el ráster:

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)

El ráster se escribirá y se cerrará al final de la declaración "con" .

Mike T
fuente
¿Por qué no podemos ver todos los valores cuando escribo print (arr)? Separa los valores con esto ..., ...,?
Mustafa Uçar
@ MustafaUçar así es como NumPy imprime las matrices, que puede modificar . O corte una ventana de la matriz para imprimir, entre muchos otros trucos de Numpy.
Mike T
Una pregunta general. Si quiero generar una única matriz con varias escenas, que son cuatro dimensiones como (escena, altura, ancho, bandas), ¿cómo debo modificar este fragmento?
Ricardo Barros Lourenço
@ RicardoBarrosLourenço Supongo que su cuarta dimensión (¿escena?) Se almacena en cada archivo. Primero llenaría una matriz vacía de 4D vacía, luego recorrería cada archivo (escena) e insertaría la porción 3D de cada uno. Es posible que deba arr.transpose((1, 2, 0))obtener (altura, ancho, bandas) de cada archivo.
Mike T
@MikeT esta población sería como np.append()?
Ricardo Barros Lourenço
3

De acuerdo, estoy leyendo una imagen png antigua y simple, pero esto funciona usando scipy ( imsaveaunque usa 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)

Mi png resultante también es 81 x 90 píxeles.

Chad Cooper
fuente
Gracias, pero estoy tratando de usar menos bibliotecas. Y por ahora puedo hacerlo con gdal + numpy ... (con suerte sin PIL).
Najuste
1
@najuste ¿Qué sistema operativo está activado? Mac y la mayoría de los sabores de Linux vienen con scipyy numpy.
Chad Cooper
Aparentemente ... estoy en Windows, varias versiones de Win. : /
najuste
2

Mi solución usando gdal se ve así. Creo que es muy reutilizable.

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
fuente
0

Estoy usando una imagen hiperespectral con 158 bandas. Quiero calcular la trama. pero consigo

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 el print(data1)obtuve solo un "1", pero los valores reales son algunos flotantes

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 píxel 0,139200

Por favor ayuda a encontrar el error

Kais Tounsi
fuente