diff options
| -rwxr-xr-x | mkindex.py | 198 |
1 files changed, 0 insertions, 198 deletions
diff --git a/mkindex.py b/mkindex.py deleted file mode 100755 index 46c5af1..0000000 --- a/mkindex.py +++ /dev/null @@ -1,198 +0,0 @@ -#!/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) |