## Archive for the ‘conversion’ tag

## MySQL Implementation of Google’s Encoded Polyline Algorithm

I just noticed someone created a MySQL implementation of Google’s Encoded Polyline Algorithm incorporating some elements of my PostgreSQL/PostGIS approach, specifically, support for multi-part, multi-ring polygon geometries.

And this is exciting, at the bottom of his post, Fabien says he’s working on a solution for consuming the Google encoded geometries in Leaflet. Nice! I love open source software!

I thought I’d mention the MySQL approach here to help expose Fabien’s solution to English searches and drive a little more traffic to his site. If you can’t read French, Google Translate is your friend!

FWIW–if anyone is interested and so inclined–there is a potential improvement that could be made to the Google Encoded Polyline implementation allowing users to specify a desired number of significant digits before rounding the Latitude and Longitude coordinate values prior to encoding. If memory serves, the documentation for Google’s algorithm allows/expects only 5 significant digits. But on some datasets with curved line features (often just heavily sampled line-segments), this limitation of 5 significant digits degrades the data, and you end up with jagged line features. So an improved solution would provide an optional parameter to specify significant digits.

All else equal–nice work, Fabien!

## Convert Google Maps Polygon (API V3) to Well Known Text (WKT) Geometry Expression

There’s dozens of reasons why you might want the Well Known Text (WKT) geometry expression for a Google Maps Polygon object.

Assuming you’re using the Google Maps API V3, and you’ve got a variable referencing your Polygon, I’ll suggest two approaches you can take to iterate over the paths and vertices in your Google Maps polygon and return the geometry expression as a Well Known Text string.

**Add a Simple Utility Method to Your Project**

Easy enough. Just add the following method to your project. Look below the method for an example of how you’d call it.

`function GMapPolygonToWKT(poly)`

{

// Start the Polygon Well Known Text (WKT) expression

var wkt = "POLYGON(";

```
```

` var paths = poly.getPaths();`

for(var i=0; i<paths.getLength(); i++)

{

var path = paths.getAt(i);

// Open a ring grouping in the Polygon Well Known Text

wkt += "(";

for(var j=0; j<path.getLength(); j++)

{

// add each vertice and anticipate another vertice (trailing comma)

wkt += path.getAt(j).lng().toString() +" "+ path.getAt(j).lat().toString() +",";

}

// Google's approach assumes the closing point is the same as the opening

// point for any given ring, so we have to refer back to the initial point

// and append it to the end of our polygon wkt, properly closing it.

//

// Also close the ring grouping and anticipate another ring (trailing comma)

wkt += path.getAt(0).lng().toString() + " " + path.getAt(0).lat().toString() + "),";

}

// resolve the last trailing "," and close the Polygon

wkt = wkt.substring(0, wkt.length - 1) + ")";

return wkt;

}

Here’s how you’d access the Well Known Text expression using the utility method:

`// Assuming you've already instantiated "myPolygon" somewhere.`

var wkt = GMapPolygonToWKT(myPolygon);

**Extend Google’s Polygon Object Prototype with a ToWKT() Method**

There’s nothing wrong with the first approach, but you might find it handy to extend Google’s Polygon object prototype, itself, to include a ToWKT() member function, which makes it even easier to get its Well Known Text. To do that, add the following JavaScript somewhere near the top of your code (caveat—this will need to be called after the Google Maps library has been loaded):

