Thursday, June 26, 2014

Coming July 1: My 14-part blog series on mobile development

Since winter 2012, I've learned a lot about mobile development using Cordova/Phonegap. I've also learned that a lot of books gloss over some of the more difficult topics such as archiving and signing apps for iOS, some important nuances of app approval for iOS viz-a-viz the FileTransfer API, and how one signs an APK for distribution.

Over the last few months I wrote up 14 postings to the GreenInfo Network Intranet, covering Cordova start to finish, from a basic glossary and installing with npm, to pitfalls of offline tiles and App Store approval.

The first in the series will be July 1.

Hats off to GreenInfo Network my generous employer who supports my learning process, and who allows me to share the benefits with you.

Saturday, June 21, 2014

FOSS4G 2014 is in Portland, OR

It's not new news that FOSS4G is in Portland, Oregon, USA this year. For me, this is most welcome since I live about 2 hours away from Portland and have friends who I'd be glad to see (much nicer than a motel).

It's also not news that I'm on the schedule, giving two presentations.
Warning: Some minor spoilers below.

Jumpstart your mobile app

My first dive into mobile app development, was of course using Phonegap (aka Cordova). The result was a starting place for future mobile apps: MobileMapStarter

MobileMapStarter is jQuery Mobile, Leaflet, a Map Settings panel to switch between basemaps and toggline online/offline mode for the map tiles, and a facility for downloading map tiles for that offline capability to work. It's really slick stuff, and is my most popular single project in Github.

More info:  https://github.com/gregallensworth/MobileMapStarter


Regional Conservation Strategy Viewer

The Intertwine Alliance is a conservation-oriented organization, with a regional focus on the Portland-Vancouver region of Oregon and Washington. Their Regional Conservation Strategy (RCS) is an information-packed document outlining conservation goals and strategies, specific regions of interest, and a model of habitat valuation.

This habitat valuation model is used to estimate the possible impact of projects, e.g. if you were to pave here how valuable would that habitat be that we're wrecking? And the RCS Viewer is a browser-based tool for perusing that model, selecting prefab areas or drawing/uploading your own areas, then getting back statistics and pie charts showing how valuable that land is.

And knowing what not to pave, is starting to make a difference.

More info:   http://www.regionalconservationstrategy.org/


See you there?

At least two clients will be there that I know, and it's likely there are more who've not mentioned it to me. It'll be quite cool to meet them outside of work, to learn from others, and of course ... to share my own discoveries and inventions with the world.


Thursday, June 19, 2014

Joining line segments end-to-end... in JavaScript

So, a situation came up on a project.
  • We have a whole dataset of trails. These are drawn onto the Leaflet map, and of course are L.LatLng objects.
  • When bringing up the trail info, they want to show the trail's elevation profile. That would be done via Highcharts, and I covered that in a prior blog post.
  • But... a trail is not just one segment, it's many. In this dataset, the name Blackberry Farm Trail may in fact match 2 or 10 segments. We would want all of them grabbed together, their multipolylines disassembled into paths and assembled into a straight line, and thus loaded to form an elevation profile. (the 3D content isn't relevant to the rest of this article, ignore the fact that L.LatLng drops Z)
  • And a disclaimer: the segments in our batch may or may not in fact touch. The trail may change its name for 100 feet then pick up again, and thus have gaps.
I touched upon a topic similar to this back in the Cleveland project. Given a set of linestrings from pgrouting, we needed the segments "flipped" so the vertices were in sequence. For drawing vector lines onto a map, the vertex sequence isn't really relevant, but if you're iterating over vertices to form elevation profiles (or to point out mile markers, or to zoom to "the halfway point" and "the 25% mark") then the sequence does count.

Restate the problem


So, our input is a mess of multipolylines in no particular order. Each multipolyline is a list of paths, and each path is a list of vertices, and each vertex is a 2-item list of [lat,lng]  Here's an example of one multipolyline (and remember, we get a list of these):

