#!/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)