[gdal-dev] RE: Compression using the create method in python and aggregation methods

[Chaitanya kumar CH]

Andy,

Try this instead.
out = outDrv.Create("outfile.tif", gscl.RasterXSize, gscl.RasterYSize, 1,
gdalconst.GDT_UInt16,  *options = *[ 'COMPRESS=LZW' ] )

However, I am not sure this makes any difference because gdalinfo did tell
that there is LZW compression.

On Wed, Sep 22, 2010 at 4:41 PM, Hartley, Andrew <
andrew.hartley at metoffice.gov.uk> wrote:

>  Hi all,
>
> I'm trying to create a Gtiff with LZW compression using python, with the
> code below, which I wrote with help from the tutorial at *
> http://www.gdal.org/gdal_tutorial.html*<http://www.gdal.org/gdal_tutorial.html>.
> Gdalinfo tells me that the resulting tif ("outfile.tif") has compression
> (Image Structure Metadata: COMPRESSION=LZW), but all my outfiles have the
> same file size and are quite large, so it seems they actually aren't
> compressed. In fact, when I tried:
>
> gdal_translate -co 'COMPRESS=LZW' outfile.tif newoutfile.tif
>
> the newoutfile.tif is considerably smaller than the original file. So, this
> leads me to think that there's a problem with my create() statement below
> (see the line in bold below). Could somebody please tell me what I have
> missed?
>
> Since I'm here, I think I will also pick your brains about cell aggregation
> methods. You'll see from the code below that I have writen a loop to
> aggregate only cells with data. I spent a bit of time considering a few
> options (for example, the excellent pages by Dr Gomez-Dans - *
> http://sites.google.com/site/spatialpython/aggregating-data-to-grid-cells*<http://sites.google.com/site/spatialpython/aggregating-data-to-grid-cells>).
> My code works reasonably well, but since I have lots of processing to do and
> it is not as fast as I would like, I was wondering if anybody could suggest
> a more efficient solution?
>
> Thanks very much in advance for any help you may be able to offer me!
>
> Kind regards,
> Andy
>
> s = (640,640)
> dt = numpy.dtype('uint16')
> # reftile is approx 1km resolution raster, with a unique ID for each cell
> for a 5 degree square window
> gscl = gdal.Open (reftile)
> tilescl = g.GetRasterBand(1).ReadAsArray().astype(numpy.uint16)
> # reftile90 is a 90m resample (using nearest neighbour) of reftile
> g90 = gdal.Open (reftile90)
> tile90 = g.GetRasterBand(1).ReadAsArray().astype(numpy.uint16)
>
> z = numpy.zeros(s, dtype=dt)
> U = unique(tile90[numpy.greater(rec90, 0)])
> lenU = len(U)
> # for each 90m cell with data, aggregate and write to low resolution output
> grid
> for u in range(lenU):
>         result = numpy.sum(rec90[numpy.equal(tile90, U[u])])
>         z[numpy.equal(tilescl, U[u])] = result
> # Write out the grid
> outDrv = gdal.GetDriverByName('GTiff')
> *out = outDrv.Create("outfile.tif", gscl.RasterXSize, gscl.RasterYSize, 1,
> gdalconst.GDT_UInt16,  [ 'COMPRESS=LZW' ] )*
> out.SetProjection(gscl.GetProjection())
> out.SetGeoTransform(gscl.GetGeoTransform())
> out.GetRasterBand(1).WriteArray(z)
> gscl = None
> G90 = None
> out = None
>
>
> --
> Andrew Hartley  Climate Impacts Risk Analyst
> *Met Office Hadley Centre*  FitzRoy Road  Exeter  Devon  EX1 3PB  United
> Kingdom
> Tel: +44 (0)1392 885720  Fax: +44 (0)1392 885681
> Email: andrew.hartley at metoffice.gov.uk  Website: *www.metoffice.gov.uk*
>
> See our guide to climate change at *
> http://www.metoffice.gov.uk/climatechange/guide/*<http://www.metoffice.gov.uk/climatechange/guide/>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>



-- 
Best regards,
Chaitanya kumar CH.
/tʃaɪθənjə/ /kʊmɑr/
+91-9494447584
17.2416N 80.1426E
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100922/01a0afd3/attachment-0001.html