metachronistic

A couple years ago we got iPhones, and one of my favorite apps is the RunKeeper app, which tracks your outdoor activities using the phone’s built-in GPS. When I first started using it I compared the results of the tracks from the phone to a Garmin eTrex, and they were so close that I’ve given up carrying the Garmin. The fact that the phone is always with me, makes keeping track of all my walks with Nika, and trips to work on my bicycle or skis pretty easy. Just like having a camera with you all the time means you capture a lot more images of daily life, having a GPS with you means you have the opportunity to keep much better track of where you go.

RunKeeper records locations on your phone and transfers the data to the RunKeeper web site when you get home (or during your trip if you’ve got a good enough cell signal). Once on the web site, you can look at the tracks on a Google map, and RunKeeper generates all kinds of statistics on your travels. You can also download the data as GPX files, which is what I’m working with here.

The GPX files are processed by a Python script that inserts each point into a spatially-enabled PostgreSQL database (PostGIS), and ties it to a track.

Summary views allow me to generate statistics like this, a summary of all my travels in 2010:

Another cool thing I can do is use R to generate a map showing where I’ve spent the most time. That’s what’s shown in the image on the right. If you’re familiar at all with the west side of the Goldstream Valley, you’ll be able to identify the roads, Creek, and trails I’ve been on in the last two years. The scale bar is the number of GPS coordinates fell within that grid, and you can get a sense of where I’ve travelled most. I’m just starting to learn what R can do with spatial data, so this is a pretty crude “analysis,” but here’s how I did it (in R):

library(RPostgreSQL)
library(spatstat)
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname="new_gps", host="nsyn")
points <- dbGetQuery(con,
    "SELECT type,
        ST_X(ST_Transform(the_geom, 32606)) AS x,
        ST_Y(ST_Transform(the_geom, 32606)) AS y
     FROM points
        INNER JOIN tracks USING (track_id)
        INNER JOIN types USING (type_id)
     WHERE ST_Y(the_geom) > 60 AND ST_X(the_geom) > -148;"
)
points_ppp <- ppp(points$x, points$y, c(min(points$x), max(points$x)), c(min(points$y), max(points$y)))
Lab.palette <- colorRampPalette(c("blue", "magenta", "red", "yellow", "white"), bias=2, space="Lab")
spatstat.options(npixel = c(500, 500))
map <- pixellate(points_ppp)
png("loc_map.png", width = 700, height = 600)
image(map, col = Lab.palette(256), main = "Gridded location counts")
dev.off()

Here’s a similar map showing just my walks with Nika and Piper:

And here's something similar using ggplot2:

library(ggplot2)
m <- ggplot(data = points, aes(x = x, y = y)) + stat_density2d(geom = "tile", aes(fill = ..density..), contour = FALSE)
m + scale_fill_gradient2(low = "white", mid = "blue", high = "red", midpoint = 5e-07)

I trimmed off the legend and axis labels:

I took Friday off from work again this week and brewed Overflow Old Ale. I drained the chilled wort directly onto the yeast cake from last week’s brewing effort, so I should get a nice fast start from the yeast. The recipe was supposed to be a high alcohol (7%) lightly hopped ale, but instead of an expected 70% efficiency at extracting sugars from the grain, I wound up at 57%. That’s the lowest efficiency I’ve ever gotten. I wouldn’t mind so much if they were consistent, but they’re all over the map, and I don’t know what is causing the variation. Same base grains, same procedure and equipment, same thermometer, etc. It’s a mystery.

The good news is that nothing went wrong, and beer is beer regardless of whether it turned into the batch you were planning to brew. I’m sure I’ll enjoy drinking it at least as much as I did brewing it.

The batch is named for the overflow that’s currently moving across the top of the frozen Creek. I took advantage of it and used it in my counterflow chiller to bring the boiling wort down to pitching temperatures in 10–15 minutes. The water is at least a foot deep, and I don’t think it’s been warm enough for it to be a natural phenomenon. Last year’s overflow event was rumored to be a large discharge from the mine several miles upstream from our house, and this event looks similar. Despite providing chiller water, it’s a drag because it’ll take at least a week for it to freeze solid again and until it snows it’ll be treacherous to walk on. I had hoped the overflow would have flooded the DNR pond again, as it was really fun ice skating on it last year, but the floodwaters aren’t high enough to cross onto the pond. The photo at the top is taken where the Creek passes the pond.

