[gdal-dev] empty Gtiff from Python

[Vincent Schut]

Tom van der Putte wrote:
> Hi List,
>
> I'm trying to save a numpy array to a GTiff, which should be quite 
> straightforward, considering all the examples available on the web. 
> However, I get an file that is not readable by Qgis, is regarded by 
> ArcMap to have only 0 and gdalinfo displays nothing on the values.
>
> I start with a ndarray, dtype = int32, shape = ( 10,10). It has 
> integer values between 1 and 10. The code is as follows:
>
> import numpy as np
> from osgeo import gdal
> from osgeo import osr
>
> #np array
> a = np.ones((10,10))
>
> #coord-system
> sr_str = 'LOCAL_CS["arbitrary"]'
> sr = osr.SpatialReference( sr_str )
>
> #driver
> driver = gdal.GetDriverByName("GTiff")
>
> #dataset
> ds = driver.Create("D:\\temp\\output.tif", 10, 10, 1, gdal.GDT_Byte)
> ds.SetProjection( sr_str )
> ds.SetGeoTransform((0, 1, 0,0, 0,1))
>
> #write and close
> ds.GetRasterBand(1).WriteArray(a)
> ds = None

Tom,

don't know if you can use WriteArray this way... what I always do, is:

from osgeo import gdal_array
# to write numpy array:
gdal_array.BandWriteAsArray(ds.GetRasterBand(1), a)

Best,
Vincent.
>
> The corresponding gdalinfo is:
> *********
> Raster dataset parameters:
>   Projection: 
> LOCAL_CS["arbitrary",UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
>   RasterCount: 1
>   RasterSize (10,10)
> Using driver GeoTIFF
>   Metadata:
>     0:  AREA_OR_POINT=Area
>
>   Image Structure Metadata:
>     0:  INTERLEAVE=BAND
>
> Corner Coordinates:
>   Upper Left (0, 0)
>   Lower Left (0, 10)
>   Upper Right (10, 0)
>   Lower Right (10, 10)
>   Center (5, 5)
>
> Coordinate System is:
> LOCAL_CS["arbitrary",
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]]]
> Band 1 :
>    DataType: Byte
>    ColorInterpretation: Gray
>    Description:
>    Size (10,10)
>    BlockSize (10,10)
> *********
>
> I have tried both version 1.6.3 (I'm trying to make a Qgis plugin, and 
> thus restricted to 1.6.3) as well as 1.8.1. I have also tried ds = 
> driver.Create("D:\\temp\\output.tif", 10, 10, 1, gdal.GDT_Int32), but 
> that didn't make much difference.
>
> I tried to make a Surfer grid (GSBG), and this somewhat worked: it 
> gave a nice picture, but the values were off the charts.
>
> Any clues as to what I'm missing here? Cheers,
>
> Tom van der Putte
>
>
>
> _______________________________________________
> 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/20111028/41a3422a/attachment.html