Andy,<br><br>Try this instead.<br><font style="font-family: courier new,monospace;" size="2">out = outDrv.Create("outfile.tif", gscl.RasterXSize, gscl.RasterYSize, 1, gdalconst.GDT_UInt16, </font><b style="font-family: courier new,monospace;"><font size="2">options = </font></b><font style="font-family: courier new,monospace;" size="2">[ 'COMPRESS=LZW' ] )</font>
<br><br>However, I am not sure this makes any difference because gdalinfo did tell that there is LZW compression.<br><br><div class="gmail_quote">On Wed, Sep 22, 2010 at 4:41 PM, Hartley, Andrew <span dir="ltr"><<a href="mailto:andrew.hartley@metoffice.gov.uk">andrew.hartley@metoffice.gov.uk</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<div>
<p><font face="Arial" size="2">Hi all,</font>
</p>
<p><font face="Arial" size="2">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" target="_blank"><u><font color="#0000ff" face="Arial" size="2">http://www.gdal.org/gdal_tutorial.html</font></u></a><font face="Arial" size="2">. 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 face="Arial" size="2">gdal_translate -co 'COMPRESS=LZW' outfile.tif newoutfile.tif</font>
</p>
<p><font face="Arial" size="2">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 face="Arial" size="2">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" target="_blank"><u><font color="#0000ff" face="Arial" size="2">http://sites.google.com/site/spatialpython/aggregating-data-to-grid-cells</font></u></a><font face="Arial" size="2">). 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 face="Arial" size="2">Thanks very much in advance for any help you may be able to offer me!</font>
</p>
<p><font face="Arial" size="2">Kind regards,</font>
<br><font face="Arial" size="2">Andy</font>
</p>
<p><font face="Arial" size="2">s = (640,640)</font>
<br><font face="Arial" size="2">dt = numpy.dtype('uint16') </font>
<br><font face="Arial" size="2"># reftile is approx 1km resolution raster, with a unique ID for each cell for a 5 degree square window</font>
<br><font face="Arial" size="2">gscl = gdal.Open (reftile)</font>
<br><font face="Arial" size="2">tilescl = g.GetRasterBand(1).ReadAsArray().astype(numpy.uint16) </font>
<br><font face="Arial" size="2"># reftile90 is a 90m resample (using nearest neighbour) of reftile</font>
<br><font face="Arial" size="2">g90 = gdal.Open (reftile90)</font>
<br><font face="Arial" size="2">tile90 = g.GetRasterBand(1).ReadAsArray().astype(numpy.uint16) </font>
<br>
<br><font face="Arial" size="2">z = numpy.zeros(s, dtype=dt)</font>
<br><font face="Arial" size="2">U = unique(tile90[numpy.greater(rec90, 0)])</font>
<br><font face="Arial" size="2">lenU = len(U)</font>
<br><font face="Arial" size="2"># for each 90m cell with data, aggregate and write to low resolution output grid</font>
<br><font face="Arial" size="2">for u in range(lenU):</font>
<br> <font face="Arial" size="2">result = numpy.sum(rec90[numpy.equal(tile90, U[u])])</font>
<br> <font face="Arial" size="2">z[numpy.equal(tilescl, U[u])] = result</font>
<br><font face="Arial" size="2"># Write out the grid</font>
<br><font face="Arial" size="2">outDrv = gdal.GetDriverByName('GTiff')</font>
<br><b><font face="Arial" size="2">out = outDrv.Create("outfile.tif", gscl.RasterXSize, gscl.RasterYSize, 1, gdalconst.GDT_UInt16, [ 'COMPRESS=LZW' ] )</font></b>
<br><font face="Arial" size="2">out.SetProjection(gscl.GetProjection())</font>
<br><font face="Arial" size="2">out.SetGeoTransform(gscl.GetGeoTransform())</font>
<br><font face="Arial" size="2">out.GetRasterBand(1).WriteArray(z)</font>
<br><font face="Arial" size="2">gscl = None</font>
<br><font face="Arial" size="2">G90 = None</font>
<br><font face="Arial" size="2">out = None</font>
</p>
<br><font color="#888888">
<br>
<p><font face="Arial" size="2">--</font>
<br><font face="Arial" size="2">Andrew Hartley Climate Impacts Risk Analyst</font>
<br><b><font face="Arial" size="2">Met Office Hadley Centre</font></b><font face="Arial" size="2"> FitzRoy Road Exeter Devon EX1 3PB United Kingdom</font>
<br><font face="Arial" size="2">Tel: +44 (0)1392 885720 Fax: +44 (0)1392 885681</font>
<br><font face="Arial" size="2">Email: <a href="mailto:andrew.hartley@metoffice.gov.uk" target="_blank">andrew.hartley@metoffice.gov.uk</a> Website: </font><a><u><font color="#0000ff" face="Arial" size="2">www.metoffice.gov.uk</font></u></a>
</p>
<p><font face="Arial" size="2">See our guide to climate change at </font><a href="http://www.metoffice.gov.uk/climatechange/guide/" target="_blank"><u><font color="#0000ff" face="Arial" size="2">http://www.metoffice.gov.uk/climatechange/guide/</font></u></a>
</p>
</font></div>
<br>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br></blockquote></div><br><br clear="all"><br>-- <br>Best regards,<br>Chaitanya kumar CH.<br>
/tʃaɪθənjə/ /kʊmɑr/ <br>+91-9494447584<br>17.2416N 80.1426E<br>