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:[email protected]] Sent: 22 September 2010 16:11 To: Hartley, Andrew Cc: [email protected] 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 <[email protected]> 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: [email protected] 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 [email protected] 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
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
