[gdal-dev] Correct way to create a new GTiff, with numpy.zeros() and python bindings (resolved)

[Dane Springmeyer]

Hi Antonio,

Thanks for the suggestion to check my numpy types. I dug into the  
various Gdal types and numpty types and saw that I'll have to be  
careful to keep sync.

However, with a little more digging I think I found the simple  
problem: my axis needed to be reversed when calling the numpy function.

This code below now works:


# create a 4 band geotiff with RGB bands all zero and Alpha band fully  
non-transparent

from osgeo import gdal
from osgeo import osr
import numpy

driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create( "test.tif", 7850, 3500, 4, gdal.GDT_Byte)
dst_ds.SetGeoTransform( [-124.75, .001, 0.0, 49.02, 0.0, -.001] )
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
dst_ds.SetProjection( srs.ExportToWkt() )
zeros = numpy.zeros( (3500, 7850), numpy.uint8 )
dst_ds.GetRasterBand(1).WriteArray( zeros )
dst_ds.GetRasterBand(2).WriteArray( zeros )
dst_ds.GetRasterBand(3).WriteArray( zeros )
opaque = numpy.ones((3500,7850), numpy.uint8 )*255
dst_ds.GetRasterBand(4).WriteArray( opaque )


Cheers,

Dane


On Jun 7, 2008, at 1:59 AM, Antonio Valentino wrote:

> Il giorno Fri, 6 Jun 2008 17:20:44 -0700
> Dane Springmeyer <blake at hailmail.net> ha scritto:
>
>> Hi list,
>>
>> Im just starting to use the python bindings and have been
>> experimenting with gdal_rasterize which great results, effectively
>> burning in multiple different colors based on some related vector
>> shapefile attributes.  The trick to getting started was that I did
>> not have any raster (at the right extent and resolution) to copy
>> from, so to get started I followed the api tutorial to create and new
>> Gtiff.
>>
>> This code works perfectly for my usecase:
>>
>> from osgeo import gdal
>> from osgeo import osr
>> import numpy
>>
>> driver = gdal.GetDriverByName('GTiff')
>> dst_ds = driver.Create( "/Users/spring/projects/wind/data/processed/
>> wpc_final.tif", 7850, 3500, 4, gdal.GDT_Byte)
>> dst_ds.SetGeoTransform( [-124.75, .001, 0.0, 49.02, 0.0, -.001] )
>> srs = osr.SpatialReference()
>> srs.ImportFromEPSG(4326)
>> dst_ds.SetProjection( srs.ExportToWkt() )
>>
>>
>> But the tutorial also has the step of writing zeros over the image
>> with numpy. Since I am already burning in the required values into
>> my raster with gdal_rasterize this seems not to matter, but I am
>> really curious for future uses how to get this step working. Starting
>> with zeros, or any default value would be useful.
>>
>> So, if I execute the following code block next:
>>
>> zeros = numpy.zeros( (7850, 3500) )
>
> Dane,
> here you are getting a floating point array. Try to print zeros.dtype.
> In my opinion the correct code should be
>
> zeros = numpy.zeros( (7850, 3500), dtype=numpy.uint8 )
>
> or something similar.
>
>> dst_ds.GetRasterBand(1).WriteArray( zeros )
>> dst_ds.GetRasterBand(2).WriteArray( zeros )
>> dst_ds.GetRasterBand(3).WriteArray( zeros )
>> dst_ds.GetRasterBand(4).WriteArray( zeros )
>>
>>
>> I get the error:
>> 'ValueError: array larger than output file, or offset off edge'
>> which obviously means that I am writing the wrong size numpy array
>> into my raster.
>>
>> So, what do I need to know to figure out the right dimensions for
>> the numpy.zeros() function? Some slick fraction of my image size in
>> pixels and my geotransform?
>>
>> Thanks,
>>
>> Dane
>>
>
>
> -- 
> Antonio Valentino
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080607/012b46f4/attachment.html