#!/usr/bin/python import sys, os, tempfile import gdal, osr from gdalconst import * try: import numpy as Numeric Numeric.arrayrange = Numeric.arange except ImportError: import Numeric from argparse import ArgumentParser parser = ArgumentParser( description = "Build an index file from a mapset." ) parser.add_argument( 'files', metavar='file', nargs='+' , help = 'A mapset tile' ) parser.add_argument( 'index', metavar='index' , help = 'The destination index' ) parser.add_argument('-q', '--quiet', action='store_true', help='Be quiet') parser.add_argument('-d', '--debug', action='store_true', help='Debug') args = parser.parse_args() progress = None if not args.quiet: try: progress = gdal.TermProgress_nocb except: progress = gdal.TermProgress # Ensure that GeoTiff is supported format = 'GTiff' dst_driver = gdal.GetDriverByName(format) if dst_driver is None: sys.stderr.write('"%s" driver not registered.\n' % format) sys.exit(1) def debug(msg): if args.debug: sys.stderr.write(msg) def temptif(): (h, f) = tempfile.mkstemp(suffix = '.tif', dir = tmpdir) debug( 'Creating temporary file "%s"\n' % f ) os.close(h) return f # Expand a palette into RBA. This code is mostly taken from pct2rgb(1). def expand_rgba( src_filename ): src_ds = gdal.Open( src_filename ) if src_ds is None: sys.stderr.write('Unable to open "%s".\n' % src_filename) sys.exit(1) if src_ds.RasterCount != 1: return src_filename out_bands = 4 # RGBA src_band = src_ds.GetRasterBand(1) ct = src_band.GetRasterColorTable() if ct is None: return src_filename n = 256 lookup = [ Numeric.arrayrange(n) , Numeric.arrayrange(n) , Numeric.arrayrange(n) , Numeric.ones(n)*(n-1) ] for i in range(min(n,ct.GetCount())): entry = ct.GetColorEntry(i) for c in range(4): lookup[c][i] = entry[c] dst_filename = temptif() dst_ds = gdal.GetDriverByName(format).Create( dst_filename, src_ds.RasterXSize, src_ds.RasterYSize, out_bands ) dst_ds.SetProjection( src_ds.GetProjection() ) dst_ds.SetGeoTransform( src_ds.GetGeoTransform() ) if src_ds.GetGCPCount() > 0: dst_ds.SetGCPs( src_ds.GetGCPs(), src_ds.GetGCPProjection() ) progress( 0.0 ) for iY in range(src_ds.RasterYSize): src_data = src_band.ReadAsArray(0,iY,src_ds.RasterXSize,1) for iBand in range(out_bands): band_lookup = lookup[iBand] dst_data = Numeric.take(band_lookup,src_data) dst_ds.GetRasterBand(iBand+1).WriteArray(dst_data,0,iY) progress( (iY+1.0) / src_ds.RasterYSize ) dst_ds = None return dst_filename # Reduce (and reproject if necessary) the given map. def reduce_map( src_filename ): src_ds = gdal.Open( src_filename ) if src_ds is None: sys.stderr.write('Unable to open "%s".\n' % src_filename) sys.exit(1) src_prj = src_ds.GetProjection() dst_prj = osr.SpatialReference() dst_prj.ImportFromEPSG ( 3067 ) tx = osr.CoordinateTransformation ( src_prj, src_prj ) src_geo = src_ds.GetGeoTransform(); (ulx, uly, ulz ) = tx.TransformPoint( src_geo[0], src_geo[3]) (lrx, lry, lrz ) = tx.TransformPoint( src_geo[0] + src_geo[1]*src_ds.RasterXSize, src_geo[3] + src_geo[5]*src_ds.RasterYSize ) res = 64; # TODO: should be computed (w, h) = ( int((lrx - ulx)/(src_geo[1]*res)) , int((lry - uly)/(src_geo[5]*res)) ) print (w,h); dst_filename = '/tmp/test.tiff' #temptif() dst_ds = gdal.GetDriverByName(format).Create( dst_filename, w, h, src_ds.RasterCount ) dst_geo = ( ulx, src_geo[1]*res, src_geo[2], uly, src_geo[4], src_geo[5]*res ) dst_ds.SetGeoTransform( dst_geo ) dst_ds.SetProjection ( dst_prj.ExportToWkt() ) res = gdal.ReprojectImage( src_ds, dst_ds, src_prj.ExportToWkt(), dst_prj.ExportToWkt(), GRA_Bilinear, 0, 0.0, progress ) if not res: sys.stderr.write('Unable to open "%s".\n' % src_filename) sys.exit(1) dst_ds = None return dst_filename tmpdir = tempfile.mkdtemp() debug( 'Temporary working directory: "%s".\n' % tmpdir ) for src_filename in args.files: expanded = expand_rgba ( src_filename ) reduced = reduce_map ( expanded ) sys.exit(1) src_filename = '/tmp/UV512_RVK_25.tiff' src_proj = osr.SpatialReference() src_proj.SetFromUserInput('epsg:3067') # http://jgomezdans.github.io/gdal_notes/reprojection.html src_filename = '/tmp/UV512_RVK_25.tiff' src_ds = gdal.Open( src_filename ) if src_ds is None: sys.stderr.write('Unable to open "%s".\n' % src_filename) sys.exit(1) gtiff_driver = gdal.GetDriverByName( 'GTiff' ) tif_ds = gtiff_driver.Create( '/tmp/zap.tiff', 300, 150, 4, GDT_Byte) #tx = osr.CoordinateTransformation ( src_ds.GetProjection(), src_ds.GetProjection() ) geo_t = src_ds.GetGeoTransform () x_size = src_ds.RasterXSize # Raster xsize y_size = src_ds.RasterYSize # Raster ysize # Work out the boundaries of the new dataset in the target projection (ulx, uly, ulz ) = tx.TransformPoint( geo_t[0], geo_t[3]) (lrx, lry, lrz ) = tx.TransformPoint( geo_t[0] + geo_t[1]*x_size, \ geo_t[3] + geo_t[5]*y_size ) new_geo = ( ulx, 147, geo_t[2], uly, geo_t[4], -147 ) tif_ds.SetGeoTransform( new_geo ) tif_ds.SetProjection ( osng.ExportToWkt() ) gdal.ReprojectImage(src_ds, tif_ds, wgs84.ExportToWkt(), osng.ExportToWkt(), GRA_Bilinear) tif_ds = None #dataset = gdal.Open( filename, GA_ReadOnly ) # #print 'Driver: ', dataset.GetDriver().ShortName,'/', \ # dataset.GetDriver().LongName #print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize, \ # 'x',dataset.RasterCount # #def get_spatialref(prj): # spatialref = osr.SpatialReference() # spatialref.SetFromUserInput(prj) # return spatialref # #sp1 = get_spatialref('epsg:3067') #sp2 = get_spatialref('+proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs') # #print sp1.IsSameGeogCS(sp2)