[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