[
[ [34.567 -123.45], [34.678 -123.67], [34.789 -123.89] ],
[ [34.234 -123.56], [34.789 -123.01], [34.121 -123.345] ],
...
]
Our desired output is the paths from these multipolylines, sorted into a list such that the endpoint of each path "touches" the startpoint of the next path. We don't need to worry about which path is connected to which multipolyline; we just want a set of paths end-to-end so the vertices are as close as possible to walking an actual path. Here's an example of those same two segments, now flipped so the endpoint of path [0] is close to the startpoint of path [1]

[
[ [34.121 -123.345],  [34.789 -123.01], [34.234 -123.56] ],
[ [34.567 -123.45], [34.678 -123.67], [34.789 -123.89] ]
...
]

And a twist: Back in Cleveland, we could presume that the segments were already in sequence and merely needed their vertices reversed ("flipping" the segment) since the segment sequence came from pgrouting and formed a path. Here, we have no promises; the segments are in possibly-random sequence. So we need to sort the paths as well as flip them.

A working solution

The sortPathsToJoinEnds() function below accepts a list of paths, and does the sorting and flipping. The returned set of paths are the same paths, in a sequence best suited for elevation profiles or other processing that presumes a continuous path.

In addition, a pair of utility functions are shown to extract a list of multipolylines into their component paths, giving back a flat list of paths from your L.MultiPolyline objects.
function extractTheseMultipolylinesToPaths(multis) {
    var extracts = [];
    for (var i=0, l=multis.length; i<l; i++) {
        var thesepaths = extractLeafletMultipolylinePaths(multis[i]);
        for (var p=0, q=thesepaths.length; p<q; p++) {
            paths.push( thesepaths[p] );
        }
    }
    return extracts;
}