if (typeof google.maps.Polygon.prototype.ToWKT !== 'function') { google.maps.Polygon.prototype.ToWKT = function() { var poly = this; // Start the Polygon Well Known Text (WKT) expression var wkt = "POLYGON("; var paths = poly.getPaths(); for(var i=0; i<paths.getLength(); i++) { var path = paths.getAt(i); // Open a ring grouping in the Polygon Well Known Text wkt += "("; for(var j=0; j<path.getLength(); j++) { // add each vertice, automatically anticipating another vertice (trailing comma) wkt += path.getAt(j).lng().toString() + " " + path.getAt(j).lat().toString() + ","; } // Google's approach assumes the closing point is the same as the opening // point for any given ring, so we have to refer back to the initial point // and append it to the end of our polygon wkt, properly closing it. // // Additionally, close the ring grouping and anticipate another ring (trailing comma) wkt += path.getAt(0).lng().toString() + " " + path.getAt(0).lat().toString() + "),"; } // resolve the last trailing "," and close the Polygon wkt = wkt.substring(0, wkt.length - 1) + ")"; return wkt; }; }

If you prefer the second approach, you can get the Well Known Text expression like this:

`// Assuming you've already instantiated "myPolygon" somewhere.`

var wkt = myPolygon.ToWKT();

## PostGREsql/PostGIS Implementation of Google’s Encoded Polyline Algorithm

**[Edit 30 Jan, 2014]**

I added an additional PostGREsql method to perform Polygon encoding by concatenating polygon geometries (delimiter: †) and their inner rings (delimiter: ‡) together into one massive encoded block of ring features. I also provided an example JavaScript method demonstrating how to bring the amalgamated polygon feature encodings into your Google Map.

By some uncanny twist of the fates, I’ve elected to use, had to use, and/or been asked to develop applications that use Google Maps ASCII Encoded Polyline expressions. In previous encounters, I’ve used a PHP class to handle the encoding task, and most recently I wrote a Python method to decode these expressions so that I could return a 3rd-party’s encoded geometries to WKT and import them into a spatially aware database.

So far so good.

However one thing has always bugged me about using the PHP solution–I don’t like using a piece of middleware to handle what I consider to be a responsibility of the data layer. Mark McClure’s page, which is basically the seminal authority on this topic, provides external links to implementations for Perl, Ruby, PHP (note: I prefer the PHP class linked, above), Java, and Mathematica. Also, by searching Stack Overflow, you can find implementations of the algorithm in both C# and VB.Net. But for all my efforts searching, I could never dredge up an implementation for either MySQL or PostGREsql/PostGIS. Bummer.

Looking up, it seems version 2.2 of PostGIS might include some built-in Google encoding conversion methods. While this is cool enough for a hat tip, unfortunately, it’s too inconvenient to wait that long, and even then, there’s no guarantee the implementation will work the way I expect with complex Polygon geometries; for instance, maybe it will encode only the exterior ring of Polygons, ignoring MultiPolygons completely, etc. For that matter, it’s equally possible there could be some bugs. So with this said, and even though the previously-mentioned PHP implementation does the job, my boss was cool-enough to let me take a crack at implementing the algorithm as a PostGREsql/PostGIS function, and then share the results with the world. Since some initial testing confirms my PostGIS implementation works, I’ll just post the code parts and hope others find it useful.

For what it’s worth, if anyone finds a bug or has recommendations for improvements, please don’t hesitate to drop me a line.

**Sample query calling the first encoding function on the EXTERIOR RING of Polygon geometries:
(Also works on single-part LINESTRING features.)**

/************************************************************************ * Note that the encoding method can accept a LINESTRING only, which * is the geometry type used to represent the ring parts of a Polygon. * To help understand this, and why, please see the trailing discussion * section, which elaborates on this situation further. ************************************************************************/ SELECT GoogleEncodeLine(ST_ExteriorRing(wkb_geometry)) as Google FROM polygons_wgs84 WHERE ST_GeometryType(wkb_geometry) = 'ST_Polygon' LIMIT 10 ;

**[Added 30 Jan, 2014]**

**Sample query calling the second encoding function on Polygon and MultiPolygon geometries:
(Preserves multi-part polygons and their inner-ring parts, a.k.a. “holes”.)
**

/************************************************************************ * This encoding method will accept Polygon and MultiPolygon geom types. * The output returned is an amalgamation of Polyline encodings, where * individual geometries and their interior rings are concatenated * together using string delimiters, †, and ‡, respectively. ************************************************************************/ SELECT GoogleEncodePolygon(wkb_geometry) as GooglePolygon FROM polygons_wgs84 LIMIT 10 ;

**Implementation functions to execute/save in your PostGREsql instance:**

**[Added 30 Jan, 2014]**

/************************************************************* * Pass in either a Polygon or MultiPolygon geometry. Returns * an array of ASCII-encoded Polygon feature parts, including * multi-part geometries and their interior rings. ************************************************************/ CREATE OR REPLACE FUNCTION GoogleEncodePolygon ( g1 GEOMETRY ) RETURNS TEXT AS $$ DECLARE ng INT; -- Store number of Geometries in the Polygon. g INT; -- Counter for the current geometry number during outer loop. g2 GEOMETRY; -- Current geometry feature isolated by the outer loop. nr INT; -- Store number of internal ring parts in the Polygon. r INT; -- Counter for the current inner-ring part. r1 GEOMETRY; -- Exterior ring part isolated BEFORE the inner loop. r2 GEOMETRY; -- Inner-ring part isolated within the inner loop. gEncoded TEXT; -- Completed Google Encoding. BEGIN gEncoded = ''; ng = ST_NumGeometries(g1); g = 1; FOR g IN 1..ng BY 1 LOOP g2 = ST_GeometryN(g1, g); if g > 1 then gEncoded = gEncoded || chr(8224); END IF; -- Get ExteriorRing now; if there are any holes, get them later in the loop.. r1 = ST_ExteriorRing(g2); gEncoded = gEncoded || GoogleEncodeLine(r1); nr = ST_NRings(g2); if nr > 1 then -- One (1) is because interior rings is one-based. -- And nr-1 is because ring count includes the boundary. FOR r IN 1..(nr-1) BY 1 LOOP r2 = ST_InteriorRingN(g2, r); gEncoded = gEncoded || chr(8225) || GoogleEncodeLine(r2); END LOOP; END IF; END LOOP; RETURN gEncoded; End $$ LANGUAGE plpgsql;

/************************************************************* * First of two methods. Pass in a geometry (LINESTRING only). * Returns ASCII-encoded point array for use in Google Maps. ************************************************************/ CREATE OR REPLACE FUNCTION GoogleEncodeLine ( g GEOMETRY ) RETURNS TEXT AS $$ DECLARE pt1 GEOMETRY; pt2 GEOMETRY; p INT; np INT; deltaX INT; deltaY INT; enX VARCHAR(255); enY VARCHAR(255); gEncoded TEXT; BEGIN gEncoded = ''; np = ST_NPoints(g); IF np > 3 THEN g = ST_SimplifyPreserveTopology(g, 0.00001); np = ST_NPoints(g); END IF; pt1 = ST_SetSRID(ST_MakePoint(0, 0),4326); FOR p IN 1..np BY 1 LOOP pt2 = ST_PointN(g, p); deltaX = (floor(ST_X(pt2)*1e5)-floor(ST_X(pt1)*1e5))::INT; deltaY = (floor(ST_Y(pt2)*1e5)-floor(ST_Y(pt1)*1e5))::INT; enX = GoogleEncodeSignedInteger(deltaX); enY = GoogleEncodeSignedInteger(deltaY); gEncoded = gEncoded || enY || enX; pt1 = ST_SetSRID(ST_MakePoint(ST_X(pt2), ST_Y(pt2)),4326); END LOOP; RETURN gEncoded; End $$ LANGUAGE plpgsql;

/************************************************************** * Second of two methods. Accepts a signed integer (LON or LAT * by 1e5) and returns an ASCII-encoded coordinate expression. *************************************************************/ CREATE OR REPLACE FUNCTION GoogleEncodeSignedInteger(c INT) RETURNS VARCHAR(255) AS $$ DECLARE e VARCHAR(255); s BIT(32); b BIT(6); n INT; BEGIN e = ''; s = (c::BIT(32))<<1; IF s::INT < 0 THEN s = ~s; END IF; WHILE s::INT >= B'100000'::INT LOOP b = B'100000' | (('0'||substring(s, 28, 5))::BIT(6)); n = b::INT + 63; e = e || chr(n); s = s >> 5; END LOOP; e = e || chr(s::INT+63); RETURN e; End $$ LANGUAGE plpgsql;

**[Added 30 Jan, 2014]**

**JavaScript method demonstrating how to add Polygon encodings to a Google Map object:
(This client implementation works for either the single or the multi-part polygons.)
**

/************************************************************* * JavaScript! Pass-in an encoded text block created by either * PostGIS method, GoogleEncodePolygon() or GoogleEncodeLine(), * and render it in your Google Map object. If you don't want * the map to zoom to each rendering, just remove the "bounds" * variable and any references to it. ************************************************************/ function renderEncoded(encoded_path) { var bounds = new google.maps.LatLngBounds(); var $encodedGeoms = encoded_path.split("†"); for (var i=0; i<$encodedGeoms.length; i++) { var encodedGeom = $encodedGeoms[i]; var $encodedRings = encodedGeom.split("‡"); var polyPaths = []; for (var j=0; j<$encodedRings.length; j++) { var ptarray = google.maps.geometry.encoding.decodePath($encodedRings[j]); polyPaths.push(ptarray); } var polygonObject = new google.maps.Polygon( { paths: polyPaths, strokeColor: '#890000', strokeOpacity: 1.0, strokeWeight: 2 }); polygonObject.setMap(map); polygonObject.getPath().forEach(function(e) { bounds.extend(e); }); } map.fitBounds(bounds); }

**And some additional discussion..**

There are two “gotchas” when it comes to implementing the encoding algorithm with respect to Polygons:

1) Polygons, as geometries, can be composed of many rings. The outer ring is considered to be the boundary, and various inner rings are often called “holes”. So this is a specified, understood, and accepted built-in many-to-one relationship between polygons and their internal ring geometries.

And 2) It’s not rare to find polygon tables containing both Polygon and MultiPolygon data types. I think this happens because ESRI allows it, and so in an effort to play well with others, other GIS systems have accommodated it. At least, I know this is true for MySQL and PostGIS.

