<html>
  <head>
    <meta content="text/html; charset=UTF-8" http-equiv="Content-Type">
  </head>
  <body text="#000000" bgcolor="#FFFFFF">
    Tom van der Putte wrote:
    <blockquote cite="mid:4335.1319790184@vdputte.nl" type="cite">
      <style> BODY { font-family:Arial, Helvetica, sans-serif;font-size:12px; }</style>Hi
      List,<br>
      <br>
      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.<br>
      <br>
      I start with a ndarray, dtype = int32, shape = ( 10,10). It has
      integer values between 1 and 10. The code is as follows:<br>
      <br>
      import numpy as np<br>
      from osgeo import gdal<br>
      from osgeo import osr<br>
      <br>
      #np array<br>
      a = np.ones((10,10))<br>
      <br>
      #coord-system<br>
      sr_str = 'LOCAL_CS["arbitrary"]'<br>
      sr = osr.SpatialReference( sr_str )<br>
      <br>
      #driver<br>
      driver = gdal.GetDriverByName("GTiff")<br>
      <br>
      #dataset<br>
      ds = driver.Create("D:\\temp\\output.tif", 10, 10, 1,
      gdal.GDT_Byte)<br>
      ds.SetProjection( sr_str )<br>
      ds.SetGeoTransform((0, 1, 0,0, 0,1))<br>
      <br>
      #write and close<br>
      ds.GetRasterBand(1).WriteArray(a)<br>
      ds = None<br>
    </blockquote>
    <br>
    Tom,<br>
    <br>
    don't know if you can use WriteArray this way... what I always do,
    is:<br>
    <br>
    from osgeo import gdal_array<br>
    # to write numpy array:<br>
    gdal_array.BandWriteAsArray(ds.GetRasterBand(1), a)<br>
    <br>
    Best,<br>
    Vincent.<br>
    <blockquote cite="mid:4335.1319790184@vdputte.nl" type="cite">
      <br>
      The corresponding gdalinfo is:<br>
      *********<br>
      Raster dataset parameters:<br>
        Projection:
      LOCAL_CS["arbitrary",UNIT["metre",1,AUTHORITY["EPSG","9001"]]]<br>
        RasterCount: 1<br>
        RasterSize (10,10)<br>
      Using driver GeoTIFF<br>
        Metadata:<br>
          0:  AREA_OR_POINT=Area<br>
      <br>
        Image Structure Metadata:<br>
          0:  INTERLEAVE=BAND<br>
      <br>
      Corner Coordinates:<br>
        Upper Left (0, 0)<br>
        Lower Left (0, 10)<br>
        Upper Right (10, 0)<br>
        Lower Right (10, 10)<br>
        Center (5, 5)<br>
      <br>
      Coordinate System is:<br>
      LOCAL_CS["arbitrary",<br>
          UNIT["metre",1,<br>
              AUTHORITY["EPSG","9001"]]]<br>
      Band 1 :<br>
         DataType: Byte<br>
         ColorInterpretation: Gray<br>
         Description:<br>
         Size (10,10)<br>
         BlockSize (10,10)<br>
      *********<br>
      <br>
      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.<br>
      <br>
      I tried to make a Surfer grid (GSBG), and this somewhat worked: it
      gave a nice picture, but the values were off the charts.<br>
      <br>
      Any clues as to what I'm missing here? Cheers,<br>
      <br>
      Tom van der Putte<br>
      <br>
      <br>
      <fieldset class="mimeAttachmentHeader"></fieldset>
      <br>
      <pre wrap="">_______________________________________________
gdal-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a></pre>
    </blockquote>
  </body>
</html>