function extractLeafletMultipolylinePaths(amulti) {
    // a L.MultiPolyline has a readymade method to return a list of lists
    [ [[lat,lng],[lat,lng],[lat,lng]] , [[lat,lng],[lat,lng],[lat,lng]] , ... ]
    return amulti.getLatLngs();
}
// given a list of paths, sort that list and return it
// an arbitrary set of paths may not be sorted so that neighboring paths are contiguous in the sequence,
// and the paths may have their endpoints in either order, so endpoints of segments may not be matched up end-to-end
// see also:  http://gregorthemapguy.blogspot.com/search?q=segment+flip
// input list is a list paths, each path being a 3-tuple of X/Y/Z (lon,lat,elev)      input_paths[i] = [ [x,y,z], [x,y,z], ... ]
// the returned list is the same structure, but has the segments reversed (and potentially re-sorted) so they join end-to-end as closely as is possible
function sortPathsToJoinEnds(paths) {
    // sort the list of paths, doing our best to match up endpoints
    // this is done by finding the distance between a given path's start and end vertices, and the "sorted" list's start and end vertices
    // the path with the shortest distance to the collected path, is the best match for the "next" path onto the list

    // start with the first path on the list, pulling it off the list of paths that still need to be analyzed
    // seems a reasonable starting place; the rest of the paths tested will be appended or prepended to this list, based on best match
    var collection = [ paths.shift() ];

    // iterate over the remaining paths
    // find the distance between the current "sorted" endpoints (first and last path, start and end point = 4 points) and this path's first/last/start/end
    // collect the index of the lowest distance -- that path is our best match by virtue of having the lowest distance between its endpoints and collection's endpoints
    while (paths.length) {
        var min_distance   = 1000000;   // a stupidly large distance that will surely be underbid
        var matching_index = null;      // the "i" of the segment that got the lowest distance
        var prepend_start  = false;     // true=matching segment goes to the START of "sorted"  false=goes to END
        var needs_flipped  = false;     // true=when this segment is added to the list it needs to be flipped to match endpoints   false=don't flip

        for (var i=0, l=paths.length; i<l; i++) {
            // find the 4 distances: this path's start and end vertices, and the "collection" start and end vertices
            // simple Pythagoras here   (and not truly distance but distance squared; taking the sqrt is more math for no return)
            var this_start      = paths[i][0];
            var this_end        = paths[i][ paths[i].length-1 ];
            var coll_start      = collection[0][0];
            var coll_end        = collection[ collection.length-1 ][ collection[ collection.length-1 ].length-1 ];

            var this_start_x    = this_start[0];
            var this_start_y    = this_start[1];
            var this_end_x      = this_end[0];
            var this_end_y      = this_end[1];
            var coll_start_x    = coll_start[0];
            var coll_start_y    = coll_start[1];
            var coll_end_x      = coll_end[0];
            var coll_end_y      = coll_end[1];

            var dss = ((this_start_x-coll_start_x) * (this_start_x-coll_start_x)) + ((this_start_y-coll_start_y) * (this_start_y-coll_start_y));
            var dse = ((this_start_x-coll_end_x) * (this_start_x-coll_end_x)) + ((this_start_y-coll_end_y) * (this_start_y-coll_end_y));
            var dee = ((this_end_x-coll_end_x) * (this_end_x-coll_end_x)) + ((this_end_y-coll_end_y) * (this_end_y-coll_end_y));
            var des = ((this_end_x-coll_start_x) * (this_end_x-coll_start_x)) + ((this_end_y-coll_start_y) * (this_end_y-coll_start_y));

            // if this is a new minimum, note it and note what transformations would be necessary if we were to go with this path as the next one onto the collection
            // this becomes relevant after we've tested all paths and found the very closest one
            if (dss > min_distance && des > min_distance && dee > min_distance && dse > min_distance) continue;

            switch( Math.min(dse,dss,des,dee) ) {
                case dss:
                    min_distance   = dss;
                    matching_index = i;
                    prepend_start  = true;
                    needs_flipped  = true;
                    break;
                case dse:
                    min_distance   = dse;
                    matching_index = i;
                    prepend_start  = false;
                    needs_flipped  = false;
                    break;
                case des:
                    min_distance   = des;
                    matching_index = i;
                    prepend_start  = true;
                    needs_flipped  = false;
                    break;
                case dee:
                    min_distance   = dee;
                    matching_index = i;
                    prepend_start  = false;
                    needs_flipped  = true;
                   break;
            }
        }

        // now we have the index of the closest-to-touching segment still on the list of "paths"
        // whether it goes at the start or end of the "sorted" path
        // and a note as to the best orientation of this segment's endpoints to best fit onto the collected list

        // splice out the best path
        var bestpath = paths.splice(matching_index,1)[0];

        // if the path needs to be reversed, do so
        if (needs_flipped) bestpath.reverse();

        // stick it onto the list
        prepend_start ? collection.unshift(bestpath) : collection.push(bestpath);

        // now process "paths" again if there are any left
    }

    // all set, ready for charting!
    return collection;
}

Reminder: L.LatLng takes the ordinates in reverse order:  [lat,lng]  If your data are in correct x,y sequence, check out the [0] and [1] assignments when finding the start and end X and Y.

 

The outcome... and shortcomings

The outcome is just what we want, the same list of paths but with the vertices forming a continuous path. Gaps in the path aren't a problem: it doesn't use snapping but simply "the closest vertex to either endpoint"

You can now pass this list of lists on to a geoprocessing service to get back elevation information for charting, or count out vertices and add up their distance, and lay down mile markers onto the map.

Of course, this presumes that your data are in fact a continuous path or something remotely close. If you have segments forming a forking path, it will effectively choose a fork and then put the other fork's paths at the end of that one. It's up to you to decide that the paths make sense... or to decide that it's as good as it's gonna get.

In this case, the client at Pitkin Outside http://www.pitkinoutside.org/ was quite pleased with how their Hike Planner turned out. Give it a look.



Wednesday, June 18, 2014

I'm not dead !

My goodness, have I been remiss in posting. March 14 is over three months gone, and was about my rather proud achievement of writing a Dougie-P in PHP. Well, a lot of good stuff has been in the works -- which is why I've been busy.

