I have been trying to check my filters on DEM raster for pattern recognition and it is always resulting in missing last rows and columns(like..20). I have tried with PIL library, image load. Then with numpy. The output is the same.

I thought, something is wrong with my loops, when checking values in array (just picking pixels with Identification in ArcCatalog) I realized that pixel values were not loaded into an array.

So, just simply opening, puting into array and saving the image from array:

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

Results in cuting away the last rows and columns. Sorry, can#t post the image

Anyone could help to understand why? And advise some solution?


So, I succeeded loading small rasters into numpy array with a help of guys, but when having a bigger image I start getting errors. I suppose it's about the limits of numpy array, and so array is automatically reshaped or smth like that... So ex:

Traceback (most recent call last):
  File "<pyshell#36>", line 1, in <module>
  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

The point is I don't want to read block by block as I need filtering, several times with different filters, different sizes.. Is there any work around or I must learn rading by blocks :O

4 Answers 4


if you have python-gdal bindings:

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

And you're done:

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)
  • Yeah, with gdal I guess I did not have problem, but I'm trying to use as less libraries... And numpy seemed so popular for that 'while googling'. Any idea, indeed, why numpy/PIL stops loading???
    – najuste
    Commented Sep 10, 2012 at 7:51
  • I don't know. PIL should robust enough so its shipped with python. But imho geotiff are more than images - they carry lots of metadata for example- and PIL is not (again imho) the right tool.
    – nickves
    Commented Sep 10, 2012 at 11:37
  • I just sometimes hate those diff quotation and slash requirements, when opening data.. But what about writing numpy array back to Raster? It works with PIL library, but using outputRaster.GetRasterBand(1).WriteArray(myarray) produces invalid raster..
    – najuste
    Commented Sep 10, 2012 at 21:44
  • don't forget to flush the data to the disk, with outBand.FlushCache() . You can find some tutorials here: gis.usu.edu/~chrisg/python/2009
    – nickves
    Commented Sep 10, 2012 at 21:53
  • 1
    Check " lists.osgeo.org/pipermail/gdal-dev/2010-January/023309.html " - it seems you've ran out or ram.
    – nickves
    Commented Sep 23, 2012 at 16:40

You can use rasterio to interface with NumPy arrays. To read a raster to an array:

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]

This will read everything into a 3D numpy array arr, with dimensions [band, row, col].

Here is an advanced example to read, edit a pixel, then save it back to the raster:

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

The raster will be written and closed at the end of the "with" statement.

  • Why can't we see all values when I write print(arr). It seperates values with this ..., ..., ? Commented Nov 22, 2017 at 14:11
  • @MustafaUçar this is how NumPy prints arrays, which you can modify. Or slice a window of the array to print, among many other Numpy tricks.
    – Mike T
    Commented Nov 22, 2017 at 22:29
  • A general question. If I want to output a single array with multiple scenes, being four dimensions such as (scene, height, width, bands), how should I modify this snippet? Commented Dec 19, 2017 at 20:15
  • 1
    @RicardoBarrosLourenço I'm guessing your fourth dimension (scene?) is stored in each file. I'd first populate an empty 4D numpy array, then loop through each file (scene) and insert the 3D portion of each. You may need to arr.transpose((1, 2, 0)) to get (height, width, bands) from each file.
    – Mike T
    Commented Dec 19, 2017 at 22:33
  • @MikeT this population would be such as np.append()? Commented Dec 19, 2017 at 22:37

My solution using gdal looks like this. I think it is very reusable.

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

Granted I'm reading a plain old png image, but this works using scipy (imsave uses PIL though):

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

My resultant png is also 81 x 90 pixels.

  • Thanks, but I'm trying to use as less libraries.. And for now I can make it with gdal+numpy...(hopefully without PIL).
    – najuste
    Commented Sep 10, 2012 at 21:47
  • 1
    @najuste What OS are on? Mac and most Linux flavors come with scipy and numpy. Commented Sep 10, 2012 at 22:18
  • Apparently... I'm on Windows, various versions of Win. :/
    – najuste
    Commented Sep 11, 2012 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.