¿Cómo cargar completamente un ráster en una matriz numpy?

26

He estado intentando verificar mis filtros en el ráster DEM para el reconocimiento de patrones y siempre es el resultado de que faltan 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 está mal con mis bucles, al verificar los valores en la matriz (solo seleccionando píxeles con Identificación en ArcCatalog) me di cuenta de que los valores de píxeles no se cargaron en una matriz.

Entonces, simplemente abriendo, colocando una matriz y guardando la imagen de una 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)

Resultados en recortar las últimas filas y columnas. Lo sentimos, no puedo # t 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 chicos, pero cuando tengo una imagen más grande, empiezo a tener errores. Supongo que se trata de los límites de la matriz numpy, por lo que la matriz se remodela automáticamente o es 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 filtrado, varias veces con diferentes filtros, diferentes tamaños .. ¿Hay algún trabajo alrededor o debo aprender a radiar por bloques: O

    
pregunta najuste 07.09.2012 - 14:46

4 respuestas

38

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 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)
    
respondido por el nickves 08.09.2012 - 23:03
18

Puede utilizar rasterio para interactuar con las matrices NumPy. Para leer un raster a 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 la trama:

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 "con" declaración .

    
respondido por el Mike T 12.01.2015 - 05:08
3

Concedido, estoy leyendo una imagen png antigua, pero esto funciona con scipy ( imsave 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.

    
respondido por el Chad Cooper 08.09.2012 - 00:28
1

Mi solución con 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
    
respondido por el MonsterMax 17.05.2018 - 15:08

Lea otras preguntas en las etiquetas