PythonScripts: gen_terrain_image_script.py

'''
First step in generating a terrain script of the entire world
Assumes you have the SRTM tiles of the entire world

Instead of loading all tiles and making a map, 
writes a geo-referenced image per slice.

Use gen_tiles_from_images.py to convert these images to web tiles.
Use convert_georef.py to convert Maperitive's georef format to GDAL's

'''
import os
import argparse

#defaults which can be modified via command-line params
images_dir = "G:\\Images"
maperitive_path = "G:\\Maperitive\\"
script_file = "gen_images.mscript"

#These are the full bounds of the SRTM data
map_bounds = {
    'W':-180,
    'E':180,
    'S':-56,
    'N':59,
}

script_head = '''
// SCRIPT TO GENERATE WORLD TERRAIN SHADING GEOREFERENCED IMAGES
// Expects {src} data in Maperitive/Cache/Rasters/{src}
// Will output images to {idir}

set-dem-source name={src}
'''

#To create tiles from these images, uncomment these lines and comment
#the clear-map call in script_gen.
#NOTE: This will only work if the images are small enough to not fill up
#   your RAM! Try a large sample rate, such as 4 or 6.
script_foot = '''
//set-geo-bounds -180,-56,179.999999,60
//generate-tiles exclude-partial=true minzoom=1 maxzoom=7 tilesdir=g:\lowtiles
'''

script_gen = '''
//Tile {tnum} of {ttot}
clear-map
set-geo-bounds {W},{S},{E},{N}
generate-relief-igor color="black" intensity=3 sample-rate={rate}
//generate-contours interval=10
save-source index=1 file={fname}
'''

def gen_bounds(args):
    #generate the maperitive script metadata
    print "Generating script starting at slice", args.index

    #Break the map into chunks, stepping an input number of degrees
    EW_step = args.step #longitude
    NS_step = args.step #latitude

    EW_range = range(map_bounds['W'], map_bounds['E'], EW_step) 
    NS_range = range(map_bounds['S'], map_bounds['N'], NS_step)

    if args.srtm3:
        mem = EW_step * NS_step * (1201*1201*2)
    else:
        mem = EW_step * NS_step * (3601*3601*2)
    total = len(EW_range) * len(NS_range)
    print "Divided world into {} slices ({:,} bytes each)".format(total,mem)

    metadata = []
    num = 1
    for EW in EW_range:
        for NS in NS_range:
            bounds = {
                'tnum': num,
                'ttot': total,
                'W':max(map_bounds['W'], EW),
                'E':min(map_bounds['E'], EW + EW_step),
                'S':max(map_bounds['S'], NS),
                'N':min(map_bounds['N'], NS + NS_step),                
            }

            if num == args.index:
                print "Start tile #{tnum}: {W},{S},{E},{N}".format(**bounds)
            if num >= args.index:
                metadata.append(bounds)

            num += 1
    return metadata

def gen_script(data, args):
    # Generate the maperitive script
    source = "SRTM3" if args.srtm3 else "SRTM1"
    with open(os.path.join(args.mdir, "Scripts", args.script), "w") as F:
        F.write(script_head.format(src=source, idir=args.idir))

        for elem in data:
            Fname = "tile_%05d.%s"%(elem['tnum'], args.format)
            F.write(script_gen.format(fname=os.path.join(args.idir, Fname), rate=args.rate, **elem))

        F.write(script_foot)

def main(args):
    if not os.path.exists(args.idir):
        os.mkdir(args.idir)

    data = gen_bounds(args)
    gen_script(data, args)

if __name__ == '__main__':
    P = argparse.ArgumentParser(
        description="Generates a maperitive script to create a set of georeferenced images of the world's terrain",
        epilog="To create a good set of world data, run with: --srtm3 -s 15 -r 6")
    P.add_argument("-i", "--index", type=int, default=1,
         help="Index to start at (helpful to continue a failed run where you left off)")
    P.add_argument("-s", "--step", type=int, default=15,
         help="Degrees to step per tile")    
    P.add_argument("-r", "--rate", type=int, default=6,
         help="Sample rate of hillshading (higher makes a smaller image)")
    P.add_argument("-f", "--format", choices=["png", "tif", "jpg"], default="png",
        help="Image format to write")
    P.add_argument("--srtm3", action="store_true",
        help="Use the lower-resolution SRTMGL3 data instead of GL1")
    P.add_argument("--mdir", default=maperitive_path, 
        help="Path to the Maperitive folder to write the script to. Default: %s"%maperitive_path)
    P.add_argument("--idir", default=images_dir, 
        help="Path to output the generated tiles. Default: %s"%images_dir)
    P.add_argument("--script", default=script_file, 
        help="File name of the script to generate. Default: %s"%script_file)
    args = P.parse_args()

    main(args)