[gdal-dev] empty Gtiff from Python

[Tom van der Putte]

 Hi Vincent,
 Thanks for the repliy; unfortunately 
 gdal_array.BandWriteArray(ds.GetRasterBand(1),a)  
 (BandWriteAsArray doesn't exist)
 has the same effect.
 Tom
 On Fri 28/10/11 10:57 , Vincent Schut schut at sarvision.nl sent:
                    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/f4601ad1/attachment.html