[gdal-dev] Creating Geotiff with WriteArray in Python

[gaffigan at sfos.uaf.edu]

Hello, Matt.  I think the array argument to the WriteArray function must
be 2D.  In your example, nxarray should have shape (1,nx), instead of
(nx,).  So, on the relevant line, try the following instead:
  dst_ds.GetRasterBand(1).WriteArray(reshape(nxarray,(1,nx)),0,j)

Hope this gets you closer.

Steve

> Dear GDAL/Python Users,
>
> I am struggling to make the following python/GDAL script work. I am
reading/converting large Palsar images to GeoTiff, and therefore I have
to be concerned with memory. I believe I have the following script very
close to working. 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
using a for loop structure as a converted float. I constructed this code
from found examples, ie, 4bandto3band.py. However,
> "WriteArray" is giving me trouble.  I believe the problem is something
very simple; the arguments in WriteArray are not quite right. I am quite
new to GDAL and Python. Some help would be very much
> appreciated. Below are the errors I receive and my script.
> thanks,
> Matt
>
>
> Traceback (most recent call last):
>   File "./test.py", line 43, 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 187, in
> BandWriteArray
>     xsize = array.shape[1]
> IndexError: tuple index out of range
>
>
>
>
> #! /usr/bin/env python
>
> import Numeric
> from Numeric import *
> import pickle, array, sys, os, shutil
> import gdal, gdalnumeric
> from gdalconst import *
> import numpy
> from math 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=zeros(nx,Float)
> 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(0,ny-1):
>   for i in range(0,nx-1):
>
>     nxarray[i]=0.0
>     tmp = array[j,i] * 1.0
>     if (tmp > 0):
>       nxarray[i] = 10.0 * log10(pow(tmp,2)) + cal_factor
>
>   dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)
>
> dst_ds = None
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>