Here’s why this makes trouble–Google’s encoding algorithm is only intended to represent individual point arrays as a singular geometry. Basically, as long as your first point equals your last point, it’s considered to be a closed geometry, and you can add it and render it in a Google Map as a polygon. The algorithm itself isn’t designed to represent nested arrays, which would be necessary to render either a Polygon with “holes” or a MultiPolygon, which could potentially define many polygons with holes of their own! As such, I suspect there could be considerable disagreement as to how a Polygon-to-Google-Encoded method should actually handle Polygons..

The only solutions I can imagine for this problem would require “faking” a one-to-many relationship by perhaps delimiting together several encodings to account for MultiPolygons and/or single feature Polygons with interior rings. But this starts to get weird. So to keep things somewhat simple for the sake of the post, I chose to stay true to the algorithm’s intent and return a single encoded geometry expression from my method. And the sample query demonstrates this by calling the method against the outermost ring (i.e. the boundary) of a Polygon geometry type, which PostGREsql regards as a LineString, anyway.

**[Added 30 Jan, 2014]**

Since I wanted to handle the more complex geometries, I wrote the wrapper method `GoogleEncodePolygon()`

, to first iterate over `ST_NumGeometries()`

and gain access to any multi-part features, then second, iterate over `ST_NRings()`

using `ST_InteriorRingN()`

