[gdal-dev] Creating Geotiff with WriteArray in Python

[Christopher Barker]

Matt Fearon wrote:
> Chris, thank you for the great detail of your response.

You're welcome.

> However,
> I still receive trouble with WriteArray and few additional things now.
> Below are my errors

>     datatype = NumericTypeCodeToGDALTypeCode( array.typecode() )
> AttributeError: 'numpy.ndarray' object has no attribute 'typecode'

Sorry, my fault -- I think what's happening is that you're passing in a 
numpy array, and it's expecting a Numeric array -- they are very 
similar, but not identical. It seems you are using GDAL 1.4, in which 
case it's probably better to just stick with Numeric all around.

Note to Frank -- we could wrap input arrays with Numeric.asarray() 
calls, so that the conversion would happen automatically.

> #! /usr/bin/env python
> 
> import pickle, array, sys, os, shutil
> import gdal, gdalnumeric
> from gdalconst import *
> import numpy
> from numpy import *

I like to do:

import numpy as N
or import Numeric as N

In this case, you're probably better off with just Numeric. And I really 
recommend against "import *". In fact, in this case, if you do that, 
you're stomping on the numpy array constructor with your "array" variable.


> gd=gdal.Open(file)
> array=gd.ReadAsArray()

I think this is now a Numeric array, check like this:

type(array) -- if it says:

<type 'array'>
it's Numeric, if it says:

<type 'numpy.ndarray'>
it's numpy.

I think you're working with Numeric here.

By the way -- what a mess! Oh well, everyone is moving to numpy now, so 
in a bit this will no longer be an issue. If you can upgrade to GDAL 1.5 
now, I would. But for now, let's go Numeric:

> nxarray=numpy.zeros(nx, dtype=float32)

You don't need that anymore, it's getting created inside the loop anyway.

> for j in range(ny):
>     nxarray=array[j,:]

I now I gave you this code, but there is an issue here, what type is 
your input data? IIUC, it's not Float32, so you need to do a conversion 
here:
     nxarray = array[j,:].astype(N.Float32)

here's the trick:
     nxarray.savespace(1)
This tells the Numeric array that you don't want it to upcast to a 
Float64 when you multiply it by a standard Python float (which is 64 
bit). One advantage to numpy is that it does this by default.

     nxarray = N.where(nxarray > 0, 10.0 * N.log10(pow(nxarray,2)) +

changed the "numpy" to N -- then if you did "import Numeric as N", 
you'll get a Numeric array here. You can just do "import Numeric", but 
then you have to type that long word everywhere...

> cal_factor, nxarray)
>     nxarray.shape = (1,-1)
>     print nxarray
>     dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)

that should do it.


-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