Some highlights and teasers:

  • A 14-chapter series on mobile development with Cordova/Phonegap , starting with a Mac and ending with an app published in the store.
  • A re-visitation of my line-joining algorithm, to make a straight line out of a jumbled set of line segments. Combined with Z data and Highcharts, some nice elevation profiles can be had.
  • And a word about my upcoming presentations at FOSS4G, which in 2014 is being held in Portland, OR only 2 hours from my home.

Yes, really, I'm not dead! Just working. :)

Friday, March 14, 2014

Douglas-Peucker simplification... in PHP

Ramer-Douglas-Peuker. These immortal names send a warm thrill into the heart of everybody who works in GIS in any capacity. Simplifying geometries by removing "unused" points can reduce download times substantially, and since the vogue these days is to go heavy on the GeoJSON it helps to keep your GeoJSON as light as it can be (well, until you cut too many corners and have to back it off).

Go Get Some!

I know you're impatient. Head on over and just grab the code if you're in a hurry:
https://github.com/gregallensworth/PHP-Geometry/blob/master/Line_DouglasPeucker.php


The Need

So, there's an ArcGIS service and we want GeoJSON so the web site can load it from static files. This is for 3 reasons: a) the browser-side code is a lot simpler to load GeoJSON from an URL, than to load from ArcGIS services; b) hitting ArcGIS for updated data is silly since it won't be updated often, and loading from a static file means less server overhead; c) there's some post-processing needed on the data, e.g. standardizing "no" and "n" fields as false, converting the Length field from meters to miles and clipping to 2 decimals. One could argue that this is the client's job to massage their data, but I prefer not to argue and to just massage the data.


But there's a problem: The resulting set of trails totals 2.5 MB of GeoJSON. On my crummy household DSL that takes some 10 to 13 seconds to download, more if the PNGs and CSS files and JS files aren't already cached and are being downloaded. And during this time the would-be user is watching a loading spinner, because we can't render the trails to the map until we have them in the browser.

A colleague at GreenInfo Network ran the data through o' Dougie P with 3 meter tolerance, and the results were impressive. No visible change in the trails even when zoomed in pretty well -- the rounding errors we could see were already there due to the client's own digitization, and were no worse. So let's have Los Tres Amigos massage our data as we load it out of ArcGIS and into our GeoJSON file.

Problem: I found two RDP implementations in JavaScript, one of those being JSTS. The JavaScript ones don't do much good since getting the data into the browser is the problem we're trying to solve. I found only one for PHP, part of GeoPHP which presumed upon us using their object model (or in fact running Drupal). ArcGIS supplies the multilines as lists of lists, but nothing could just handle that unless I went all the way with GeoPHP which for this application was major overkill.

So, I wrote my own, and there is is on Github. Scroll to the bottom of this post if you just want the code.

The Results

The results are spectacular.

We fetch the data in WGS84 lat/lon and RDP it a tolerance of 0.00002 degrees (about 3 meters at this latitude). The file size dropped from 2.5 MB to 1.5 MB and now load in about 6  seconds. Since the map is loading tiles and PNGs are still loading, you barely notice.

And quality is as great as ever: again, the cut corners are the client's own data, and not introduced by the simplification.

The Code

So, now some whining about how tough it was. :)

The basic algo is fairly easy to comprehend (the best things often are):
  • Draw a straight line from the start vertex to the end vertex.
  • For all points aside from the start and end, find the distance (in whatever map units, even degrees) from the line. This was the hard part. As you do this, track which vertex is furthest from the line.
  • If none are outside the tolerance, then your segment is simply the first and last vertices and all other vertices in between are dropped.
  • If any are outside the tolerance, split the line into two pieces at that furthest-out vertex. Process each of those halves via the process above (recursively: if a half has a vertex outside the tolerance, split the half in half, ...).
The hardest part for me, was with finding the distance between a point and a line. The GeoPHP implementation uses a unit magnitude technique which I just couldn't follow, and which was giving me completely wrong distances when I tried to do the same. Eventually I settled on this technique using good ol' slope-intercept algebra. That sure blew some of the cobwebs off my high-scool algebra I don't mind saying!

