<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=us-ascii">
<META NAME="Generator" CONTENT="MS Exchange Server version 6.5.7654.12">
<TITLE>RE: Compression using the create method in python and aggregation methods</TITLE>
</HEAD>
<BODY>
<!-- Converted from text/rtf format -->
<P><FONT SIZE=2 FACE="Arial">Hi all,</FONT>
</P>
<P><FONT SIZE=2 FACE="Arial">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 </FONT><A HREF="http://www.gdal.org/gdal_tutorial.html"><U><FONT COLOR="#0000FF" SIZE=2 FACE="Arial">http://www.gdal.org/gdal_tutorial.html</FONT></U></A><FONT SIZE=2 FACE="Arial">. 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:</FONT></P>
<P><FONT SIZE=2 FACE="Arial">gdal_translate -co 'COMPRESS=LZW' outfile.tif newoutfile.tif</FONT>
</P>
<P><FONT SIZE=2 FACE="Arial">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?</FONT></P>
<P><FONT SIZE=2 FACE="Arial">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 - </FONT><A HREF="http://sites.google.com/site/spatialpython/aggregating-data-to-grid-cells"><U><FONT COLOR="#0000FF" SIZE=2 FACE="Arial">http://sites.google.com/site/spatialpython/aggregating-data-to-grid-cells</FONT></U></A><FONT SIZE=2 FACE="Arial">). 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?</FONT></P>
<P><FONT SIZE=2 FACE="Arial">Thanks very much in advance for any help you may be able to offer me!</FONT>
</P>
<P><FONT SIZE=2 FACE="Arial">Kind regards,</FONT>
<BR><FONT SIZE=2 FACE="Arial">Andy</FONT>
</P>
<P><FONT SIZE=2 FACE="Arial">s = (640,640)</FONT>
<BR><FONT SIZE=2 FACE="Arial">dt = numpy.dtype('uint16') </FONT>
<BR><FONT SIZE=2 FACE="Arial"># reftile is approx 1km resolution raster, with a unique ID for each cell for a 5 degree square window</FONT>
<BR><FONT SIZE=2 FACE="Arial">gscl = gdal.Open (reftile)</FONT>
<BR><FONT SIZE=2 FACE="Arial">tilescl = g.GetRasterBand(1).ReadAsArray().astype(numpy.uint16) </FONT>
<BR><FONT SIZE=2 FACE="Arial"># reftile90 is a 90m resample (using nearest neighbour) of reftile</FONT>
<BR><FONT SIZE=2 FACE="Arial">g90 = gdal.Open (reftile90)</FONT>
<BR><FONT SIZE=2 FACE="Arial">tile90 = g.GetRasterBand(1).ReadAsArray().astype(numpy.uint16) </FONT>
<BR>
<BR><FONT SIZE=2 FACE="Arial">z = numpy.zeros(s, dtype=dt)</FONT>
<BR><FONT SIZE=2 FACE="Arial">U = unique(tile90[numpy.greater(rec90, 0)])</FONT>
<BR><FONT SIZE=2 FACE="Arial">lenU = len(U)</FONT>
<BR><FONT SIZE=2 FACE="Arial"># for each 90m cell with data, aggregate and write to low resolution output grid</FONT>
<BR><FONT SIZE=2 FACE="Arial">for u in range(lenU):</FONT>
<BR> <FONT SIZE=2 FACE="Arial">result = numpy.sum(rec90[numpy.equal(tile90, U[u])])</FONT>
<BR> <FONT SIZE=2 FACE="Arial">z[numpy.equal(tilescl, U[u])] = result</FONT>
<BR><FONT SIZE=2 FACE="Arial"># Write out the grid</FONT>
<BR><FONT SIZE=2 FACE="Arial">outDrv = gdal.GetDriverByName('GTiff')</FONT>
<BR><B><FONT SIZE=2 FACE="Arial">out = outDrv.Create("outfile.tif", gscl.RasterXSize, gscl.RasterYSize, 1, gdalconst.GDT_UInt16, [ 'COMPRESS=LZW' ] )</FONT></B>
<BR><FONT SIZE=2 FACE="Arial">out.SetProjection(gscl.GetProjection())</FONT>
<BR><FONT SIZE=2 FACE="Arial">out.SetGeoTransform(gscl.GetGeoTransform())</FONT>
<BR><FONT SIZE=2 FACE="Arial">out.GetRasterBand(1).WriteArray(z)</FONT>
<BR><FONT SIZE=2 FACE="Arial">gscl = None</FONT>
<BR><FONT SIZE=2 FACE="Arial">G90 = None</FONT>
<BR><FONT SIZE=2 FACE="Arial">out = None</FONT>
</P>
<BR>
<BR>
<P><FONT SIZE=2 FACE="Arial">--</FONT>
<BR><FONT SIZE=2 FACE="Arial">Andrew Hartley Climate Impacts Risk Analyst</FONT>
<BR><B><FONT SIZE=2 FACE="Arial">Met Office Hadley Centre</FONT></B><FONT SIZE=2 FACE="Arial"> FitzRoy Road Exeter Devon EX1 3PB United Kingdom</FONT>
<BR><FONT SIZE=2 FACE="Arial">Tel: +44 (0)1392 885720 Fax: +44 (0)1392 885681</FONT>
<BR><FONT SIZE=2 FACE="Arial">Email: andrew.hartley@metoffice.gov.uk Website: </FONT><A HREF="file://www.metoffice.gov.uk"><U><FONT COLOR="#0000FF" SIZE=2 FACE="Arial">www.metoffice.gov.uk</FONT></U></A>
</P>
<P><FONT SIZE=2 FACE="Arial">See our guide to climate change at </FONT><A HREF="http://www.metoffice.gov.uk/climatechange/guide/"><U><FONT COLOR="#0000FF" SIZE=2 FACE="Arial">http://www.metoffice.gov.uk/climatechange/guide/</FONT></U></A>
</P>
</BODY>
</HTML>