Now that I have sample tiles from my SRTM data, what format do I want to store them in?

An 8 bit paletted image is sufficient. I really only need 2 colours: transparent (below 8000′) and white (above 8000′). All the other colouring is just so it looks nice. (May save space using 64 or 128 colours). Even GIF would work, but I think PNG-8 (color type 3) is a better choice. Interestingly, even palettized PNGs support a full alpha channel for all 256 colours, see the tRNS chunk.

Do all the tiles need to have the same palette? They can, since it’s the same colormap for all the images. Not sure how to convince GDAL / other tools to reuse the palette, though, particularly with the transparency data. I’m guessing I’m going to need a post-processing step after GDAL generates the final slippy map tiles. Maybe that final PNG baking step is the time to enforce a single palette. The other processing I want to do is optimize via one of pngslim pngnq pngcrush optipng advpng. A quick test shows optipng saves 50%, I think by translating 24+8 bit images to 8 bit paletted images. I swear I saw a map hacker’s notes on comparing these tools once.

One related idea I’ve had for awhile: it’d be awesome to make a set of SRTM tiles that were a simple greyscale, 0-255, and then remap them to full colour + transparency in Javascript, in the client. That way the user could quickly choose what altitudes to highlight. I think it’s barely possible using the canvas API, not sure if it’s a good idea.

Working on a static render of my SRTM tiles. Decided to use gdaldem to do the color mapping, it does all I need.

First I needed gdal 1.7 to use gdaldem, but Debian still packages older versions. There’s a lot of instructions for compiling. Long story short, ./configure; make works fine. It’s 1 million lines of code, so I had a nice cup of coffee.

Here’s how to use gdaldem to turn an SRTM .hgt file into a PNG: gdaldem color-relief -alpha -of png N39W120.hgt nelson.cmap gdaldem.png. There’s a variety of color map options, I hand-coded my own:

30000   180     0       0       255
4572    180     0       0       255
3049    130     80      80      255
3048    80      80      80      255
2438    0       0       0       255
2437    0       0       0       0
0       0       0       0       0
nv      0       255     0       255

That maps everything from 0-2437 meters (8000′) to transparent, then a grey ramp up to 10000′, then a red ramp up to 15000′, then solid red. Missing values (NV) are coloured green; need to deal with them later. The result looks like this

Now a bit of shell scripting to run the whole set:

for f in *.hgt; do
  ../experiment/gdaldem color-relief -alpha -of png $f ../png/nelson.cmap ../png/${f%.hgt}.png;
done

It’s awfully fast, I’m impressed with GDAL speed. The output PNG isn’t so bad, either, the empty tiles end up being 2k or so each.

Next up: gdal2tiles.py to turn these into tiles I can serve for a slippy map. I might need to gdalbuildvrt first to combine all the files into one big virtual map.

Make a copy of the SRTM data set, but at 1/16th resolution

for f in *.hgt; do
gdal_translate -outsize 25% 25% $f small/$f;
done

I took a quick trip through this gdal tutorial on raster data, albeit with SRTM data instead of the example data. Things I learned:

• gdalinfo has useful flags for finding min/max values, stats, etc.
• gdal_translate is how you render a dataset to an image
  gdal_translate -ot byte -scale 1152 3184  -of PNG  N39W120.hgt test.png
  Input range doesn’t seem to autoscale, despite what the docs say. Haven’t figured out how to apply a custom colormap yet. (The FAQ is discouraging, but gdaldem color-relief might do what I need.) The PNG output from GDAL doesn’t seem to have any geodata stored in PNG chunks.
• gdalwarp is how you reproject data
• gdal_merge is how you mosaic images. You can use gdalwarp to do it, too, but it’s overkill.
  gdalwarp N39W120.hgt N39W121.hgt N38W120.hgt N38W121.hgt big.tiff
