• Tag Archives QGIS
  • Experiments in Routing, Part 1

    This is the introductory post for a hopefully four-part series about using QGIS to find the shortest path between two points, not shortest as the crow flies, but following a given network of roads. This is called routing, it’s what’s Google Maps and other mapping software uses, and it relies on graph theory and network analysis to do its job. I’ll talk about the what and the why of this little experiment here; the how (for three different versions of how) will be the subject of subsequent posts.

    UPDATE: Part 2 can be found here.

    The reason I’m looking at all this goes back to my interest in cycling tourism, and my attempts to identify cycling-accessible amenities — convenience stores, restaurants, hotels, that sort of thing — along the Lehigh Towpath. My first attempt (you can find it here) basically looked at a region, within a mile (as the crow flies) of one section of the Lehigh River, and searching within that region for the amenities I was interested in. That was an interesting project in its own right, but, as I said in my earlier post, it didn’t really solve the right problem:  there are many places within a mile, or even a quarter mile of the river, that are not anywhere near accessible from the towpath: they may be on the wrong side of the river, say, or not near a towpath access point. To be considered accessible, the points of interest would need to be within a mile (or whatever arbitrary distance I end up choosing), by road, of an access point on the towpath.

    map of roads near lehigh river
    Data gathered, and almost ready for routing.

    I didn’t really have a plan to make this happen yet, but with or without a specific plan, I figured my first order of business was to get the information I would use. That previous analysis used Google Maps, but I felt that their data was a bit encumbered (in terms of my rights to it), and it seemed that Google didn’t play as well as I’d like with QGIS anyway, so I decided to use the data available through Open Streetmap, for both the road network and the set of amenities. (I already had a collection of the towpath access point locations left over from a previous experiment.) I got those sets of data, and massaged them so that I only had the parts that fell within a mile of the bike paths in the Lehigh Valley. This gave me the data seen to the right, where the aqua lines are the road network, the red lines are bike trails (the towpath, plus the Palmer Bike Path), the red stars are trail access points, and the orange dots are the amenities (restaurants, fast food etc).

    (One note about the road network: You probably can’t see it at this resolution, but I made a point of excluding roads that are not practical/legal/safe for cycling, like US-22, I-78 and a few others. There are also a number of places, like the New Street and Hill-to-Hill Bridges, where roads or the trail are connected via stairways to the bridges above; after our own struggles, a few years ago, with stairs and fully loaded touring bikes at the Ben Franklin Bridge, I decided to also exclude stairways from my network.)

    So that gets us the data, what about the analysis? My first thoughts were to see if I could find all the points on the road network that were a mile away from an access point, then connect the dots to define a region, and then find all the amenities within that region. My second thoughts were that this approach would put me back in the same situation as my first attempt, since I could easily find roads that were not reachable within that region, such as bridges. (Bridges became my nemesis for a while.) I eventually decided that my best strategy would be to find the shortest route between each access point and each amenity, and select from the amenities based on the lengths of the routes I found.

    To perform the actual routing analysis, I have three options:

    In terms of a learning curve, I have some experience with networks in GRASS, and I feel at least a little comfortable with Python (and copy-paste, with scripts I find online), so pgRouting will probably be the most difficult for me to pick up. Meanwhile, the Network Analysis library can use the data I already have, but Open Streetmap deals with road networks in a way that’s not directly compatible with either GRASS or pgRouting — their topological models are different, but that’s an issue for a future post. I would have to either re-import the road network to get it to work with pgRouting, or further process the one I have for GRASS.

    Each one of these approaches will be the subject of its own post. Given that the Python approach is not the hardest, and my data is already in the form I’d need, I am going to try my hand with the Network Analysis library first. Stay tuned for Part 2, whenever…


  • Befuddle Oneself Methodically

    I’ve been thinking a lot about the economic impact of the towpath lately, and have been looking to the Great Allegheny Passage as a model, where many food and lodging places have popped up to cater to the cycle touring crowd.  I know that the D&L Corridor people are also looking at various businesses and how the towpaths might impact them, but I believe that they are looking at it from a county-wide perspective, when they should probably really be looking at impact within a few blocks of the path, and or at least within about a mile of the river — would you decide to take a fully loaded touring bike miles out of your way, and probably up some hill to get away from the river, just for say, lunch, if you didn’t absolutely have to? So, that got me thinking about the question: what restaurants, hotels, bike shops, and other amenities are actually within a mile of the relevant sections of the  Delaware and Lehigh Rivers?

    (This was a good first approximation, but it’s surely a naive way of looking at the problem, since there are many places within a mile of the river as the crow flies, that are not actually within a mile, or maybe even many miles, of anyplace accessible from the towpath — and places that will see towpath business will need to be located within a matter of blocks, not miles, from towpath access points. But I realized all that later as I thought more about the overall situation, and my first analysis of towpath business prospects was what I worked on first.)

    The way I looked at it, my original problem broke down into two parts. First, what is the region within one mile (or whatever distance) from the river, and second, what are the amenities within that region? The first part was fairly straightforward, but the second, which looked like it would involve some kind of Google Maps search — and eventually it did — turned out to be more complicated than I thought…

    Partial map of Lehigh County and River
    Nifty QGIS: Five minutes of work to get the buffer zone.

    I used QGIS to deal with the first part. I took as my reference some Pennsylvania aerial photos, plus a property map of Lehigh County, and created a new line vector (in a projection that uses feet as a unit of measure), following what looked like the middle of the Lehigh River from about Laurys Station to just past Bethlehem, and then I used the “Create Buffer” geoprocessing tool to create a vector polygon buffer region around that section of river, whose distance from my line vector was 5280 feet, in other words, one mile. That part worked great, but what to do with my buffer region?

    My first thought was to take the buffer vector and export it to a KML file, import that KML file into a custom Google Map (using Google’s “My Maps” personal map creation/editing feature), and then “do something” with it. That all worked great as well, up until the “do something” part — the KML file, and the personal map, were not much use when it came to customizing a map search.

    I did find online, however, that there were some things you could do with polygonal regions, and testing locations (such as the ones returned from search results) to see if they fell within those regions, using the Google Maps API. This added two new steps: first I had to re-export my buffer region, this time as a GeoJSON file because that was what the API would accept, and I also had to sign up for an API key from Google Maps. Both of these were also straightforward and easy to do.

    The final step was to put it all together: make a web page, and (with some javascript), load and draw the GeoJSON file, run a search (for restaurants, in my experimental code), and then find and display results that fell within my region. Code code code, give it a try… nothing. I was able to load the file and see my region, but no place results would be displayed.

    Turns out, there is more than one Polygon type in the API, and the one created by loading a GeoJSON file is different than the one you can test locations against; I would have to convert my polygon from one form to another. (This seemed to me a bit much, especially since I thought I should have been able to load the original KML file and it would “just work.” After all, isn’t KML a Google thing, and kind of a standard?) No matter, the conversion process from one polygon to the other looked as straightforward as every other step so far, so I just added it to the end of the task chain. Code code code, give it a try… nothing, and here is where it started to get really frustrating.

    I couldn’t for the life of me figure out what was going wrong, it looked like I did things exactly the way I was supposed to but my new, converted polygon could not be made, and it looked like the original polygon actually was empty, even though it was drawn on screen. I eventually used a callback routine from the GeoJSON loading function to get the polygon coordinates, and for some reason that worked.

    That gave me my clue: the “some reason” was that the callback was not executed until after the file was done loading, so the conversion routine had something — a non-empty original polygon — to work with, while in my original code the rest of the script wouldn’t wait for the file to finish loading before continuing, so there really was nothing to work with yet when I tried to do the conversion. That took three paragraphs to write, but more than a day to work out…

    I didn’t really like my solution: if you’re forced to use callbacks like that, you end up going down the rabbit hole, callback after callback after callback, just to get some semblance of sequential execution. (Meantime, I found that some methods did not suffer from these kinds of problems, they seemed to wait for the data to load before trying to work on it. Strangely enough, all the simple API examples I found at Google used these methods instead of the one I needed.) Eventually I set up a wrapper function to hide the messy details and just get me my goddamned polygon from the goddamned GeoJSON file.

    Anyway, here is my demo map:

    UPDATE (7/26/2018): This map will stop working after July 30th, because the “radar search” function (see script below) has been deprecated by Google Maps. I may take some time to update the script — which I’ll mark with another update — but then again I may not, because this is a low-usefulness, low-visibility, low-priority thing for me, and also because fuck Google.

    And here’s my script. Most of this is based on Google Maps API examples, but the function getBuffer() loads the data, and createBufferPolygon() is the wrapper that creates the polygon object:

      var myNewMap;        // the google map
      var myPlaceService;  // object for google places api
      var myBufferPoly;    // the polygon that holds the buffer region
      var myInfoWindow;    // info window for the selected place
          
      // the callback function from loading the API, where everything actually happens
      function initMap() {
        myNewMap = new google.maps.Map(document.getElementById('map'), {
          center: {lat: 40.672628, lng: -75.422778 },
          mapTypeId: google.maps.MapTypeId.TERRAIN,
          zoom: 11
        });
            
        var bikeLayer = new google.maps.BicyclingLayer();
        bikeLayer.setMap(myNewMap);
        myBufferPoly = createBufferPolygon(
          'lbuf2.geojson',
          'lehigh',
          myNewMap,
          true,
          'green'
        );
        myPlaceService = new google.maps.places.PlacesService(myNewMap);
        myInfoWindow = new google.maps.InfoWindow();
    
        getSearch();
      }
          
      // this is the wrapper function, which calls the function that loads the GeoJSON file
      // after creating the polygon to hold the buffer region 
      function createBufferPolygon(url, featureName, map, isVisible, polyFillColor) {
        var bufPoly = new google.maps.Polygon({
          map: map,
          clickable: false,
          visible: isVisible,
          fillColor: polyFillColor
        });
        getBuffer(url, featureName, bufPoly);
        return bufPoly;
      }
          
      // this function loads a GeoJSON file containing a named polygon
      // then adds it to the given polygon object
      function getBuffer(url, featureName, poly) {
        var bufGeom;
        var bufferData = new google.maps.Data();
        bufferData.loadGeoJson(
          url, 
          {idPropertyName: 'name'},
          function(featarr) {
            bufGeom = bufferData.getFeatureById(featureName).getGeometry();
            poly.setPaths(bufGeom.getAt(0).getArray());
          });
      }
          
      // finds all restaurants within 15km of a certain location 
      function getSearch() {
        var request = {
          location: {lat: 40.703117, lng: -75.416561 },
          radius: 15000,
          keyword: 'restaurant'
        };
        myPlaceService.radarSearch(request, displayResults);
       }
          
        // displays search results that fall within the buffer region
        function displayResults(searchResults, searchStatus) {
          if (searchStatus !== google.maps.places.PlacesServiceStatus.OK) {
            console.error("getSearch error: " + searchStatus);
            return;
          }
          for (var i=0, result; result=searchResults[i]; ++i) {
            if (google.maps.geometry.poly.containsLocation(
              result.geometry.location, myBufferPoly)) {
                addMarker(result);
              }
            }
          }
          
      // adds marker for selected places
      function addMarker(place) {
        var marker = new google.maps.Marker({
          map: myNewMap,
          position: place.geometry.location,
          icon: {
            url: 'http://maps.gstatic.com/mapfiles/circle.png',
            anchor: new google.maps.Point(10, 10),
            scaledSize: new google.maps.Size(10, 17)
          }
        });    
        google.maps.event.addListener(marker, 'click', function() {
        myPlaceService.getDetails(place, function(result, status) {
          if (status !== google.maps.places.PlacesServiceStatus.OK) {
            console.error(status);
            return;
          }
          myInfoWindow.setContent(result.name);
          myInfoWindow.open(map, marker);
          });
        });   
      }
    

    That’s a lot of work for something that solves the wrong problem! My next look at this will likely just involve finding the access points, and doing Google searches near each one — soooo much less elegant…