The photo on the right shows the Google Maps image of the location I took the photo using the iPhone web app I wrote. This morning I added some code to draw a blueish buffer around the GPS location that’s sized using the reported GPS accuracy. It also scales the map so that the entire property fits on the screen (I zoomed in before taking the photo on the right). I think the app is pretty much finished, except I want to add Fish and Game Management Units to it. Unfortunately, they don’t have GIS files for the smaller Management Areas, so it won’t be as useful as it could be. You can see I’ve drained my battery playing with it!

Now that my two year web hosting contract with BlueHost is coming up, I’ve decided to migrate over to Linode. It’s about twice as expensive, but you get your own Xen virtual machine to do with as you want. With BlueHost I wound up needing to compile a whole bunch of software packages (gnuplot, remind, spatialite, a variety of Python modules, etc.) just to get my site to work. And I couldn’t upgrade PostgreSQL, which left me unable to use PostGIS. With my Linode, I put the latest Ubuntu on it, and was able to install all the software I need.

As mentioned, PostGIS is one of the benefits of this. Instead of loading the Fairbanks North Star Borough tax parcel database into spatialite and using a CGI program to extract data, I’m able to query the same data in PostGIS directly from PHP. That means I can do all sorts of cool stuff with the data I’m extracting.

Here, I am using the Google Maps API to draw the property polygon and put a marker at the location from the iPhone’s GPS. You can see the result on the right, taken while eating lunch at the office. If you’ve got a GPS-enabled smart device and you’re in the Fairbanks North Star Borough, the link is http://swingleydev.com/gis/loc.html.

The code that drives the map is at the bottom of this post. The complicated query in the middle of it gets the outermost ring from the polygon, and dumps each vertex as a separate row. This section is in the middle of the Google Maps JavaScript that builds the overlay polygon. Note that this doesn’t handle the few properties in the database that have multiple polygons in their geometry, and it ignores the holes that might be present (if a small chunk out of the middle of a property was sold, for example). I expect these cases are pretty infrequent.

One enhancement would be to include a circular buffer around the point marker so you could get a sense of where the GPS accuracy places you on the map. This is complicated with Google because of the custom projection they’re using. I think it’s called “Google Mercator,” but spatialreference.org lists several projections that look like they might be compatible. I’ll have to give those projections a try to see if I can make a projected circular buffer that looks like a circle. It would also be nice to set the map extent to match the extent of the parcel you’re inside. I don’t know if that’s part of the Google Maps API or not.

Here’s the code:

<script type="text/javascript">
function initialize() {
  if (GBrowserIsCompatible()) {
      var map = new GMap2(document.getElementById("map_canvas"));
      var mapControl = new GMapTypeControl();
      map.addControl(mapControl);
      map.setCenter(new GLatLng(<?php echo $lat; ?>, <?php echo $lon; ?>), 14);
      var polygon = new GPolygon([
        <?php
            $query = "
SELECT
    ST_X(ST_PointN(ST_ExteriorRing(ST_GeometryN(wkb_geometry, 1)),
        generate_series(1, ST_NPoints(ST_ExteriorRing(ST_GeometryN(wkb_geometry, 1))))))
        AS x,
    ST_Y(ST_PointN(ST_ExteriorRing(ST_GeometryN(wkb_geometry, 1)),
        generate_series(1, ST_NPoints(ST_ExteriorRing(ST_GeometryN(wkb_geometry, 1))))))
        AS y
FROM tax_parcels
WHERE ST_Intersects(wkb_geometry, ST_SetSRID(ST_MakePoint($lon, $lat), 4326));";
            $result = pg_query($query);
            $num_rows = pg_num_rows($result);
            for ($i = 0; $i < $num_rows - 1; $i++) {
                $row = pg_fetch_row($result);
                list($vlon, $vlat) = $row;
                echo "new GLatLng($vlat, $vlon),";
            }
            $row = pg_fetch_row($result);
            list($vlon, $vlat) = $row;
            echo "new GLatLng($vlat, $vlon)";
            pg_free_result($result);
        ?>
      ], "#f33f00", 2, 1, "#ff0000", 0.2);
    map.addOverlay(polygon);
    var point = new GLatLng(<?php echo $lat; ?>, <?php echo $lon; ?>);
    map.addOverlay(new GMarker(point));
  }
}
</script>