gdal_translate -outsize 50% 50% -scale 1000 3000 -of PNG big.tiff test.png
• gdal_rasterize lets you render vector data (via OGR) into raster data
• gdaladdo creates downsampled overviews for data, for faster access. In my previous experiments the overviews didn’t make mapserver any faster on my SRTM data.
• gdal has virtual formats which are a text file that describes multiple raster files as a single dataset. mapserver uses something different to accomplish this, a shapefile.

I do all my Javascript development in Chrome. But I never really sat down to learn about the Chrome developer tools other than learning how to open the console. Here’s some bits I learned reviewing the docs.

• Ctrl-shift-C lets you inspect DOM elements with mouseover. This is the same as clicking the little magnifying glass in the tools.
• In Elements, the panels on the right show the live state of the DOM, CSS, etc. Most bits are editable as well as inspectable.
• In Scripts, you can set breakpoints by clicking on the line numbers.
• The docs describe a cool-sounding network page that I don’t have in my Chrome/stable. May be only available in developer for now?
• There’s ridiculously powerful performance tools built in. Learning it all would take hours, but quickly: go to Timeline, click the round record button, reload your page, then click record again to turn off recording. There’s a profiler too, for statistical profiles. Both seem to work better in Incognito mode so you don’t have any extensions.
• The console scripting prompt has some limited support for jQuery-like DOM identification of elements.

No hacking today, went flying instead. Ended up being a passenger the whole 3 hour flight, crappy IFR weather. But as always I captured the flight on my little AMOD GPS tracker and mapped the track through my custom GPS track mapper project. The result is quite surprising; not only is data missing, but the track is wrong! In particular when we left Modesto we did not fly south, we flew more or less west, and the GPS recorded the wrong position. What happened?

My GPS tracker logs NMEA sentences once a second and records all sorts of neat data. Unfortunately I’m not graphing it yet, so I turned instead to GPS visualizer to visualize some auxilliary data: specifically the # of satellites we had for the fix as well as the horizontal dilution of precision, a measure of fix accuracy.

As you can see in the graph, the HDOP spiked really high for about 10 minutes around 23:20, or exactly when my position data was bad and then disappeared. I’d naively assumed the GPS tracker wouldn’t even record data if it didn’t have a good fix. Apparently it does. The Wikipedia article suggests I should be ignoring fixes with HDOP > 10 or so. My tracker also records vertical DOP; in general GPS isn’t as good at vertical position as horizontal, so that data may be even worse.

One other graph: # of satellites as a function of HDOP. My GPS apparently needs to see 5 satellites to give an accurate position. No idea if that’s true for all GPS chipsets (I think 4 is minimal), but this tracker uses the decent SiFR III chipset so I imagine it’s doing about as good as it can.

I’ve got a long list of things to do with my GPS track mapping program. I hope to make it a web site for others to use some day.

I wanted to re-use my useless CapsLock key to search Google. I started by modifying WinUrl, an old program I’d used forever. Seemed easy enough, just RegisterHotKey(..., VK_CAPITAL). It even seemed to somehow suppress the capital locking function, all done! Except not, some confusing thing happened where caps would be on in some windows, off in others.

I went down a rabbit hole trying to disable caps locking. It’s buried deep in the drivers somewhere in Windows. The recommended work around is to use SendInput() to fake a new CapsLock keystroke to toggle it back, or maybe SetKeyboardState() to obliterate the CapsLock bit. But about the time I was putting 2 second delays in my code to prevent recursion from triggering my own functions, I gave up on that approach.

I ended up using AutoHotKey instead with this script. I don’t like AutoHotKey, it’s way overengineered. Custom scripting language, custom compiler to native code?! But it’s also very capable. And it knows how to override CapsLock.

Awhile back I made a slippy map of parts of the US over 8000′. I did it using SRTM elevation data turned into images with mapserver, a dynamic program. Easy to get going but awfully slow to serve. I’d like to make the same map but statically. Inspired by Matt’s success with gdal2tiles.py I thought I’d try again.

My source data is the .hgt files for North America. They’re bundled together into a single index.shp file using gdaltindex, as described in the mapserver docs. A quick and dirty render is accomplished with

