summaryrefslogtreecommitdiffstats
path: root/mkindex.py
diff options
context:
space:
mode:
Diffstat (limited to 'mkindex.py')
-rwxr-xr-xmkindex.py198
1 files changed, 198 insertions, 0 deletions
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)