Today I went for a walk on the Creek with Nika and Piper and it occurred to me that the iPhone ought to be able to tell me not only where I am (which it does quite well with the Google-driven map application), but on whose property I’m walking through. It took me a little time, but I now have a webapp (http://swingleydev.com/gis/loc.html) that can do this (don’t bother clicking on the link unless you’re on an iPhone or other GPS-enabled device).

The result is a pair of pages. The first one shows you your current location, speed, and heading. With a press of a button (and a delay while the database finds the property owner) you can see the Borough information on the property you’re currently inside (assuming you’re in the Fairbanks North Star Borough—this page doesn’t do non-Borough residents any good).

Here’s how I set it up:

My web hosting provider doesn’t have a new enough version of PostgreSQL to run PostGIS, which means I need to use the spatially-enabled version of SQLite3, called Spatialite. To get the parcel database from the shapefile to Spatialite requires the following steps. Some parcels are eliminated in this process, but it’s a small fraction of the parcels in the database:

Download the parcels shapefile:

$ wget ftp://co.fairbanks.ak.us/GIS/tax_parcels.zip

Unzip it and insert it into a spatially enabled PostgreSQL database using ogr2ogr. Crash the process immediately:

$ unzip tax_parcels.zip
$ ogr2ogr -f "PostgreSQL"
    -t_srs EPSG:4326
    -overwrite -skipfailures
    PG:"dbname='test'"
    tax_parcels.shp

Fix the column definitions:

sql> DELETE FROM tax_parcels;
sql> ALTER TABLE tax_parcels ALTER COLUMN
    sqft_calc TYPE numeric (22,12);

Re-import the data:

$ ogr2ogr -f "PostgreSQL"
    -t_srs EPSG:4326
    -append -skipfailures
    PG:"dbname='test'"
    tax_parcels.shp

Convert the geometry column to MULTIPOLYGON:

sql> ALTER TABLE tax_parcels DROP CONSTRAINT "enforce_geotype_wkb_geometry";
sql> ALTER TABLE tax_parcels ADD CONSTRAINT "enforce_geotype_wkb_geometry"
     CHECK (geometrytype(wkb_geometry) = 'MULTIPOLYGON'::text OR
            wkb_geometry IS NULL OR
            geometrytype(wkb_geometry) = 'POLYGON'::text);
sql> UPDATE tax_parcels SET wkb_geometry = ST_Multi(wkb_geometry);
sql> ALTER TABLE tax_parcels DROP CONSTRAINT "enforce_geotype_wkb_geometry";
sql> ALTER TABLE tax_parcels ADD CONSTRAINT "enforce_geotype_wkb_geometry"
     CHECK (geometrytype(wkb_geometry) = 'MULTIPOLYGON'::text OR
            wkb_geometry IS NULL);
sql> UPDATE geometry_columns SET type='MULTIPOLYGON'
     WHERE f_table_name='tax_parcels' AND f_geometry_column='wkb_geometry';

Re-import the data (there will be thousands of errors, but this insert should add the MULTIPOLYGON rows that weren’t inserted the first time around):

$ ogr2ogr -f "PostgreSQL"
    -t_srs EPSG:4326
    -append -skipfailures
    PG:"dbname='test'"
    tax_parcels.shp

Get rid of the illegal polygons:

sql> DELETE FROM tax_parcels WHERE NOT ST_IsValid(wkb_geometry);

Convert to spatialite:

$ ogr2ogr -f "SQLite"
    tax_parcels.sqlite
    PG:"dbname='test'"
    -dsco SPATIALITE=YES

With the data in the proper format, I built a mobile web page that can pull location information from the mobile device, display it on the screen, and pass this information to a page that finds and displays the parcel information from the Spatialite database. JavaScript handles gathering the location information, but because I don’t control the web server, I had to write a CGI script to pass a query to spatialite and format the results as a list. This CGI is called by the main page from a <form> element. The spatialite query looks like this:

sqlite> SELECT street_add, owner_firs, owner_last,
           owner1, owner2, owner3, mail_add, ci_st_zip,
           round(acres_calc, 1) as acres, ’$’ || total as total,
           ’$’ || land as land, ’$’ || improvemen as improvements,
           pan, sub, block, lot, road_water, lot_size, units,
           neighborho, primary_us, tax_status, tax_year, mill_rate,
           business, year_built, situs_numb, situs_name
        FROM tax_parcels
        WHERE Intersects(GEOMETRY, SetSRID(MakePoint(lon, lat), 4326));

Again, implementing this would be a lot easier if I could install PostGIS on the server and use PHP to access the data directly. Because I can’t, spatialite and CGI do the trick.

Update: I added a few more steps to convert the initially imported POLYGON layers to MULTIPOLYGON, which then allows us to include the MULTIPOLYGON rows from the shapefile.

We got 2 inches of snow yesterday (October 26th), so the wait is finally over.

I made the mistake of riding my bicycle to work yesterday, as the snow was falling. It wasn’t too bad on my way in to work, but by the time I left, more than an inch of snow had fallen and the roads hadn’t been plowed. I do have studded, knobby tires on my bicycle, but they’re don’t work very well in situations where the snow is deeper than the tread. I managed to stay upright the whole way home, but it was some white-knuckle, one-wheel drive bicycling.

Note: Yesterday’s first real snowfall was the 8th latest in the 62 year historical record I have access to for the Fairbanks airport station. I’m not sure where the statistics reported in Tuesday’s newspaper came from.

Another Toy Camera photo. It’s a shot down our driveway from the red cabin.

Having a camera attached to the iPhone really is an amazing thing. Inasmuch as I carry the thing with me most places, that means I have a camera most of the time. The first shot is of the sunrise I saw this morning riding my bicycle in to work (using the Pano app to make a panorama).

This is a shot of the coffee table we’ve got under the west window. It’s an underutilized spot, but the presence of the window and the Monitor heater next to it means there’s not a lot we can do with it. I used the Toy Camera app to take this shot, which takes a photo and then applies some sort of filter to the image. I’m not sure what filter was used here, but I like the way it looks. Vintagey.

Another iPhone panorama, this time with the Pano app. I had been using AutoStitch, which makes panoramas from a series of existing photos. Pano takes a different approach: you shoot a series of photos, in order from left to right, from inside the app. After each image, the app asks if you’re happy with the shot, or if you’d like to retake it. If you’re happy with it, it stores it, and then shows you a semi-transparent slice of the right edge of the photo superimposed over the left side of the current camera viewscreen. This makes it fairly easy to line up each shot as you pan across your scene. When you’ve taken all the photos you want, the app joins the images together and saves it to the Camera Roll on the iPhone.

The upside to Pano is that it’s much easier to get well aligned images, as long as there’s enough contrast in the individual pictures to allow you to line them up as you’re shooting. The down side is that the only layout the app can handle is a single row of landscape or portrait shots. AutoStitch can join photos in any combination. The panorama at the bottom of my previous post (our back yard) was built from two rows of four photos (8 images total). The top row included a nicely exposed blue sky, and the bottom row was primarily the tussock–permafrost landscape of our backyard. Even though there are some obvious artifacts in the final image, it would be hard to get such a nice overall exposure with Pano and the iPhone camera.

We’ve had quite a few moose encounters over the last few days. This morning I felt what I thought was an earthquake, but when I looked out the window to see what was happening, there was a bull moose next to the deck trying to work the velvet off his antlers by rubbing them on the deck stair railing. He didn’t completely tear off the railing, but he came close. Along the way he stabbed our water tank several times (you can see the damage to the foam insulation in the photo on the right), knocked over a wheelbarrow, and ripped my National Weather Service rain gage off the dog yard fence. We’d seen the same moose yesterday afternoon along the power line using a power pole support wire on his antlers.

I shot a couple videos of this morning’s action. The first one is after he started ramming the water tank. I decided I needed to go outside and scare him away before he did any more damage to the house. You can see that he thinks twice about charging me before he runs off. I thought for sure he was going to barrel right into the Jeep. Here’s the video:

The second video shows him going to town on the stair railing. This was was shot through a window, so it’s not quite as sharp:

It’s only a couple days from moose season. We live in a bow-only hunting area, but the regulations here allow bow hunters to take any male with more than a “spike/fork” set of antlers (meaning a little antler spike on one side of his head and a forked antler on the other side). This moose is much larger than that. I’m not sure if he’d be legal in a typical hunting area where hunters can shoot any moose with more than a 50” antler spread, but he’s close.

Hunting season means cool temperatures, lots of sun, and gorgeous colors. Here’s an iPhone panorama shot of our back yard: