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

[Hartley, Andrew]

Chaitanya,
Thanks for the help, but it didn't change anything. If it helps, I'm using gdal 1.7.2 with python 3.1 and numpy 1.5.
Cheers,
Andy

________________________________

From: Chaitanya kumar CH [mailto:chaitanya.ch at gmail.com] 
Sent: 22 September 2010 16:11
To: Hartley, Andrew
Cc: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] RE: Compression using the create method in python and aggregation methods


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/e3a4abc9/attachment.html