gdal2tiles.py N37W119.hgt ~/public_html/tmp/gdal2tiles

That does all the expected tiling and rendering at multiple scales, but doesn’t have a meaningful colourmap. But in my work with Mapserver I created a meaningful colormap:

  LAYER
    NAME         srtm
    TILEINDEX    "index.shp"
    TYPE         RASTER
    STATUS       DEFAULT
    PROCESSING   "SCALE=2438,4572"    # 914m = 3000', 2438m = 8000', 4572m = 15000'
    PROCESSING   "BANDS=1,1,1,1"
    PROCESSING   "LUT_1=0:0,1:0,73:80, 74:130, 255:180"
    PROCESSING   "LUT_2=0:0,1:0,73:80, 74:0, 255:0"
    PROCESSING   "LUT_3=0:0,1:0,73:80, 74:0, 255:0"
    PROCESSING   "LUT_4=0:0,1:255,255:255"

What that says, in mapserver-eze, is:

1. For each input point from SRTM, map the range 2438:4572 to 0:255.
   This maps everything < 8000′ to 0, everything over 15000′ to 255, and a linear range between.
2. Map the input band “1” (the altitude) to the output channels RGBA. Ie, a full 32 bit image with 8 bits of alpha.
3. LUT_4 maps everything that is 0 (< 8000′) to transparent, everything > 0 to fully visible (alpha = 255)
4. LUT_1-3 apply a colourmap. 1 (8000′) is mapped to 0,0,0 (black), 73 (10000′?) maps to 80,80,80 (grey), and 255 (15000′) maps to red. This gives a color ramp: black to grey to red.

Now I need to figure out how to generate the same colourmap using solely GDAL. And then convince gdal2tiles to do that colourmapping, probably by first turning my HGT files into GeoTIFF and then tiling the GeoTIFF.

 

Now that I have GeoDjango up and running I was able to import a bunch of airport data using Adam’s code. Back to PostGIS to play!

select facility_name,
  floor(ST_Distance_Sphere(point,
    ST_GeomFromText('POINT(-122.4522 37.7555)')))
    as distance
  from airports_airport
  where control_tower = true
  order by distance
  limit 10;

ERROR:  geometry_distance_spheroid: Operation on two GEOMETRIES with different SRIDs

Frowny. Face. What’s that mean?

Long story short, an SRID identifies which Spatial Reference System the coordinates are with reference to. It’s a bit like units: it’s not enough to say something “weighs 100”, you have to specifiy if you mean 100 lbs, or 100 kgs, or whatever. For more theoretical background, see Geodesy. There’s a giant database of SRIDs at http://spatialreference.org/. (Related keyword: EPSG). Greatest hits:

• 4326: WGS 84, the reference most GPSes use. I think many people assume this is the SRID by default.
• 4269: NAD83, a popular North American standard
• 54004: World Mercator
• 900913: Google Maps’ modified Mercator. (It spells “google” in l33t).

So the problem is I’m trying to do arithmetic with points in two different SRIDs.

select ST_SRID(point) from airports_airport limit 1;
4326

select ST_SRID(ST_GeomFromText('POINT(-122.4522 37.7555)'));
-1

select ST_SRID(ST_GeomFromEWKT('SRID=4326; POINT(-122.4522 37.7555)'));
4326

GeoDjango defaults to SRID 4326; you can use others, details here. PostGIS ST_GeomFromText defaults to an SRID of -1, which I guess means unknown. But if you’re explicit then it works fine. You can also use PostGIS to transform from one SRID to another, I believe underneath it uses PROJ.4

SELECT ST_AsEWKT(
  ST_Transform(
    ST_GeomFromText('POINT(743238 2967416)', 2249),
  4326
));
SRID=4326;POINT(-71.1776848522251 42.3902896512902)

Personally I find this all frustrating; can’t we just use WGS 84 everywhere and be done with it? Sure, you still need to project when rendering a flat map, but consider that a separate thing entirely.

