Thursday 25 August 2016

python - Fully load raster into a numpy array?


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)
Image.fromarray(a).save(outraster)

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?



EDIT:


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 "", 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: 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



Answer



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:


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)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...