From 549390f42b7cce09ad59c6882258b73f593cbbd6 Mon Sep 17 00:00:00 2001 From: Guilhem Moulin Date: Sat, 10 Aug 2013 18:05:20 +0200 Subject: Initial commit --- mkindex.py | 198 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100755 mkindex.py (limited to 'mkindex.py') diff --git a/mkindex.py b/mkindex.py new file mode 100755 index 0000000..46c5af1 --- /dev/null +++ b/mkindex.py @@ -0,0 +1,198 @@ +#!/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) -- cgit v1.2.3