–you could also do this using `ST_DumpRings()`

–and gain access to any interior rings of the Polygon geometries, themselves. Then, for each ring part, I call `GoogleEncodeLine()`

, and concatenate together all those expressions into one massive block of “compound” expressions. I chose to delimit each geometry encoding using an extra-special character that would never be used by Google’s algorithm; for example `chr(8224)`

, which corresponds to “†”. I then further delimit the internal ring parts using another special character, `chr(8225)`

, which corresponds to “‡”, and return all these concatenated together as a compound encoding expression. Then, on the client-side (a JavaScript example is provided above), I merely split the compound expression against my delimiters, loop over the various expressions, and add them to the map individually. Note if you are attaching attributes to your features, you’ll need to remember to include them explicitly to each unique Polygon added to your map.

## Decode Google Map encoded points as Well Known Text (WKT) with Python

I had close encounter of the 5th kind yesterday.. here’s the gist..

It started when someone “gave” me a GIS dataset (..of polygons, kind of..) that a colleague of theirs, way back in ancient history, chose to pre-cook as ASCII-encoded point pairs. Their intention was almost certainly to use the pre-cooked data in Google Maps. Anyway, being arguably sane, I wanted to return this data to a more GIS-normal format so I could put it in a database like MySQL or Post and use it for other stuff.

