#!/usr/bin/env python

# topomap: use data files from Topo! CD to generate a map at
# the specified coordinates.

# Copyright 2005 by Akkana Peck.
# You may use, distribute or modify this program under the terms of the GPL.

# Topo! maps are USGS quadrants scanned to 262x328 GIF files.
# They live in filenames like: q37122c2/012t0501.gif
#
# The directory name corresponds to one 7.5-minute USGS map.
# Starting at 37 degrees latitude, 122 degrees longitude,
# it's the third map up (c) and second over (2) at 7.5 minutes each;
# in other words, 122^07'30" 37^15' in the lower right corner
# (coordinates increase going left and up).
#
# The filename represents a piece of the map in that directory:
# 012t means the 7.5 minute series, 024t is the 15-minute series.
# 0503 is the fifth map over, third map down, out of ten in each direction.
# Hence each map for the 7.5-minute series spans .75 minutes.

# We will use coordinates specified as deg.decimal_minutes

import sys, os

def Usage() :
    print "Usage: topomap start_lat start_long [long_width lat_height]"
    print "Run from inside the Topo! data directory (e.g. sfr_data on the CD)."
    print "Use degrees.decimal_minutes format for coordinates."
    sys.exit(1)

# Good grief: python thinks 15.0 - 15 is 5.68434188608e-13
# and entering .nn30 gets rounded to 29 seconds.
# (Discussion at http://www.network-theory.co.uk/docs/pytut/tut_86.html
# of binary floating point arithmetic precision.)
# So if we want the seconds right, every time we need an integer
# truncation, we'd better add a little bit of slop to make sure
# python isn't erring on the low side.
def intTrunc(num) :
    return int(num + .00001)

# Translate .mmss into decimal minutes.
# def dotMMSStoDecimalMinutes(dotMMSS) :
#     dotMMSS = dotMMSS * 100
#     decMin = intTrunc(dotMMSS)     # truncated minutes
#     dotMMSS = dotMMSS - decMin
#     sec = intTrunc(dotMMSS * 100)
#     decMin = decMin + sec / 60.
#     return decMin

# Add angles in deg.mmss
def AddAngle(angle1, angle2) :
    #print "AddAngle(", angle1, angle2, ")"
    ang1 = intTrunc(angle1)
    ang2 = intTrunc(angle2)
    #print "trunc:", ang1, ang2
    minutes1 = (angle1 - ang1) * 5 / 3
    minutes2 = (angle2 - ang2) * 5 / 3
    ang1 = ang1 + ang2 + minutes1 + minutes2
    #print "Decimal sum:", ang1
    # convert back to deg.minutes
    ang2 = intTrunc(ang1)
    minutes1 = (ang1 - ang2) * .6
    return ang2 + minutes1

def ohstring(num) :
    if (num < 10) : return "0" + str(num)
    return str(num)

# Given a pair of coordinates in deg.mmss, map to the closest filename
# e.g. q37122c2/012t0501.gif.
def CoordsToFilename(longitude, latitude, series) :
    latDeg = intTrunc(latitude)
    longDeg = intTrunc(longitude)
    latMin = (latitude - latDeg) * 100
    longMin = (longitude - longDeg) * 100
    longMinOrd = intTrunc(longMin / 7.5)
    latMinOrd = intTrunc(latMin / 7.5)
    #print "longMin =", longMin, "latMin =", latMin
    #print "Ords:", longMinOrd, latMinOrd

    dirname = "q" + str(latDeg) + str(longDeg) + \
              chr(ord('a') + latMinOrd) + str(longMinOrd+1)

    # Find the difference between our desired coordinates
    # and the origin of the map this directory represents.
    latMinDiff = latMin - (latMinOrd * 7.5)
    longMinDiff = longMin - (longMinOrd * 7.5)

    latOffset = intTrunc(latMinDiff * 10 / series)
    longOffset = intTrunc(longMinDiff * 10 / series)

    #print "Diffs:", longMinDiff, latMinDiff
    #print "Offsets:", longOffset, latOffset

    # Now calculate the current filename.
    if (series > 13) :
        fileprefix = "024t"
        numcharts = 5
    else :
        fileprefix = "012t"
        numcharts = 10
    filename = fileprefix + ohstring(numcharts-longOffset) + \
               ohstring(numcharts-latOffset) + ".gif"

    return dirname + "/" + filename

def MakeMap(start_long, start_lat, end_long, end_lat, series) :
    outfile = topodir + "topo" + str(series) + "_" + \
              str(start_long) + "_" + str(start_lat) + "_" + \
              str(end_long) + "_" + str(end_lat) + ".gif"
    nofiles = 1

    cmdstr1 = "montage -geometry 262x328"
    cmdstr2 = ""

    if (series > 13) : dAngle = .015
    else: dAngle = .0075

    ny = 0
    curlat = start_lat
    while curlat <= end_lat :
        curlong = start_long
        nx = 0
        while curlong <= end_long :
            #print "current:", curlat, curlong
            filename = CoordsToFilename(curlong, curlat, series)
            #print "7.5:", CoordsToFilename(curlong, curlat, 7.5)
            #print " 15:", CoordsToFilename(curlong, curlat, 15)
            if os.access(filename, os.R_OK) :
                nofiles = 0
            else :
                print "Can't open", filename
                filename = topodir + "blank.gif"

            cmdstr2 = filename + " " + cmdstr2
            curlong = AddAngle(curlong, dAngle)
            nx += 1
        curlat = AddAngle(curlat, dAngle)
        ny += 1

    if nofiles :
        print "No files!"
        return

    cmdstr1 += " -tile " + str(nx) + "x" + str(ny) + " " + \
               cmdstr2 + " " + outfile
    #print cmdstr1
    os.system(cmdstr1)
    if (os.access(outfile, os.R_OK)) :
        print "Created file:", outfile

# main()

# Parse flags:
series7_5 = 0
series15 = 0
while len(sys.argv) > 3 and sys.argv[1][0] == '-' :
    if sys.argv[1] == "-15" :
        series15 = 1
    elif sys.argv[1][0:2] == "-7" :
        series7_5 = 1
    sys.argv = sys.argv[1:]

if len(sys.argv) < 3 :
    Usage()

# Default is to make both 7.5 minute and 15 minute maps:
if series7_5 + series15 <= 0 :
    series7_5 = 1
    series15 = 1

try :
    start_long = float(sys.argv[1])
    start_lat = float(sys.argv[2])
    if len(sys.argv) > 3 :
        end_long = AddAngle(start_long, float(sys.argv[3]))
    else :
        end_long = start_long + .0625
    if len(sys.argv) > 4 :
        end_lat = AddAngle(end_lat, float(sys.argv[4]))
    else :
        end_lat = start_lat + .0625
except ValueError, e :
    Usage()

print start_long, start_lat, end_long, end_lat

topodir = os.environ["HOME"] + "/Topo/"

if series7_5 :
    MakeMap(start_long, start_lat, end_long, end_lat, 7.5)
if series15 :
    MakeMap(start_long, start_lat, end_long, end_lat, 15)