summaryrefslogtreecommitdiffstats
path: root/mkindex.py
diff options
context:
space:
mode:
Diffstat (limited to 'mkindex.py')
-rwxr-xr-xmkindex.py198
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)