Mitch Hedberg and GIS

So, a recent post on Reddit highlighted a Mitch Hedberg joke.

“La Quinta” is Spanish for “next to Denny’s.”

Thinking about this, I realized we could use GIS to find the number of La Quintas that are next to Denny’s. Last night after the kids were asleep, I sat down with a beer (Sierra Nevada) and figured out how I could get the data and perform the calculation.

First, I visited La Quinta’s website and their interactive map of hotel locations. Using Firebug, I found “hotelMarkers.js” which contains the locations of the chain’s hotels in JSON. Using a regular expression, I converted the hotel data into CSV.

A sample of the data, before and after conversion to CSV.

Sample of the data. Click for full size.

Next, I went to Denny’s website and their map. Denny’s uses Where2GetIt to provide their restaurant locations. They provide a web service to return an XML document containing the restaurants near the user’s location (via GeoIP) or by a specified address. Their web service also has a hard limit of 1,000 results returned per call. Again, I used Firebug to get the URL of the service, then changed the URL to search for locations near Washington, DC and Salt Lake City, Utah, with the result limit set to 1000 and the radius set to 10000 (miles, presumably). From this, I was able to get an “east” and “west” set of results for the whole country. These results were in XML, so I wrote a quick Python script to convert the XML to CSV.

#!/usr/bin/env python

import sys
from xml.dom.minidom import parse
dom = parse(sys.argv[1])
tags = ("uid", "address1", "address2", "city", "state", "country", "latitude", "longitude")
print ",".join(tags)
for poi in dom.getElementsByTagName("poi"):
    values = []
    for tag in tags:
        if len( poi.getElementsByTagName(tag)[0].childNodes ) == 1:
            values.append( poi.getElementsByTagName(tag)[0].firstChild.nodeValue )
    print ",".join( map(lambda x: '"'+str(x)+'"', values) )

I then loaded the CSVs into PostgreSQL. First, I needed to remove the duplicates from my two Denny’s CSVs.

CREATE TABLE dennys2 (LIKE dennys);
DROP TABLE dennys;
ALTER TABLE dennys2 RENAME TO dennys;

Then I added a Geometry column to both tables.

SELECT AddGeometryColumn('laquinta', 'shape', 4326, 'POINT', 2);
UPDATE laquinta SET shape = ST_SetSRID(ST_MakePoint(longitude, latitude),4326);
SELECT AddGeometryColumn('dennys', 'shape', 4326, 'POINT', 2);
UPDATE dennys SET shape = ST_SetSRID(ST_MakePoint(longitude, latitude),4326);

From there, finding all of the La Quinta hotels that live up to their name quite easy.

SELECT, d.state, d.shape <-> l.shape as distance,
      ST_MakeLine(d.shape, l.shape) as shape
FROM dennys d, laquinta l 
WHERE (d.shape <-> l.shape) < 0.001 -- 'bout 100m

Here are the 29 cities that are home to La Quinta – Denny’s combos:

  • Mobile, AL
  • Phoenix, AZ
  • Bakersfield, Irvine, & Tulare, CA
  • Golden, CO
  • Orlando & Pensacola, FL
  • Augusta & Savannah, GA
  • Lenexa, KS
  • Metairie, LA
  • Amarillo, Austin, Brenham, College Station, Corpus Christi, Dallas, El Paso (2!), Galveston, Irving, Killeen, Laredo, Lubbock, McAllen, San Antonio, & The Woodlands, TX.

You might have noticed El Paso, Texas, has two La Quinta – Denny’s combos. Both happen to also be on the same street, just a few miles apart.

Gateway Boulevard to the West...

Gateway Boulevard to the West…

and further east on Gateway...

and further east…

I love that the second one even shares a post for both of their highway-scale signage.

So out of the 833 La Quintas and 1,675 Denny’s, there are 29 that are very close (if not adjacent) to one another. So, only 3.4% of the La Quintas out there live up to Mitch Hedberg’s expectations.

GIS: coming up with solutions for the problems no one asked!

Update: Chris in the comments made a good point about projection. So here’s the data reprojected into US National Atlas Equal Area and then limited to 150 meters distance between points.

SELECT, d.state, ST_Transform(d.shape,2163) <-> ST_Transform(l.shape,2163) as distance 
FROM dennys d, laquinta l 
WHERE (ST_Transform(d.shape,2163) <-> ST_Transform(l.shape,2163)) < 150

This yields 49 pairs (or 5.8% of all La Quintas):

  • Alabama: Huntsville, Huntsville, Mobile
  • Arizona: Phoenix, Tempe, Tucson
  • California: Bakersfield, Bakersfield, Irvine, South San Francisco, Tulare
  • Colorado: Golden
  • Florida: Cocoa Beach, Orlando, Pensacola, St Petersburg
  • Georgia: Augusta, Savannah
  • Illinois: Schaumburg
  • Indiana: Greenwood
  • Kansas: Lenexa
  • Louisiana: Metairie
  • New Mexico: Albuquerque
  • Oregon: Salem
  • Texas: Amarillo, Austin, Brenham, College Station, Corpus Christi, Dallas, Dallas, El Paso, El Paso, Galveston, Irving, Killeen, Laredo, Live Oak, Lubbock, McAllen, San Antonio, San Antonio, San Antonio, Victoria, The Woodlands
  • Utah: Midvale
  • Virginia: Virginia Beach
  • Washington: Auburn, Seatac

This modification to the query introduced a new oddity: Huntsville, AL has a Denny’s that is between two La Quintas, which is why it’s on the revised list twice.

Between two La Quintas

A Denny’s adjacent to one La Quinta and about 400 feet from another…

Update: I uploaded a dump of the data to Github, if you want to explore the data on your own.

Posted in Data, Google Maps, In the News, Teaching, Technology, Tools and Scripts, Web Mapping | 4 Comments

Using a TMS in JOSM’s Download Data window

Now that our crowdsourced update to New Jersey’s land use map has had some serious contributions – we’re just shy of 20,000 user-contributed points of new (2007-2012) urbanization – I realized that many of these locations are in need of updating in OpenStreetMap.

JOSM is my preferred editor for OpenStreetMap. I find it to be incredibly robust and powerful. I also love that it’s extensible; there are plenty of great plugins and it integrates with web services well. One thing that JOSM supports is WMS and TMS background imagery. I often use the freely-available 2012 imagery for New Jersey as a base for my edits. We’re also using the 2012 WMS for our NJ MAP “Growth” crowdsourcing web app to help identify areas of recent development in NJ. If you haven’t seen the app, check it out. It’s an easy to use app where 2012 imagery is presented along with a black mask derived from the 2007 urban lands in the Land Use/Land Cover data. Simply put, if you see a building on the aerials, click on it and tell us what it is.

Because those clusters of single family housing built post-2007 are not likely reflected in OpenStreetMap, I wanted to see if I could provide some base roadways for the new subdivisions as well as clean up the land use imported into OSM, which was from 2002.

I have a WMS service of the points in the Growth app, served up by GeoServer. GeoServer is also capable of providing the same data in TMS. While I was going to simply add a link to the WMS to JOSM, so that I could see the project’s contributions along side the aerial photo and the OSM data, I realized that it wouldn’t be as useful in the main map interface, because I’d only see the points after I downloaded some OSM data.

Screen Shot 2013-11-18 at 4.26.56 PM To help you find your area of interest, JOSM includes several OSM-derived map services  through the Download Data window. I did not realize that the list of layers included any TMS layers you added to JOSM. So instead of adding the new urban points from Growth as a WMS, I added them as a TMS. Now, instead of browsing for locations using Mapnik tiles, I can look for clusters of development points.

Screen Shot 2013-11-18 at 4.27.04 PMGranted, it’s not too meaningful – it’s just the points without any other background information – but it’s easy to find areas that likely need attention.

Screen Shot 2013-11-18 at 4.27.32 PMThis area was forest in 2002, but in 2012 is a new housing development. I was able to add the roads of the new development, as well as clean up the surrounding land use.

Screen Shot 2013-11-18 at 4.32.14 PMNow, I’ll still need to refer to another available source (such as NJ’s road network or our parcel data – both freely available under OSM-friendly licenses) for things like the road names, being able to focus in on areas that need updating in OSM will help us all improve the New Jersey portion of the world’s best free map.

If you want to help update OSM using NJ MAP: Growth as a guide, add the following TMS to JOSM:

Posted in Data, OpenStreetMap, Web Mapping | Tagged , , , , , , | Comments Off on Using a TMS in JOSM’s Download Data window

Portfolio Workshop

Yesterday, I held a Portfolio Workshop for our Rowan Geography and Environment students. The audio of my talk along with the slides are available on the workshop web page. While I think the talk and discussion was well-received, the audio is unedited, so you’ll hear plenty of ums and ahs before I hit my stride. I’ll probably repost this to YouTube with a clean narration at some later date.

Posted in Teaching | Tagged , , , , , | Comments Off on Portfolio Workshop

Rendering weirdness with reference data

Continued from yesterday’s post: Bizarre GeoServer WMS Rendering

Both the NJ Landscape Project habitat data and the Municipalities data are complex. It is hard to pinpoint a rendering issue when the data is very complex. I prepared a reference data set composed of a few multipart features. You can download the reference shapefile here.

Reference multipart polygons. Click for full size.

I created the polygon features in ArcGIS Desktop as a shapefile. I then imported the shapefile into my PostGIS database twice, once using PG_GEOMETRY & ArcSDE, the other using OGR. I then prepared two WMS layers in GeoServer, one based off of the SDE import and one from the OGR import. The SDE import causes GeoServer to improperly render the data. GeoServer can render the OGR import without issue.

Improperly rendered multipart features.

Improperly rendered multipart features.

I exported the geometries of each table as WKT: hexagons imported using SDE; hexagons imported using OGR. Note that the coordinates are essentially the same, but the “wrappers” around the coordinates are different. You can also tell that it is a rendering issue; in GeoServer’s layer preview mode, you can click on the OpenLayers map to identify a feature. Clicking on one of the erroneous polygons yields no results.

It appears that GeoServer is improperly handling MULTIPOLYGON features. POLYGON features (with multiple parts) are properly rendered. This is also what is being suggested in the support forum.

One other issue I encountered was that OGR would not import the shapefile without the -nlt flag set.

$ ogr2ogr -f "PostgreSQL" PG:"pg connection string" ./testing.shp -nln test_hexagons_ogr 
Warning 1: Geometry to be inserted is of type Multi Polygon, whereas the layer geometry type is Polygon.
Insertion is likely to fail
ERROR 1: INSERT command for new feature failed.
ERROR:  new row for relation "test_hexagons_ogr" violates check constraint "enforce_geotype_wkb_geometry"

Command: INSERT INTO "test_hexagons_ogr" ("wkb_geometry" , "id") VALUES ('snip'::GEOMETRY, 3) RETURNING "ogc_fid"
ERROR 1: Terminating translation prematurely after failed
translation of layer testing (use -skipfailures to skip errors)

$ ogr2ogr -f "PostgreSQL" PG:"pg connection string" ./testing.shp -nln test_hexagons_ogr -nlt POLYGON
Posted in Data, Web Mapping | Tagged , , | 5 Comments

Bizarre GeoServer WMS rendering possibly due to ArcGIS handling of PostGIS data

I just submitted a problem to the OpenGeo support forum regarding an unusual bug I’ve discovered in GeoServer. Multipart polygon geometries are producing “unintentional art” when rendered through GeoServer.

NJ Landscape Project data improperly rendered.

Several GIS data sets that have been brought into PostGIS using ArcGIS 10.0 (with the PG_GEOMETRY keyword) have then exhibited some bizarre rendering issues when rendered using GeoServer. Above, the Landscape Project data appears to have some random triangular polygons that stretch beyond New Jersey’s western border. We used ArcGIS to process this layer as I wanted a Union (ESRI not PostGIS terminology) between Landscape and NJ’s municipality layer. ArcGIS Desktop’s method of storing multipart polygons in PostGIS seems valid – the table renders just fine in ArcGIS Desktop and QGIS – but is a little different than what GeoServer’s rendering engine is expecting.

Here’s a “fixed” version that enables GeoServer to render the data correctly.

Proper way to render Landscape Project data.

The “fix” comes in the form of a simple table view. By forcing the geometry field through ST_Force_Collection and then setting the table view up as a new WMS layer, the problem polygons disappear. Frankly, I was surprised this worked, but I assumed that if ArcGIS and QGIS could handle the data, the problem wasn’t the data, but GeoServer’s ability to process it. ST_Force_Collection is just enough to massage the data into a GeoServer-friendly form.

I’m not sure if it’s something to do with ArcGIS’s PostGIS support directly, but I have already experienced issues with weirdness due to inconsistent handling of circular features between ArcGIS, PostGIS and OGR. I’m going to generate a shapefile of some bizarre multipart shapes and see how the shapes are handled via SDE & PG_GEOMETRY and well as OGR to PostGIS. I’ll share my findings here along with the source data.

Posted in Data, Technology, Visualization, Web Mapping | Tagged , , , , , | 1 Comment

Beach erosion in Ocean City

While not as hard hit as the towns near Seaside Heights, Ocean City felt the effects of Sandy and is still in the process of cleanup. Back in July, I took a Gigapan photo of the Ocean City beach from the south side of the Music Pier.

View on

Two weeks after Sandy, I took another shot from roughly the same spot.

View on

The dunes present in the top photo are completely gone. Many of the shops fared well, a little broken glass and best signs, but overall, intact.

Here’s another panorama of the boardwalk in front of the Music Pier.

View on

Posted in In the News, Technology, Visualization | Tagged , , , | Comments Off on Beach erosion in Ocean City

From NOAA imagery, extent of Seaside Heights damage can be seen

The amusement piers in Seaside Heights have been significantly damaged.

Funtown Pier partially destroyed.

Funtown Pier. Click through for full size.

A roller coaster at Casino Pier can be seen sitting in the surf. Two buildings at the end of the Pier are entirely gone.

Casino Pier. Click through for full size image.

For those that are unfamiliar with the two piers, here’s a before and after.

Funtown Pier, before and after.

Casino Pier, before and after

Posted in Aerial Photography, In the News | Tagged , , , | Comments Off on From NOAA imagery, extent of Seaside Heights damage can be seen

Post Sandy imagery available from NOAA

NOAA‘s web viewer of post-Hurricane Sandy imagery is now live. Many portions of the Jersey Shore received incredible damage. Below are some before and after pictures compiled from the viewer. I’m hoping to get NJ’s March 2012 imagery available as a standard baseline for comparison. Click for full size images.

Brick Beach 3

Mantoloking Bridge

Lyman Street, now an inlet

Atlantic City Boardwalk, by the casinos

abandoned portion of the AC Boardwark, now destroyed

One bright note: Lucy the Elephant survived!

Posted in Aerial Photography, In the News | Tagged , , , , | Comments Off on Post Sandy imagery available from NOAA

Hurricane Sandy live data resources

As Hurricane Sandy is about to make landfall somewhere in NJ or DE, keep safe and dry. Move to a safe location and stay there. Don’t attempt to travel after dusk due to reduced visibility and the intensifying storm.

Once you’re somewhere safe and have power and Internet access, you can keep apprised of the storm through the following links. I’ll keep the list below here updated throughout the storm.

Live Maps & GIS Data


Crowdsourced Reporting

Offices of Emergency Management

 Static Data and Maps

Posted in Data, Events, In the News, Technology | Comments Off on Hurricane Sandy live data resources

Pulling data from OSM using XAPI, bringing it to ArcGIS

I was recently asked by a friend if I knew of an ArcGIS-friendly data set representing the highway interchanges within New Jersey. I knew the information was already available in OpenStreetMap, but how could I quickly pull just the highway interchanges out of OSM and bring them into ArcGIS?

Using the XAPI, I was able to pull just the relevant information from OSM. In this case, I only wanted nodes that contained a highway key with the value of motorway_junction. I used MapQuest‘s XAPI service and quickly generated a URL using a nifty map mashup. I had to manipulate the XAPI URL manually to change the bounding box; New Jersey was too large to search with a statewide bbox, so I broke the state into four portions by changing the latitude ranges. Making smaller requests for the part of the state between latitude 39 and 40, then 40 and 41, allowed my request to go through without timing out.

I was then left with four XML files that contained the highway interchanges. ArcGIS doesn’t know what to do with this, so I whipped up some Python to convert all four tables into one tab-delimited file. This code (included below) can then be used to add the points to ArcMap using the Add XY Data functionality.

All said and done, this project took me about an hour. I was able to send the data off in about 30 minutes, but then spent another thirty minutes cleaning up the script, modifying it to spit out any available key, not just the “ref” tag. The script is GPL, free for you to use as you see fit. Just remember to give credit where it’s due, attributing your data back to OSM.

Taking the time to do this for a friend is a good reminder of why I push our students to do more than just the GIS curriculum. Being familiar with OSM gives you a different point of view on GIS data, schema, and community than what you would gain from most GIS courses and entry-level positions. GIS is already about problem solving; being able to do some basic programming will enable you to tackle much larger problems efficiently. And providing technical assistance to others just makes you a better person.

Posted in OpenStreetMap, Technology, Tools and Scripts | Tagged , , , , , | 1 Comment