Tag Archives: Python

Map Tiles with Python + GDAL

After getting GPSDrive running on OS X I started looking into different sources for maps. The GPSDrive wiki has a couple of interesting pages about this: one on creating maps and one on map sources. The second page pointed me to the LibreMap Project which has a complete collection of  USGS 1:24000 (large scale) topographic maps for the entire US. These files are distributed as high-resolution GeoTIFF files, along with world files, so the maps are fully geotagged. This means that they contain all of the information needed to generate GPSDrive-friendly map tiles, but GPSDrive will not read them directly. The creating maps page gives some (heavily out-of-date) advice on how to use the gdal_slice.sh script distributed with GPSDrive to create map tiles using the tools from GDAL. Sadly, this script is completely unusable on OS X, because of incompatibilities in the command line tools that are used by the shell script. After spending about five minutes trying to tune it to work with the tools that Apple ships, I completely gave up and looked into other approaches.

Luckily, I quickly found out that GDAL ships with Python bindings. I installed GDAL from source using the instructions from their wiki. To make sure it built with Python modules I ran ./configure –with-python. After that make and sudo make install ran without a hitch.

I then set about re-implementing gdal_slice.sh in Python. Once I got a handle on how to use the Python bindings it was fairly easy to write. You can find the script here.

Using gdal_slice.py is very similar to using gdal_slice.sh:

$ ./gdal_slice.py -h
Usage: gdal_slice.py [options] FILENAME

Options:
  -h, --help            show this help message and exit
  -o OVERLAP, --overlap=OVERLAP
                        percentage tiles will overlap. should be at least 20%
  -a, --add             write map info to map_koord.txt in current working
                        directory
  -m, --map             use *_map folders for output; use if input image is
                        UTM. Default behavior
  -t, --topo            use *_top folders for output; use if input image is
                        not UTM
  -v, --verbose         enable debugging output

The biggest difference is that gdal_slice.py does not perform any format conversion; it will write maps to TIFF files. Although the file size is decidedly larger than that of PNGs, my version of GDAL is unable to create PNGs, and disk space is (fairly) cheap. If I’m motivated I’ll get around to adding a format option, at least to allow conversion to PNG.

Running gdal_slice.py filename.ext will create a folder named filename_map (or _top, if you use the -t flag) that contains a set of 1280×1024 TIFF files, as well as a map_koord_draft.txt file. Moving this folder into GPSDrive’s map directory, and merging filename_[map|top]/map_koord_draft.txt with map_koord.txt will make the new maps available to GPSDrive. Note that filename.ext must be in the current working directory, as the script handles filenames/paths somewhat naivëly.

So far I’ve only tested this with GeoTIFFs from LibreMaps, but it seems to work great, with the caveat that LibreMaps’ files include the borders of map, and my script doesn’t do any cropping, so you see borders if you’re near the edges of quads. This doesn’t bother me terribly much, but you may need to turn mapsets on and off in the Map Control dialog to see the appropriate images.

Update 7-14-11: I’ve modified the script with an additional flag (-a) that will automatically merge the generated map_koords file with map_koord.txt in the current working directory. This means that if you run the script from inside .gpsdrive/maps with the -a flag the generated tiles will be automatically added to GPSDrive’s database. The script avoids creating duplicate entries in map_koord.txt, as well as alphabetizing all entries.

I’ve modified the documentation here, as well as the links, to reference the latest version.

Sign your email like Randall Waterhouse

Neal Stephenson’s novel Cryptonomicon is a fascinating book; I’m re-reading it for the umpteenth time and loving it every bit as much as when I first read it in 7th grade. One thing that caught my attention on this reading (especially in light of my recent GPS post) is the way that one of the main characters, Randall Lawrence Waterhouse, signs his email:

Randall Lawrence Waterhouse

Current meatspace coordinates, hot from the GPS receiver card in my laptop:

8 degrees, 52.33 minutes N latitude 117 degrees, 42.75 minutes E longitude

Nearest geographical feature: Palawan, the Philippines

(from http://www.euskalnet.net/larraorma/crypto/slide51.html)

Using GPSd’s python bindings, and the free, CC-licenced RESTful API from GeoNames.org, I wrote a short python script that will produce a similar output:

Seth Just

Current meatspace coordinates, hot from the GPS receiver card in my laptop:

41.751 N latitude, 111.807 W longitude

Nearest geographical feature: Logan Canyon, Utah, United States

The core of the script is two functions — one that gets GPS coordinates from GPSd, and a second that does a lookup for geographical features (mountains, canyons, islands, etc) on GeoNames.org. For more details, see the final section of this post.

Using gpssig.py

To use gpssig.py, first download it from here, and change the extension to .py.

In order to be useful as an email signature, I had to produce output in HTML format, which allows almost any mail client to use the output as a signature. I’ve gone through the effort of setting it up with the two clients I use the most: Mozilla Thunderbird and Apple’s Mail.app. Of these two, Thunderbird has documented support for HTML signatures (see here for details), and should be the most like other email clients. This script has been tested on Mac OS X Snow Leopard, and should run on Linux without problems (although I haven’t tried it on versions of python higher than 2.6.1).

Using gpssig.py takes two steps: first configure your Mail client to take a signature from a file, and then configure gpssig.py to run automatically to update that file.

Thunderbird is relatively easy to set up with an HTML signature. Because it only supports one signature per account, you simply need to tell it to draw that signature from a file. Under each account’s settings there’s an option to “Attach the signature from a file” (see here for a detailed howto). I chose to use ~/.signature.html, but you can use any file you’d like to, just remember the location you choose.

Mail.app is a more finicky beast — although the script supports it, I recommend using Thunderbird, mainly because Mail only refreshes its signatures from files on every launch, whereas Thunderbird will read the file each time you compose a message. However, if that limitation is acceptable (if you don’t move very much, perhaps), it will work just fine. The main issue with Mail.app is that it doesn’t store signatures in an HTML file, but instead in the proprietary .webarchive format. These files can be found in ~/Library/Mail/Signatures. To use gpssig.py with Mail, you need to follow steps 4&5 from http://bytes.com/topic/macosx/insights/825900-setup-html-signatures-apple-mail-mail-app. First create a new signature (under Preferences->Signatures). This will create a new file in the signature folder — take note of the file name. You will use this file name to configure gpssig.py.

gpssig.py is configured with command line options:

Usage: gpssig.py [options] NAME FILENAME
Note: NAME should be enclosed in quotes, unless it contains no spaces.

Options:
  -h, --help   show this help message and exit
  --no-gps     don't attempt to get coordinates from GPS
  -a, --amail  generate files for Apple's Mail.app (.webarchive)
  --lat=LAT    default latitude, if GPS is unavailable
  --lon=LON    default longitude, if GPS is unavailable

The simplest way to use it is manually — simply running it whenever you want to update your signature. For example, I use ./gpssig.py –lat=41.752 –lon=-111.793 ‘Seth Just’ ~/.signature.html to generate a signature for Thunderbird.

If you want to make this automatic, you’ll need to have the script run automatically. On OS X you’ll need to use launchd, while on linux you need cron. Both of these make setting up repeating actions very easy — see here for information on launchd and here for information on cron.

Sadly I haven’t yet figured out a way to run the script whenever you write an email — if I figure it out, I’ll be sure to post something.

The Script

gpssig.py is built around three major functions — one gets coordinates from GPSd, one gets geographical feature information from

The function I use to get GPS information is a bit hacky — there doesn’t seem to be terribly much documentation on GPSd’s python bindings, so I made do with what I could figure out:

class NoGPSError (Exception): pass

def get_lat_lon():
  # GPSd Python bindings
  import gps

  # Create GPS object
  session = gps.gps()

  # Set GPS object as an iterator over reports from GPSd
  session.stream(gps.WATCH_ENABLE|gps.WATCH_NEWSTYLE)

  # Loop until we get a report with lat and lon. The limit of 5 loops should be more than enough -- my gps never takes more than three when it has a lock
  i = 0
  while (1):
    try:
      session.next()
      lat, lon = session.data['lat'], session.data['lon']
      break
    except:
      if (i > 5): raise NoGPSError 
      i += 1

  return lat, lon

The other function is more straightforward — it makes two http requests to GeoNames.org to get the nearest geographical feature and locality. It then combines those responses in a way that should provide a good response almost anywhere in the world (provided that GeoNames has good, local data). I’ve tested it with a variety of coordinates (provided by http://www.getlatlon.com/), and the responses seem to be pretty good.

def get_geo_string(lat, lon):
  # Modules for REST/XML request from GeoNames
  import urllib
  from xml.etree import ElementTree as ET
  
  # Request the name of the nearest geographical feature
  placename = ET.parse(
      urllib.urlopen("http://api.geonames.org/findNearby?lat=" + str(lat) + "&lng=" + str(lon) + "&featureClass=T&username=sethjust")
      )
  subdiv = ET.parse(
      urllib.urlopen("http://api.geonames.org/countrySubdivision?lat=" + str(lat) + "&lng=" + str(lon) + "&username=sethjust")
      )
  
  # Format and return place, region, country
  try:
    country = subdiv.getiterator("countryName")[0].text
  except:
    try:
      country = placename.getiterator("countryName")[0].text
    except:
      country = ''
  try:
    region = subdiv.getiterator("adminName1")[0].text
  except:
    region = ''
  try:
    name = placename.getiterator("name")[0].text
  except:
    name = ''
  
  if country:
    if region:
      if name:
        result = "%s, %s, %s" % (name, region, country)
      else:
        result = "%s, %s" % (region, country)
    else:
      if name:
        result = "%s, %s" % (name, country)
      else:
        result = country
  else:
    result = "unknown"
  
  return result

Finally, the main method of the program parses arguments, determines what coordinates to use, gets the geographical feature string, and then generates the proper HTML and copies it into the correct place.

if (__name__ == "__main__"):
  # Parse command line options
  from optparse import OptionParser
  parser = OptionParser("Usage: %prog [options] NAME FILENAME\nNote: NAME should be enclosed in quotes, unless it contains no spaces.")
  parser.add_option("--no-gps", action="store_true", dest="nogps", default=False, help="don't attempt to get coordinates from GPS")
  parser.add_option("-a", "--amail", action="store_true", dest="amail", default=False, help="generate files for Apple's Mail.app (.webarchive)")
  parser.add_option("--lat", action="store", type="float", dest="lat", help="default latitude, if GPS is unavailable")
  parser.add_option("--lon", action="store", type="float", dest="lon", help="default longitude, if GPS is unavailable")

  options, args = parser.parse_args()
  
  try: assert(len(args) == 2)
  except AssertionError:
    print "You must provide NAME FILENAME\n"
    parser.print_help()
    exit()

  try: assert(((options.lat!=None) & (options.lon!=None)) | ((options.lat==None) & (options.lon==None)))
  except AssertionError:
    print "You must provide both lat and lon arguments\n"
    parser.print_help()
    exit()

  # Modules for creating signature file
  import tempfile, os

  try:
    if (options.nogps): raise NoGPSError
    lat, lon = get_lat_lon()
    print "Got %f, %f from gps" % (lat, lon)
  except NoGPSError:
    if (not options.nogps): print "Failed to get coordinates from GPS"
    try:
      assert((options.lat!=None) & (options.lon!=None))
    except AssertionError:
      print "No default coordinates provided; exiting"
      exit()
    lat, lon = options.lat, options.lon
    print "Using %f, %f as coordinates" % (lat, lon)

  geostring = get_geo_string(lat, lon)

  lac, loc = 'N', 'E'
  if lat < 0:
    lat *= -1
    lac = 'S'
  if lon < 0:
    lon *= -1
    loc = 'W'

  file = tempfile.NamedTemporaryFile()

  file.write("\n

") file.write(args[0]) file.write("

\n

") file.write("Current meatspace coordinates, hot from the GPS receiver card in my laptop:") file.write("

\n

") file.write("%.3f %s latitude, %.3f %s longitude" % (lat, lac, lon, loc)) file.write("

\n

") file.write("Nearest geographical feature: " + geostring) file.write("

\n") file.write("\n") file.flush() if (options.amail): os.system("textutil -convert webarchive -output " + args[1] + " " + file.name) else: os.system("cp " + file.name + " " + args[1])

More Mandelbrot

I just recently revisited the M-Set code from my Perl Snippets post. The code I had was pretty ugly, so I decided to rewrite it in Python. The result is not only a lot cleaner and easier to understand, but it’s also a lot faster:

$ time python mandel.py > \dev\null
real	0m0.051s
user	0m0.036s
sys	0m0.010s
$ time perl mandel.pl > \dev\null
real	0m3.518s
user	0m3.463s
sys	0m0.029s

You can find the code here.

This script works well for zooms, as long as you stay below a few thousand iterations. The following picture was generated with x=-1.1887204, y=-0.3032472, width=0.01 and 150 iterations.