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

[Hartley, Andrew]

> 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. 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-ce
> lls). 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/
> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100922/10fd2e07/attachment.html