I considered a few different approaches to this problem, including creating a Google Map that could load-in all of the polygons from their encodings, then iterate over the polygons, interrogate the polygon point pairs, and finally concatenate WKT features from the points and save the geofeatures into a MySQL table. This approach offered the advantage of using Google’s existing Maps API to do the decoding for me. But let’s face it, that’s lame, uninspired, and not inventive.. it wasn’t even interesting.

Besides.. I wanted to use Python.

I expected to find a Python recipe for this looming in the misty www, but I didn’t. However, I did find a JavaScript recipe by trolling around in Mark McClure’s website. Specifically, he provides a Polyline Decoder utility, and when I viewed the page source, I found the JavaScript code that actually does the decoding (opening in FireFox will show you the code, or IE should prompt you to download the file).

Long story short, the following Python method is an adaptation of Mark McClure’s JavaScript method (twisted a little to return WKT features rather than the pure point array). If you’re somewhat comfortable with Python, you should be able to copy/paste the method right into your Python file and start calling it; just pass-in the encoded point string and let the method do the rest.

Best / Elijah

———

def decodeGMapPolylineEncoding(asciiEncodedString): print "\nExtrapolating WKT For:" print asciiEncodedString strLen = len(asciiEncodedString) index = 0 lat = 0 lng = 0 coordPairString = "" # Make it easy to close PolyWKT with the first pair. countOfLatLonPairs = 0 firstLatLonPair = "" gotFirstPair = False while index < strLen: shift = 0 result = 0 stayInLoop = True while stayInLoop: # GET THE LATITUDE b = ord(asciiEncodedString[index]) - 63 result |= (b & 0x1f) << shift shift += 5 index += 1 if not b >= 0x20: stayInLoop = False # Python ternary instruction.. dlat = ~(result >> 1) if (result & 1) else (result >> 1) lat += dlat shift = 0 result = 0 stayInLoop = True while stayInLoop: # GET THE LONGITUDE b = ord(asciiEncodedString[index]) - 63 result |= (b & 0x1f) << shift shift += 5 index += 1 if not b >= 0x20: stayInLoop = False # Python ternary instruction.. dlng = ~(result >> 1) if (result & 1) else (result >> 1) lng += dlng lonNum = lng * 1e-5 latNum = lat * 1e-5 coordPairString += str(lonNum) + " " + str(latNum) if gotFirstPair == False: gotFirstPair = True firstLatLonPair = str(lonNum) + " " + str(latNum) countOfLatLonPairs += 1 if countOfLatLonPairs > 1: coordPairString += "," # The data I was converting was rather dirty.. # At first I expected 100% polygons, but sometimes the encodings returned only one point. # Clearly one point cannot represent a polygon. Nor can two points represent a polygon. # This was an issue because I wanted to return proper WKT for every encoding, so I chose # To handle the matter by screening for 1, 2, and >=3 points, and returning WKT for # Points, Lines, and Polygons, respectively, and returning proper WKT. # # It's arguable that any encodings resulting in only one or two points should be rejected. wkt = "" if countOfLatLonPairs == 1: wkt = "POINT(" + coordPairString + ")" elif countOfLatLonPairs == 2: wkt = "POLYLINE(" + coordPairString + ")" elif countOfLatLonPairs >= 3: wkt = "POLYGON((" + coordPairString + "," + firstLatLonPair + "))" return wkt