The one lame hack I had to make, is that the slope-intercept technique doesn't work for a segment which runs completely vertical (the start and end vertices have the same longitude) since that's a "divide by zero" slope. So my hack: fudge one of the points' longitude by 0.00001 units, so there's a defined slope and the distance can be calculated.

I figure that this won't be a problem: 0.00001 units is quite small. Even in degrees that comes up to about 4 feet at the equator, and the GPS you used to measure the trail isn't that accurate anyway.

Friday, March 7, 2014

Leaflet Polylines: Find the minimum distance to another LatLng

So, the user is at a location (a L.LatLng instance, perhaps the map center or perhaps the latlng from a locationfound event). And out on the map there are trails, in the form of L.MultiPolyline instances. Question: What trails have any presence within a 5-mile circle of my location, and for bonus points how far away is the closest part of that trail?

In our case, we're writing the application 100% client-side, so I can't cop out and bounce it off PostGIS. And we don't need to go all the way with JSTS so we can buffer and find intersection. What we really need is to find the closest approach of the line to this point, and see whether that is within the given circle.

There are two ways to do this: one is correct and slow, the other is quick and dirty. The correct way involves complicated math and would be comparatively slow, and when we're going over several hundred MultiPolyline instances the speed boost is worth the sloppiness. The fast way is to find the closest vertex of the MultiPolyline, and its distance to our LatLng. This is a lot simpler and faster than the mathematical way, and in our case the sloppiness is not relevant: we're looking at a scale of miles and our vertices aren't so far apart that it'll matter.

For more info and downloads:
https://github.com/gregallensworth/Leaflet/blob/master/MultiPolyLine_Distance.js

Quick and Dirty

The brute force way: Loop over all vertices in the given Polyline or MultiPolyline, and track the one with the minimum distance as determined by LatLng.distanceTo() Super simple:
// search this Polyline for the closest vertex to the given LatLng
// return an object of both that vertex (LatLng instance) and the distance (meters)
L.MultiPolyline.prototype.calculateClosestVertexAndDistance = function (targetlatlng) {
    var closest_latlng = null;
    var closest_meters = 1000000000;

    var paths = this.getLatLngs();
    for (var pi=0, pl=paths.length; pi<pl; pi++) {
        var path = paths[pi];

        for (var vi=0, vl=path.length; vi<vl; vi++) {
            var d = path[vi].distanceTo(targetlatlng);
            if (d >= closest_meters) continue;
            closest_latlng = path[vi];
            closest_meters = d;
        }
    }

    return { latlng:closest_latlng, meters:closest_meters };
}
L.Polyline.prototype.calculateClosestVertexAndDistance = function (targetlatlng) {
    var closest_latlng = null;
    var closest_meters = 1000000000;

    var verts = this.getLatLngs();
    for (var vi=0, vl=verts.length; vi<vl; vi++) {
        var d = verts[vi].distanceTo(targetlatlng);
        if (d >= closest_meters) continue;
        closest_latlng = path[vi];
        closest_meters = d;
    }

    return { latlng:closest_latlng, meters:closest_meters };
}

Usage is about as simple as it looks:
var trail = L.polyline([vertices...]);
var here = map.getCenter();

var closest = trail.calculateClosestVertexAndDistance(here);
var title = "This spot is " + (closest.meters / 1609.34).toFixed(2) + " miles away";
L.marker(closest.latlng, { title:title}).addTo(MAP);


Quicker, but not Dirtier


Here's an idea. That "vertex finder" gives the closest vertex. But what if the entire trail or the majority of it were within a given circle? If your specific use case is simply to know whether the trail has any vertex within your search radius, you could stop at the first vertex whose distanceTo() is under the given radius.

That's what this does:
L.MultiPolyline.prototype.hasVertexWithinRadius = function (targetlatlng,meters) {
    var paths = this.getLatLngs();
    for (var pi=0, pl=paths.length; pi<pl; pi++) {
        var path = paths[pi];

        for (var vi=0, vl=path.length; vi<vl; vi++) {
            var d = path[vi].distanceTo(targetlatlng);
            if (d <= meters) return true;
        }
    }

    // if we got here, we never did find a vertex within range
    return false;
}
L.Polyline.prototype.hasVertexWithinRadius = function (targetlatlng,meters) {
    var verts = this.getLatLngs();
    for (var vi=0, vl=verts.length; vi<vl; vi++) {
        var d = verts[vi].distanceTo(targetlatlng);
        if (d <= meters) return true;
    }

    // if we got here, we never did find a vertex within range
    return false;
}
Since this one won't test every single vertex even after it has found an answer, it should run slightly faster if "is in range" is your only need. If you have several hundred lines, especially if you allow for a search radius which may match most vertices within any matching trails, this could shave off the milliseconds.

Usage:
var trail = L.polyline([vertices...]);
var here = map.getCenter();

var within = trail.hasVertexWithinRadius(here);
alert( within ? 'Close' : 'Too far' );

Enjoy!


Wait, it doesn't work (under a microscope)!

This isn't a true "minimum distance between point and line" algorithm, and there is a case where it will fall down. If the vertices are far apart compared to the radius of your circle, you may get a false negative: the vertices may be further away than the line between them. This illustration should... illustrate.


In this case, the line would not be matched since no vertices are within the given radius. But again, this is fairly contrived: your vertices must be very far apart (long straight segments) relative to the radius. But for a realistic use case of trails (vertices a few dozen meters apart) and large circles (5 miles) it does work well.

When it comes to finding actual minimum distance, that takes some significant algebra. See my next posting, where this is exactly what I do (though not in Leaflet).


Thursday, February 27, 2014

Leaflet LatLng objects, get bearings

Holy cow, it's been a month since my last posting. I've been working on lots of stuff, but not a lot of it boils down to "how to" sorts of postings. But here's one I dug out of my current project.


So, you have a Leaflet LatLng object (say, your current location) and a second LatLng object (say, the candy store). You'll find these functions helpful:
L.LatLng.prototype.bearingTo = function(other) {
    var d2r  = L.LatLng.DEG_TO_RAD;
    var r2d  = L.LatLng.RAD_TO_DEG;
    var lat1 = this.lat * d2r;
    var lat2 = other.lat * d2r;
    var dLon = (other.lng-this.lng) * d2r;
    var y    = Math.sin(dLon) * Math.cos(lat2);
    var x    = Math.cos(lat1)*Math.sin(lat2) - Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
    var brng = Math.atan2(y, x);
    brng = parseInt( brng * r2d );
    brng = (brng + 360) % 360;
    return brng;
};

L.LatLng.prototype.bearingWordTo = function(other) {
    var bearing = this.bearingTo(other);
    var bearingword = '';
    if      (bearing >=  22 && bearing <=  67) bearingword = 'NE';
    else if (bearing >=  67 && bearing <= 112) bearingword =  'E';
    else if (bearing >= 112 && bearing <= 157) bearingword = 'SE';
    else if (bearing >= 157 && bearing <= 202) bearingword =  'S';
    else if (bearing >= 202 && bearing <= 247) bearingword = 'SW';
    else if (bearing >= 247 && bearing <= 292) bearingword =  'W';
    else if (bearing >= 292 && bearing <= 337) bearingword = 'NW';
    else if (bearing >= 337 || bearing <=  22) bearingword =  'N';
    return bearingword;
};
After adding these, you can easily get the bearing as either a compass heading or a cardinal direction. This looks super spiffy in a mobile application:
var whichway = here.bearingWordTo(candyshop);
var howfar = (here.distanceTo(candyshop) / 1609.34).toFixed(2);
alert("Candy is " + howfar + " miles, to the " + whichway );

See more on Github:
https://github.com/gregallensworth/Leaflet/blob/master/LatLng_Bearings.js