> 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: [email protected] Website: www.metoffice.gov.uk > > See our guide to climate change at > http://www.metoffice.gov.uk/climatechange/guide/ >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
