[gdal-dev] Creating Geotiff with WriteArray in Python

[Christopher Barker]

Matt Fearon wrote:
> The problem is related to the line
> "dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)" - I am trying to
> write an imported integer array(y,x) , 13700x13300, to a tif file row
> by row 

> BandWriteArray
>     xsize = array.shape[1]
> IndexError: tuple index out of range

It looks like RasterBand.WriteArray is expecting a 2-d array, and you 
are passing in a 1-d array. This makes sense as you are trying to pass 
in a single row. Do  you need to do it row by row due to memory limitations?

> import Numeric

> import numpy

which version of GDAL are you using? if you use 1.5, you can just use 
numpy, rather than both Numeric and numpy.

> from math import *

Just so you know both Numeric an numpy are essentially supersets of 
math, so you usually don't need the math module. Also, as math and 
Numeric have a lot of overlapping symbols, you generally don't want to 
"import *" both of them

> nxarray=zeros(nx,Float)

> dst_datatype = gdal.GDT_Float32

Another issue here -- I think Numeric.Float is a 64bit float (C double), 
same as the built-in python float. If you want a Float32 array, you need 
something else. In numpy, it's "float32", in Numeric, it's "Float32".

> for j in range(0,ny-1):
>   for i in range(0,nx-1):

I think you want "range(ny)" here:

 >>> range(4)
[0, 1, 2, 3]

There is also room to "vectorize" this operation -- you should be able 
to process a whole row at once, and only loop through y, but not x -- 
that will speed things up a lot (see below). If you have enough memory, 
you could process the whole thing at once with numpy operations -- 
that's the whole point of numpy(and Numeric).

>     nxarray[i]=0.0
>     tmp = array[j,i] * 1.0

why multiply by 1.0 here? that shouldn't do anything.

>     if (tmp > 0):
>       nxarray[i] = 10.0 * log10(pow(tmp,2)) + cal_factor
> 
>   dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)

Here's the issue. If I have this right, WriteArray is expecting a 2-d 
array, and writing to the RasterBand at the given offsets. So what you 
need to do is convert your (nx,) array to a (nx,1) array, and it should 
work:

nxarray.shape = (-1,1)

This tells Numeric or numpy to set the second dimension to 1, and the 
first to whatever it needs to be to fit. (that could be (nx, 1), but I 
like to be generic when I can.

Now you can call WriteRaster:

dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)

Also:

dst_ds = None

I think that will delete your destination raster -- has it been written 
to file yet?

Here is my (untested) suggestion to eliminating a loop:

for j in range(ny): # loop though the rows
     nxarray=array[j,:] # note - this is a view into your original array
     nxarray = numpy.where(nxarray > 0,
                           10.0 * log10(pow(tmp,2)) + cal_factor,
                           nxarray)
    nxarray.shape = (-1, 1)
    dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)


Hope that helps.

-Chris

-- 
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov