Moving
This blog has moved to indicatrix.org. For those with feed readers, the feed is now here.
I presented a poster (above) at the California Geographical Society, hosted by UC Davis, this weekend demonstrating some of the analytics possible with OpenTripPlanner Analyst. The full poster and maps are available here.
Elevators in OpenTripPlanner
I’ve been working on the OpenTripPlanner project quite a bit lately. One thing I did a month or so ago was to implement elevator support in the routing engine. I decomposed OSM nodes tagged highway=elevator by their constituent levels, and built edges between them to represent boarding, riding on, and alighting from an elevator. I was very impressed with the friendliness and responsiveness of the community.
One challenge was parsing OSM levels. They can come from multiple sources—level_map relations, level tags or layer tags. I wanted to support all of these, or any combination (on a single elevator). I originally did this by noting the source and adding 0, 1000 or 2000 to the level, but Andrew Byrd has made an OSMLevel class which handles this much more neatly. Level maps allow levels to be named, which is quite nice: “take elevator to garage” instead of “take elevator to -1.” So, for all OSM mappers out there, here is a quick guide to making routable elevators:
- Tag your elevator nodes highway=elevator
- Add access restrictions—wheelchair=yes and bicycle=yes tags
- Use an OSM Level Map relation if there is any possibility that level names would not be clear.
EDIT: To clarify, the level annotations are always on the ways connecting to the elevators, not the elevator nodes themselves.
BART shocked Bay Area transit enthusiasts this morning with an announcement that it plans to build a new infill station, to serve all San Francisco lines, at Treasure Island.
“Treasure Island has long been underserved by transit,” said a BART employee, who is heading the project and acting as a liaison between the agencies involved. “The residents have long been frustrated that they live within a mile of rail rapid transit but cannot access it.”
The announcement has met with mixed results. While most residents of Treasure Island are ecstatic, East Bay commuters are not so pleased. One person, who lives in Castro Valley but commutes via BART daily to his job as a web designer in San Francisco’s Mission, complained that it will slow service. “BART travels at up to 80 mph in the Transbay Tube, crossing the Bay in just a few minutes. By adding a stop in the middle of that, you not only add the 20-30 seconds of time at the stop, but also the time it takes to accelerate and decelerate to that stop.”
Cyclists, however, are happy about the project. The new eastern span of the Bay Bridge will sport a bike path, but the western span, not scheduled for replacement, does not. By allowing a quick connection between the Yerba Buena Island and San Francisco, a cyclist can ride to Yerba Buena, then catch BART for a quick ride into the city. Responding to the head of the Greater Golden Gate, Gough and Geary Boulevard Cyclists’ Association, a BART employee confirmed that the cycle link is a critical part of the project, and that the station will have new electronic bike lockers, as well as easy bicycle access from the bridge.
For anyone who hasn’t caught it yet, note date of issue.
And if you still haven’t caught it, read this article.
Happy April Fools’ Day!
Conditionals in the QGIS raster calculator
I needed to do some conditionals in the QGIS raster calculator, but it doesn’t support that—or at least doesn’t seem to. But it does support logical operators, with a result of either 0 or 1. For instance, here’s the script I wrote:
# Subtract them ((DavisQuad2012-02-25T16_00_00Z@1 - DavisQuad2012-02-29T16_00_00Z@1)* # Multiply by 1 if neither is 255 (NoData), 0 otherwise (DavisQuad2012-02-25T16_00_00Z@1 != 255 AND DavisQuad2012-02-29T16_00_00Z@1 != 255)) # Subtract 32768 if either one was NoData, giving us -32768 for NoData. - (32768*(DavisQuad2012-02-25T16_00_00Z@1 = 255 OR DavisQuad2012-02-29T16_00_00Z@1 = 255))
Of course, you can’t actually put the comments in. But what it does is this: First, I subtract one raster from the other and multiply that by the logical operation that neither one contains NoData. That gives me the difference of the rasters, or 0 if either one contains NoData. Then I subtract 32768 multiplied by the inverse of the aforementioned logical operation, so any pixel with a NoData value in either of the original rasters is -32768 in the new one.
Note: I expanded on an idea from the GDAL gdal_calc.py script.
Conditional Labels in QGIS
I fairly commonly find myself in a situation where I would like to display one label for certain features and another for other features in the same layer. QGIS doesn’t have an official way to split labels up into categories, and until now I’d resorted to having two layers to render otherwise identical features. But, in the new, excellent expression based labeling from Nathan Woodrow, I realized one can use an SQL CASE statement. For instance, one time I need to use two types labels is when labeling roads in OpenStreetMap: I want to use the name tag, unless the feature has a ref tag defined (a name tag might be ‘Capital City Freeway’, while the ref tag would be ‘US 50’). Here’s how I solved that particular problem:
CASE
WHEN ref IS NOT NULL THEN ref
ELSE name
END
For now, anyhow, you’ll need to be running the dev build of QGIS. Happy GISing!
More Basemaps in QGIS
One of the more popular posts on this blog has been my piece on adding basemaps to QGIS. While the OpenLayer plugin is great, one of the things that I find dissatisfying is that it requires reprojecting your data to match the EPSG:3857 basemap. I often work in State Plane, and I’d just as soon have my data stay in that projection, which will also minimize local distortion. Well, as it turns out, one can add tiled map services as GDAL raster layers, with all the benefits that entails (e.g. reprojection). What you need to do is create an XML file like the following (which is lifted almost verbatim from the GDAL website, specifically this file):
<GDAL_WMS>
<Service name="TMS">
<!-- note: if you use this file verbatim, you *must* credit MapQuest and OpenStreetMap! -->
<ServerUrl>http://otile1.mqcdn.com/tiles/1.0.0/osm/${z}/${x}/${y}.png</ServerUrl>
</Service>
<DataWindow>
<UpperLeftX>-20037508.34</UpperLeftX>
<UpperLeftY>20037508.34</UpperLeftY>
<LowerRightX>20037508.34</LowerRightX>
<LowerRightY>-20037508.34</LowerRightY>
<TileLevel>18</TileLevel>
<TileCountX>1</TileCountX>
<TileCountY>1</TileCountY>
<YOrigin>top</YOrigin>
</DataWindow>
<Projection>EPSG:3857</Projection>
<BlockSizeX>256</BlockSizeX>
<BlockSizeY>256</BlockSizeY>
<BandsCount>3</BandsCount>
<Cache/>
</GDAL_WMS>
Change the ServerUrl to your Tiled Map Service server (this one is for MapQuest Open Tiles), then go into QGIS, Layer->Add Raster Layer and select the XML file.
A few caveats:
- The OpenLayers plugin automatically adds the required attributions, at least for OSM. This can be nice or not—nice in that you don’t have to remember to add the attribution, not so nice in that you can’t choose where to place the attribution; it’s always in the lower-right.
- Reprojecting layers with text may yield strange skewing and distortion.
- You need to be sure that you use these services legally (as you did with OpenLayers plugin).
Have fun!
You can also convert tiles to GeoTIFF using gdal_translate, but I suspect most TMS providers would prefer you didn’t.
Transit to Everywhere
Data courtesy MapQuest and OpenStreetMap CC-BY-SA, the City and County of San Francisco, and Bay Area Rapid Transit
This is an overlay of the transit and walking trip plans generated by OpenTripPlanner from Powell and Market to every other intersection in San Francisco, after Eric Fischer’s map of walking routes to every intersection in San Francisco. It brings out the transit routes but also shows well-used walking routes. The lines do not vary in width (don’t let Market Street fool you, it’s actually several lines—BART, MUNI rail in 2 directions, Muni bus, walking—very near each other). The lines fade where there are fewer routes using them, because they are rendered as black set at 10% opacity. Where there are more lines overlapping, the lines become darker, in what I believe is a log (or log-like) scale. It ended up just mostly being a map of San Francisco, with transit routes emphasized. It doesn’t show potential utilization of the transit system, because the routes are not weighted (it would probably be wise to weight the routes by the density of the block they terminate in and by their service area; i.e., estimate the number of people within the Thiessen polygon of each intersection and weight the route by that). Also, I had difficulty finding an opacity level where the usage of transit routes fades towards the end (as it clearly should) but still shows the streets that walked down by just one or two trip plans.
I think the data I used to make this map could possibly be better utilized to make a cartogram of San Francisco transit times (like another of Eric Fischer’s maps, but including transfers and walking times).
I’d also like to make a companion map using the OTP bike router. I think it could look really interesting in San Francisco, because the router will try to avoid hills.
Methodology
I set up an instance of OpenTripPlanner using a graph built from OpenStreetMap data for the San Francisco area, as well as GTFS data from BART and San Francisco Muni. I used the pre-built binaries of OTP. I then used a Python script to request directions from Market and Powell to every other intersection in San Francisco, as defined in the StIntersections dataset from here. I stored the directions in a PostGIS database. I used one machine as the OTP server, and ran the script and PostGIS on another machine, but I see no reason why they couldn’t be on the same machine. I used QGIS to render the map. For what it’s worth, I’ve open-sourced the script I wrote. It may provide a good example of how to use the OTP JSON API in Python.
GitHub Image Diffs
As you may have gathered, I like Git and GitHub. Today, I ran across a GitHub feature this is really cool and above and beyond the call of duty: not only do they produce and display diffs on text files, but also on image files! You can see an example in one of my repositories. Added points if you can figure out where the map tile is from!
Shapefiles in OpenLayers
Update 2011-12-14: It seems that a lot of people are coming here from web searches with phrases like “shapefile openlayers.” If all you want to do is display your data in OpenLayers, I’d highly recommend using a program like Quantum GIS to convert your Shapefile to a more web-friendly format like KML or GeoJSON. Both of these formats can be read by OpenLayers directly, and you’ll see faster performance and more browser compatibility than if you were to load your Shapefiles directly.
Over the last few days, I’ve been using Tom Carden’s shapefile-js library that reads ESRI Shapefiles in JavaScript, which I found via a post on the Prodevelop blog. The library is quite incredible, but his samples use a simple canvas for display. I thought it would be really cool if this could be integrated with OpenLayers, so I created a bit of JavaScript to do so. You can give it a test drive here, and look at the code on GitHub.
Basically, the library does all the heavy lifting. My code converts the shapefile shapes to WKT, which is passed to OpenLayers. Ultimately, I’d like to see an OpenLayers plugin so that you can use Shapefiles directly (i.e., an OpenLayers.Format.Shapefile). The main issue I see is that there needs to be a new strategy as well as a new format, because a) Shapefiles are made up of multiple pieces and b) we need to use the BinaryAjax loader since Shapefiles are binary.
My code seems to work well with points, lines and polygons, including the donut polygon case (to see for yourself, look at South Africa). (I did not test the donut polygon case, but I think it should work). More eyes are of course welcome! Also, the shapefile-js library can only handle pretty small Shapefiles. If I integrate this into OpenLayers, I think, long term, using a Web Worker thread to parse the Shapefile would be wise (which is another challenge to direct OpenLayers integration).
EDIT 2011-12-13 22:36 -0800: I tested the donut polygon case.
Indicatrix 