[gdal-dev] numpy array to raster

questions anon questions.anon at gmail.com
Wed Feb 1 16:50:08 EST 2012


Thanks for responding Anton, I have been playing with the datatype but so
far it hasn't fixed the problem. Also I do have access to ArcGIS to check
the output.
I will keep trying, this is where I am at currently:

from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr

myarray=myarray.astype(N.float32)
xmin,ymin,xmax,ymax=[139.8,-39.2,150.0,-33.6]
ncols,nrows=[193,106]
xres=(xmax-xmin)/float(ncols)
yres=(ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
dst_ds =
gdal.GetDriverByName('GTiff').Create('E:/test/rasterise/myraster.tif',ncols,
nrows, 1 ,gdal.GDT_Float32)
dst_ds.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) #WGS84 lat long.
dst_ds.SetProjection( srs.ExportToWkt() )
dst_ds.GetRasterBand(1).WriteArray(myarray)

dst_ds=None


On Wed, Feb 1, 2012 at 7:12 PM, Anton Korosov <anton.korosov at nersc.no>wrote:

> Probably the problem is with data type. You obviously have data of float
> type and you give GDT_Byte to the gdal.GetDriverByName('GTiff').**
> Create().
> Either try to create array with Byte data or give GDT_Float32. In the
> latter case the produced geotiff will have 32 bits fo each pixel. It is not
> supported by simple viewers but should be read with ArcGIS.
>
>
> On 02/01/2012 03:31 AM, questions anon wrote:
>
>> no I quickly checked that!
>> myarray is:
>> [[ 16.15035553  16.14380074  16.15581551 ...,  18.06388149  18.08930645
>>    18.08825245]
>>  [ 16.2154911   16.21180592  16.23977184 ...,  18.1085537   18.12040272
>>    18.12342682]
>>  [ 16.32851467  16.29202938  16.28964043 ...,  18.16753635  18.14905453
>>    18.16632977]
>>  ...,
>>  [ 30.50812759  30.58384018  30.66707535 ...,  26.31020527  26.47789459
>>    27.11495361]
>>  [ 30.76499577  30.76497536  30.79138317 ...,  26.45928288  26.64059887
>>    27.03641129]
>>  [ 31.01263275  30.96269417  30.9933857  ...,  26.78247185  26.77845631
>>    26.97975636]]
>>
>> On Wed, Feb 1, 2012 at 1:25 PM, Kyle Shannon <ksshannon at gmail.com
>> <mailto:ksshannon at gmail.com>> wrote:
>>
>>    Is 'myarray' full of zeros?
>>
>>    /**
>>      *
>>      * Kyle Shannon
>>      * ksshannon at gmail.com <mailto:ksshannon at gmail.com>
>>
>>      *
>>      */
>>
>>
>>
>>    On Tue, Jan 31, 2012 at 19:20, questions anon
>>    <questions.anon at gmail.com <mailto:questions.anon at gmail.**com<questions.anon at gmail.com>>>
>> wrote:
>>
>>        thanks Frank, following your instructions with:
>>
>>        src_ds=gdal_array.OpenArray(**myarray)
>>        dst_ds =
>>        gdal.GetDriverByName('GTiff').**Create('E:/test/rasterise/**
>> mynewraster.tif',ncols,
>>        nrows, 1 ,gdal.GDT_Byte)
>>        dst_ds.SetGeoTransform(**geotransform)
>>        dst_ds.GetRasterBand(1).**WriteArray(myarray)
>>
>>        I do not receive any error messages but the tiff produced are
>>        all just zeros. Is there a step I am missing?
>>        thanks
>>
>>
>>        On Wed, Feb 1, 2012 at 12:07 PM, Frank Warmerdam
>>        <warmerdam at pobox.com <mailto:warmerdam at pobox.com>> wrote:
>>
>>            On Tue, Jan 31, 2012 at 4:38 PM, questions anon
>>            <questions.anon at gmail.com <mailto:questions.anon at gmail.**com<questions.anon at gmail.com>
>> >>
>>
>>            wrote:
>>             > I need to output my numpy array as a raster so that
>>            someone else can access
>>             > the data in ArcGIS. So basically the steps I need are:
>>             > read numpy array into gdal
>>             > convert to raster
>>             > use latitude and longitude and array size to set projection
>>             >
>>             > I am really struggling with gdal because I can't seem to
>>            find enough
>>             > documentation about each step to understand what it is
>> doing.
>>             > Here are some of the steps I think I need:
>>             >
>>             > myarray=myarray
>>             > #the extent and shape of my array
>>             > xmin,ymin,xmax,ymax=[139.8,-**39.2,150.0,-33.6]
>>             > ncols,nrows=[193,106]
>>             > xres=(xmax-xmin)/float(ncols)
>>             > yres=(ymax-ymin)/float(nrows)
>>             > geotransform=(xmin,xres,0,**ymax,0, -yres)
>>             >
>>             > from osgeo import gdal
>>             > from osgeo import gdal_array
>>             >
>>             > src_ds=gdal_array.OpenArray(**myarray)
>>             >
>>             > dst_ds =
>>             >
>>            gdal.GetDriverByName('GTiff').**Create('E:/test/rasterise/**
>> mynewraster.tif',ncols,
>>             > nrows, 1 ,gdal.GDT_Byte)
>>             > dst_rb = dst_ds.GetRasterBand(0)
>>             > dst_ds.SetGeoTransform(**geotransform)
>>             > output = gdal.RasterizeLayer(dst_ds)
>>
>>            Dear "Questions Anon",
>>
>>            There are a variety of Python related information at:
>>
>>            http://trac.osgeo.org/gdal/**wiki/GdalOgrInPython<http://trac.osgeo.org/gdal/wiki/GdalOgrInPython>
>>
>>            One serious issue with the above is that
>>            "gdal.RasterizeLayer()" is used to
>>            turn vector data into raster data - for instance to
>>            rasterize polygon features
>>            into an existing raster file.  I think you want to write
>>            your array into a
>>            raster file.  I think you can just call
>>            "dst_ds.WriteArray(myarray)".
>>
>>            Alternatively if that method does not exist on the dataset, you
>>            can write the one band like:
>>
>>              dst_ds.GetRasterBand(1).**WriteArray(myarray)
>>
>>            Some of the samples referenced from the GdalOgrInPython
>>            should be
>>            helpful though I understand it can be hard to know where to
>>            look.
>>
>>            Best regards,
>>            --
>>            ------------------------------**---------+--------------------
>> **------------------
>>            I set the clouds in motion - turn up   | Frank Warmerdam,
>>            warmerdam at pobox.com <mailto:warmerdam at pobox.com>
>>
>>            light and sound - activate the windows |
>>            http://pobox.com/~warmerdam <http://pobox.com/%7Ewarmerdam**>
>>
>>            and watch the world go round - Rush    | Geospatial Software
>>            Developer
>>
>>
>>
>>        ______________________________**_________________
>>        gdal-dev mailing list
>>        gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.**org<gdal-dev at lists.osgeo.org>
>> >
>>        http://lists.osgeo.org/**mailman/listinfo/gdal-dev<http://lists.osgeo.org/mailman/listinfo/gdal-dev>
>>
>>
>>
>>
>>
>>
>> ______________________________**_________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/**mailman/listinfo/gdal-dev<http://lists.osgeo.org/mailman/listinfo/gdal-dev>
>>
> ______________________________**_________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/**mailman/listinfo/gdal-dev<http://lists.osgeo.org/mailman/listinfo/gdal-dev>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20120202/7a45f9a0/attachment.html


More information about the gdal-dev mailing list