[Gdal-dev] How to efficiently read in a GDAL image into a matrix

[Ryan Lewis]

Hi,

I would like to read in an entire ENVI image into a GDAL Dataset, and
then, create a matrix the size of which is pixels X number of raster
bands
So each column of the matrix is a raster band of the image, and each
row is a pixel of the image. so the idea is I would linearize the 2d
raster band with something like the pixel at (x,y) in a particular
rasterband would be at the position
x*width + y in the column.

Apparently it is a really bad idea to try and read one value at a
time out of each raster band, so could someone just tell me the best
method for doing this.

(If you have the code already that would be cool, but it is not
necessary, just an outline of the best way to do this in GDAL would
be awesome).

Thanks,

- Ryan H. Lewis

[Nate Jennings]

Ryan,

You might try my website below, then go to Skillset--Application
Development. I just put up a python script to do a mean filter on a
multiband image. It uses the python gdal build of 1.4.2 that is posted on
the gdal site and also requires the old Numeric libraries as well. I based
my code on sample code from the gdal site as well as my past experience in
accessing pixels within bands from a graduate remote sensing algorithms
course I took during grad school. I don't know if this is what you need,
but I do access individual pixels. The code runs very slow (~5 min for a
400 x 600 image). I am looking to reproduce a similar function using c++ if
and when I can understand how the gdal libraries are used within c++.

Nate

[Christopher Barker]

Post by Nate Jennings
You might try my website below, then go to Skillset--Application
Development. I just put up a python script to do a mean filter on a
multiband image. It uses the python gdal build of 1.4.2 that is
posted on the gdal site and also requires the old Numeric libraries
as well.

You should be able to do it with the latest GDAL and numpy as well.

Post by Nate Jennings
I do access individual pixels. The code runs very slow (~5 min for a
400 x 600 image).

Is it necessary to do by individual pixels? numpy gives you a lot of
tools to work with whole arrays at a time, which would make it MUCH
faster. There are also a fair number of image processing algorithms out
there that use numpy -- ask on the numpy list if you can't find what you
need.

Post by Nate Jennings
I am looking to reproduce a similar function using c++ if and when I
can understand how the gdal libraries are used within c++.

If you do really need to go to C/C++, another option is to write code
that works with numpy arrays as a python extension -- that way you can
still write everything else in Python. There are lots of way to write
C/C++ for use with numpy arrays:

SWIG
Pyrex
weave.blitz
f2py

See: http://www.scipy.org/PerformancePython

and the scipy and numpy lists for examples and discussion.

-Chris
--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R [phone] voice
7600 Sand Point Way NE [phone] fax
Seattle, WA 98115 [phone] main reception

***@noaa.gov

[Frank Warmerdam]

Post by Ryan Lewis
Hi,
I would like to read in an entire ENVI image into a GDAL Dataset, and
then, create a matrix the size of which is pixels X number of raster bands
So each column of the matrix is a raster band of the image, and each row
is a pixel of the image. so the idea is I would linearize the 2d raster
band with something like the pixel at (x,y) in a particular rasterband
would be at the position
x*width + y in the column.

Ryan,

I don't understand how you would use x*width+y to get to a particular pixel.
Perhaps you mean y*width + x?

If you have an image that is 200 width, 300 height and 3 bands you could
read it into a single "band sequential" array like this:

GByte data[200*300*3];

...

GDALDatasetRasterIO( hDS, GF_Read, 0, 0, 200, 300,
data, 200, 300, GDT_Byte,
3, NULL, 1, 200, 200*300 );

If you want to have an array of 2D arrays, one for each band, you might do:

GByte *bands[3];
GByte band1[200*300];
GByte band2[200*300];
GByte band3[200*300];
GByte *bands[3] = {band1,band2,band3};

for( int band = 1; band <= 3; band++ )
GDALRasterIO( GDALGetRasterBand(hDS,band), GF_Read,
0, 0, 200, 300,
bands[band-1], 200, 300, GDT_Byte, 0, 0 );

Key things to keep in mind.

o If you want to read more than one band at a time, you need to use
GDALDatasetRasterIO() (or the RasterIO() method on GDALDataset).
o The final 2/3 values in the RasterIO() call control the packing
of the buffer. Values of 0,0,0 are default packing, otherwise it
is "pixel step", "line step" and "band step".

I'm open to contributed improvements to the GDAL raster api tutorial,
and the appropriate method docs to make this stuff clearer.

Best regards,

[Jose Luis Gomez Dans]

Hi Ryan,

Post by Ryan Lewis
I would like to read in an entire ENVI image into a GDAL Dataset, and
then, create a matrix the size of which is pixels X number of raster
bands.

I do this sort of thing in Python (+numpy et al.) fairly often. Unfortunately, I haven't got gdal, python or anything installed over here, but I think the code goes like this:

g=gdal.Open ( fname )
for band in xrange g.GetRasterCount):
b = g.GetRasterBand( band ).ReadAsArray()
#b is a 2D array with your band.
#Use numpy's reshape to stack it however you like
#In particular, you can use slice notation for this sort of thing

Using large-ish images (Landsat, around 6400x6500 pixels and a few bands), it is fairly confortable to read the data and apply algorithms using numpy and scipy. If you apply this to whole images, it is very fast. If you're using a few pixels at a time, use a masked array to improve performance.

Finally, another way to speed things up is to read N rows from B bands at a time and work on smaller chunks. If you have dual cores or several machines, you can try to process "several smaller chunks" in parallel, but I haven't tried that yet.

Post by Ryan Lewis
Apparently it is a really bad idea to try and read one value at a
time out of each raster band, so could someone just tell me the best
method for doing this.

You might as well read a line (or several lines into a buffer, and deal with them. It is significantly faster (you make use of IO cache, processor cache, etc.) than going pixel by pixel.

I am sorry I can't provide any code (I have done the approaches mentioned here in both python and C. In the end, when factoring the debugging and me screwing up with things, python ends up being 100x faster than C :D).

Hope that helps,
Jose