¿Cómo cargar por completo un raster en un arreglo de numpy?

Question

He estado intentando verificar mis filtros en raster DEM para el reconocimiento de patrones y siempre resulta en la falta de las últimas filas y columnas (alrededor de 20). He probado con la biblioteca PIL, carga de imagen. Luego con numpy. La salida es la misma.

Pensé que algo estaba mal con mis bucles, al verificar los valores en un array (simplemente seleccionando píxeles con Identificación en ArcCatalog) me di cuenta de que los valores de píxeles no se cargaron en un array.

Entonces, simplemente abriendo, poniéndolo en un array y guardando la imagen desde el array:

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

Produce como resultado cortar las últimas filas y columnas. Lo siento, no puedo publicar la imagen

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

EDITAR:

Así que logré cargar pequeños rasters en un array numpy con la ayuda de algunos amigos, pero cuando tengo una imagen más grande comienzo a recibir errores. Supongo que se trata de los límites del array numpy, y el array se redimensiona automáticamente o algo así... Así que por ejemplo:

Traceback (most recent call last):
  File "", line 1, in 
    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: el tamaño total del nuevo array debe permanecer inalterado

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 a leer por bloques :O

Answer 1

Si tienes las conexiones 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)

Answer 2

Puedes usar rasterio para interactuar con arrays de NumPy. Para leer un raster y guardarlo en un array:

import rasterio

with rasterio.open('/path/to/raster.tif', 'r') as ds:
    arr = ds.read()  # lee todos los valores del raster

print(arr.shape)  # este es un array de NumPy en 3D, con dimensiones [banda, fila, columna]

Esto leerá todo en un array de NumPy en 3D arr, con dimensiones [banda, fila, columna].

---

Aquí tienes un ejemplo avanzado para leer, editar un píxel y luego guardarlo de nuevo en el raster:

with rasterio.open('/path/to/raster.tif', 'r+') as ds:
    arr = ds.read()  # lee todos los valores del raster
    arr[0, 10, 20] = 3  # cambia el valor de un píxel en la banda 1, fila 11, columna 21
    ds.write(arr)

El raster se escribirá y se cerrará al final de la declaración "with".

Answer 3

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

Answer 4

Concedido que estoy leyendo una simple imagen png, pero esto funciona usando scipy (imsave utiliza PIL sin embargo):

>>> 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 imagen png resultante también es de 81 x 90 píxeles.