[gdal-dev] Creating Geotiff with WriteArray in Python

[Matt Fearon]

Thanks Chris, Steve, and Ivan. I have included comments from both of
you. Chris, thank you for the great detail of your response. However,
I still receive trouble with WriteArray and few additional things now.
Below are my errors and updated code. Any thoughts would be
appreciated.

thanks,
Matt


Traceback (most recent call last):
  File "./test.py", line 38, in <module>
    dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)
  File "/usr/lib/python2.5/site-packages/gdal.py", line 876, in WriteArray
    return gdalnumeric.BandWriteArray( self, array, xoff, yoff )
  File "/usr/lib/python2.5/site-packages/gdalnumeric.py", line 193, in
BandWriteArray
    datatype = NumericTypeCodeToGDALTypeCode( array.typecode() )
AttributeError: 'numpy.ndarray' object has no attribute 'typecode'



#! /usr/bin/env python

import pickle, array, sys, os, shutil
import gdal, gdalnumeric
from gdalconst import *
import numpy
from numpy import *


file = "china/hangzhou/processed/ALPSRP021230580/IMG-HH-ALPSRP021230580-H1.5GUA"

gd=gdal.Open(file)
array=gd.ReadAsArray()
dims = array.shape
nx=dims[1]
ny=dims[0]

nxarray=numpy.zeros(nx, dtype=float32)
cal_factor=-83.0

dst_format = 'GTiff'
dst_datatype = gdal.GDT_Float32
dst_options = ['COMPRESS=LZW']
dst_file='foo.tif'
dst_xsize = nx
dst_ysize = ny
dst_nbands = 1

driver = gdal.GetDriverByName( dst_format )
dst_ds = driver.Create(dst_file, dst_xsize, dst_ysize, dst_nbands,
dst_datatype, dst_options)


for j in range(ny):
    nxarray=array[j,:]
    nxarray = numpy.where(nxarray > 0, 10.0 * log10(pow(nxarray,2)) +
cal_factor, nxarray)
    nxarray.shape = (1,-1)
    print nxarray
    dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)




On Thu, Mar 20, 2008 at 12:52 PM, Christopher Barker
<Chris.Barker at noaa.gov> wrote:
> 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
>
>
> _______________________________________________
>  gdal-dev mailing list
>  gdal-dev at lists.osgeo.org
>  http://lists.osgeo.org/mailman/listinfo/gdal-dev
>