Fun fact: the reference points for WGS 84 move. And are remeasured roughly yearly. Apparently Greenwich is not at 0 longitude in WGS 84.

I’ve been trying to develop code on my Windows box. It’s the machine I sit in front of, it’s fast, and I like WingIDE. I have a Linux box for “real work”, but think there’s some discipline in moving back and forth between the two systems.

For GeoDjango hacking I decided to run the PostGIS database on my Linux box and the Django webapp on my Windows box. I’d hoped this would avoid having to install a bunch of funky binaries on Windows. No such luck. The GeoDjango docs are a bit out of date and the official GeoDjango Windows installer binary is for an old version.  Here’s what I’ve learned doing it myself.

The first step is getting regular Django working. Here’s what I did:

• Installed PostGIS + dependencies on Debian (testing) via
  aptitude install postgresql-8.4-postgis
• Installed stock Python2.5 on Windows from Python.org
• Installed Django 1.2 on Windows. I forget how I did this, I think using easy_install, which in turn requires installation as well. It wasn’t hard, just followed the tutorial. Django 1.2 includes GeoDjango Python libraries.
• Installed psychopg2, the database driver, on Windows. There’s binaries at http://stickpeople.com/projects/python/win-psycopg/

GeoDjango is harder. It relies on a lot of C libraries, some of which don’t support Windows well. The main hurdle is a working GEOS DLL. The GEOS page itself doesn’t have binaries. I found binaries from three places: OSGeo4W,  the GeoDjango Windows Installer, and PostGIS. Here’s a quickie test program I wrote to try them out:

geosPath = r'.\libgeos_c-1.dll'

from django.conf import settings
settings.configure(GEOS_LIBRARY_PATH=geosPath)
from django.contrib.gis.geos import tests
tests.run()

Unfortunately, both the OSGeo4W and GeoDjango versions of the DLL don’t pass the tests. Version skew? Bug? The underlying error seems to be in memory management, something like

ERROR: test02_wktwriter (django.contrib.gis.geos.tests.test_io.GEOSIOTest)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "c:\Python25\lib\site-packages\django\contrib\gis\geos\tests\test_io.py", line 30, in test02_wktwriter
    self.assertEqual(ref_wkt, wkt_w.write(ref))
  File "c:\Python25\Lib\site-packages\django\contrib\gis\geos\prototypes\io.py", line 147, in write
    return wkt_writer_write(self.ptr, geom.ptr)
  File "c:\Python25\Lib\site-packages\django\contrib\gis\geos\prototypes\threadsafe.py", line 55, in __call__
    return self.cfunc(*args)
  File "c:\Python25\Lib\site-packages\django\contrib\gis\geos\prototypes\errcheck.py", line 87, in check_string
    free(result)
WindowsError: exception: access violation reading 0x02F12494

The library that did work for me was the version that came with PostGIS. I didn’t even install PostGIS, just grabbed the DLLs out of their zipfile. Unfortunately GEOS is spread into two files: libgeos_c-1.dll and libgeos-3-2-2.dll. The latter needs to be in a directory where Windows can find it as a normal DLL, I put it in the working directory of my own code.

GEOS was enough for me to get the tests working and my simple GeoDjango program working correctly. There’s more work to be done: I believe you want both GDAL and PROJ.4 available to GeoDjango on the Windows side. Or does PostGIS do the work on the database server? I haven’t tested that yet.

In addition to basic GeoDjango on Windows setup, here’s notes specific to my code:

• Followed instructions on getting WingIDE to work with Django
• Sweated a lot over environment considerations, specifically Python path. I hate setting up new development environments, and Django with its mix of servers and applications and developer code and third party libraries is complicated.
• I’m trying to run a script directly, not via manage.py. To make that work you need to set the DJANGO_SETTINGS_MODULE system environment variable to point to the Python path of your settings.py
• settings.py needs to specify a special database engine:
  ‘ENGINE’: ‘django.contrib.gis.db.backends.postgis’

There’s some more useful notes about Windows and GeoDjango in this blog post by Reinout van Rees.