Cartometric Blog

A scrapbook of GIS tricks, with emphasis on FOSS4G.

ArcPy CalculateField_management Complex Python Expression to Return a String

without comments

My coworker was having some trouble implementing a CalculateField expression in Python so I wanted to document the solution for anyone with a similar issue. (Note: If you  just want the answer, see the What ultimately worked.. heading, below.)

Being unfamiliar with ArcPy, I was intrigued by this idea of defining a function as a string value, then passing that string as a variable into another function. I can see how that opens a powerful door—but it also strikes me as ultra weird, and because Python has some very particular spacing/indentation rules (as compared to say, JavaScript), I figured this technique would be ultra-prone to syntax issues ..and thus could be especially difficult to troubleshoot. So in this case, I also wanted to demonstrate how I ultimately thought through the problem.

After some quick Googling, we found two relatively useful pages of documentation in the ESRI support ecosystem, and both pages demonstrated different ways to formulate the function-as-string parameter necessary for the CalculateField_management() function.

The biggest difference between these examples is how they portrayed line breaks in the function-as-string; we’ll consider them first.


The first example..

In the first ESRI example, the function-as-string was concatenated together across multiple lines using a backslash \ like this..

codeblock = "def getclass(area):\
        if area <= 1000:\
            return 1\
        if area > 1000 and area <= 10000:\
            return 2\
            return 3" 

print codeblock

Realizing Python is very picky about syntax, I wanted to see what happened if I print codeblock  to the console–notice that spacing is preserved, but there are no line breaks:

Frankly, it’s difficult for me to believe this was ever a working example, and surprise—this approach didn’t end up working for us..


And the second example..

In the second ESRI example, the function-as-string was concatenated together across multiple lines using triple-quotes """ to declare a multiline string, like this..

codeblock = """def getClass(area):
    if area <= 1000:
        return 1
    if area > 1000 and area <= 10000:
        return 2
        return 3"""

print codeblock

Similarly, I wanted to codeblock to the terminal to see how the string was formatted before going into the ArcPy function; this time, the linebreaks were preserved:

..that looks more like a Python function, so that was a good start. But things still weren’t working..


How to test our custom function..?

We were still getting errors, and the stack dump claimed they were syntax errors. But because of the unusual nature of this software design, I wasn’t sure if they were really syntax errors, or some kind of misleading, catch-all error ..perhaps our code-in-a-string function was flawed?

So my next step was to test our function, as a legit Python function, to see if it would accept and return a value like I expected.

def getClass(area):
    if area <= 1000:
        return 1
    if area > 1000 and area <= 10000:
        return 2
        return 3

print getClass(1500)

As you can see, the code ran without any syntax errors..

Also, 2 is the correct return value for the input we used in the test. So I felt comfortable that the function itself—as written—worked. So now I was perplexed as to why the ESRI CalculateField_management() function was rejecting it ..with a syntax error.


What ultimately worked..

My co-worker’s function had another gotcha—unlike the two ESRI examples, which return numeric values, we needed our function to return a string. My Python was rusty, and I was forgetting if a string could be declared equally between single quotes i.e. 'MyString' and double quotes, i.e.  "MyString". Plus, I still new we were looking for an alleged syntax error.

I finally Googled a third page in the ESRI documentation and noticed this little tidbit in some ESRI documentation:

Python enforces indentation as part of the syntax. Use two or four spaces to define each logical level. Align the beginning and end of statement blocks, and be consistent.

My coworker had originally used one and two spaces, respectively, to control indentation in his custom function, along with the first technique for declaring a multiline string. So I finally rewrote our function using the triple-quotes technique as well as 2 and 4 space indents like ESRI said to, and this is what we ended up with..

For the record, while a large if/elseif block may not be the best way to write this function, it is a very readable way to write this function, not to mention the most intuitive from my coworker’s perspective.

codeblock = """def getclass(SPEEDLIM):
  if SPEEDLIM >= 60:
    return 'A10'
  elif SPEEDLIM == 55:
    return 'A30'
  elif SPEEDLIM == 50:
    return 'A20'
  elif SPEEDLIM == 45:
    return 'A30'
  elif SPEEDLIM == 40:
    return 'A30'
  elif SPEEDLIM == 35:
    return 'A40'
  elif SPEEDLIM == 30:
    return 'A00'
  elif SPEEDLIM == 20:
    return 'A00'
  elif SPEEDLIM == 15:
    return 'A00'
  elif SPEEDLIM == 10:
    return 'A61'
  elif SPEEDLIM == 5:
    return 'A62'
  elif SPEEDLIM <= 5:
    return 'A71'"""

And this worked. I think the fix resulted from using triple-quotes to declare the function string, rather than the backslash approach we started with. While we also changed the indentation to use 2/4-spaces, as mentioned above, it’s possible we unwittingly fixed a different syntax error in the function along the way. Since we didn’t make these changes independently, it’s impossible to know which exactly made the difference.



So this is a summary that might help you if you’re having difficulty with this kind of ArcPy data management exercise:

  • Try writing a basic version of your function in a throw-away .py file and test it to make sure it definitely works—if nothing else, this is a good sanity check.
  • Make sure to use the triple-quote technique to properly preserve line breaks, and declare your function as a string variable that you can pass into your call to the geoprocessor. (You may need to use single-quotes to declare any strings within your larger, multiline string, which was the technique we used.)
  • If things still aren’t working, try using increments of 2 spaces to control various levels of indentation within your custom function, i.e. 2/4/6/8/etc..

Hopefully that helps you get where you’re trying to go.


- elrobis

Written by elrobis

November 25th, 2014 at 3:50 pm

Posted in ArcPy,Conversion,Python

Create UTFGrid Tiles from PostGIS Tables

without comments

I assume you’re where I was about a week ago. That is, you’ve heard of UTFGrid, and now you want to render your own UTFGrid tiles. Perhaps you found yourself looking at this thread over at GIS.SE, but you don’t want to jump into TileStache just now. Anyway if you’ve got a working installation of GDAL/OGR and Mapnik 2+, complete with Python bindings, I’ll show you what worked for me..

Because this is merely an adaptation of Matthew Perry’s original solution, I highly recommend considering his original blog entry on the topic, complete with discussion points and caveats, before proceeding further!

Once you’re ready, head over to GitHub and download some code, specifically,, and [1] You don’t actually need, but I snagged all three. To keep things simple put these files in the same folder.

Next, in that same directory, create a new Python file, I called mine

Where works entirely on shapefiles, accepts an OGR PostgreSQL connection string in place of a shapefile path. Fortunately the original OGR code didn’t change, but to use the Mapnik PostGIS driver I had to iterate over the PostgreSQL connection string and store the connection parameters so I could supply them a differently to Mapnik.

Finally, copy the following code and paste it into your new file. It’s basically the same as, but I changed shppath to pgconn and added some code to use a mapnik.PostGIS() datasource in place of a mapnik.Shapefile() datasource.

If you’re curious, the comment # ELR 2014.9.26: flags the few places I changed the original code.

#!/usr/bin/env python
# -*- coding: utf-8  -*-
Author: Matthew Perry
License: BSD

Creates utfgrid .json tiles for the given polygon shapefile

Thx to Dane Springmeyer for the utfgrid spec and mapnik rendering code
and to  Klokan Petr Přidal for his MapTiler code

import globalmaptiles
import mapnik
import ogr
import os
from optparse import OptionParser, OptionError
    import simplejson as json
except ImportError:
    import json

def create_utfgrids(pgconn, minzoom, maxzoom, outdir, fields=None, layernum=0):

    # ELR 2014.9.26:
    # Original implementation pushed in a shapefile path.
    #ds = ogr.Open(shppath)
    ds = ogr.Open(pgconn)

    # ELR 2014.9.26:
    # Iterate over the PostgreSQL connection string and pull out values we need
    # to use Mapnik's PostGIS datasource constructor.
    pgConnARR = pgconn[3:].split(' ')
    for kvPair in pgConnARR:
        if kvPair.split('=')[0] == "host":
            nikHost = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "port":
            nikPort = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "user":
            nikUser = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "password":
            nikPass = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "dbname":
            nikDB = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "tables":
            nikTable = kvPair.split('=')[1]

    print "WARNING:"
    print " This script assumes a polygon shapefile in spherical mercator projection."
    print " If any of these assumptions are not true, don't count on the results!"
    # TODO confirm polygons
    # TODO confirm mercator
    # TODO get layernum from command line
    layer = ds.GetLayer(layernum)
    bbox = layer.GetExtent()
    print ""
    print str(bbox)

    mercator = globalmaptiles.GlobalMercator()

    m = mapnik.Map(256,256)

    # Since grids are `rendered` they need a style
    s = mapnik.Style()
    r = mapnik.Rule()
    polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#f2eff9'))
    line_symbolizer = mapnik.LineSymbolizer(mapnik.Color('rgb(50%,50%,50%)'),0.1)
    m.append_style('My Style',s)

    print ""
    # ELR 2014.9.26:
    # Original implementation using shapefile..
    #ds = mapnik.Shapefile(file=shppath)

    # ELR 2014.9.26:
    # Parameterized PostGIS implementation..
    ds = mapnik.PostGIS(host=nikHost,port=nikPort,user=nikUser,password=nikPass,dbname=nikDB,table=nikTable)

    mlayer = mapnik.Layer('poly')
    mlayer.datasource = ds
    mlayer.styles.append('My Style')

    print ""
    if fields is None:
        fields = mlayer.datasource.fields()
        print "Fields were NONE. Using.."
        print fields
        print "Fields are USER PROVIDED. Using.."
        print fields
    print ""

    for tz in range(minzoom, maxzoom+1):
        print " * Processing Zoom Level %s" % tz
        tminx, tminy = mercator.MetersToTile( bbox[0], bbox[2], tz)
        tmaxx, tmaxy = mercator.MetersToTile( bbox[1], bbox[3], tz)
        for ty in range(tminy, tmaxy+1):
            for tx in range(tminx, tmaxx+1):
                output = os.path.join(outdir, str(tz), str(tx))
                if not os.path.exists(output):

                # Use top origin tile scheme (like OSM or GMaps)
                # TODO support option for TMS bottom origin scheme (ie opt to not invert)
                ymax = 1 << tz;
                invert_ty = ymax - ty - 1;

                tilefilename = os.path.join(output, "%s.json" % invert_ty) # ty for TMS bottom origin
                tilebounds = mercator.TileBounds( tx, ty, tz)
                #print tilefilename, tilebounds

                box = mapnik.Box2d(*tilebounds)
                grid = mapnik.Grid(m.width,m.height)
                utfgrid = grid.encode('utf',resolution=4)
                with open(tilefilename, 'w') as file:

if __name__ == "__main__":
    usage = "usage: %prog [options] shapefile minzoom maxzoom output_directory"
    parser = OptionParser(usage)
    parser.add_option("-f", '--fields', dest="fields", help="Comma-seperated list of fields; default is all")
    (options, args) = parser.parse_args()

    if len(args) != 4:
        parser.error("Incorrect number of arguments")

    pgconn = args[0]
    minzoom, maxzoom = int(args[1]), int(args[2])
    outdir = args[3]

    if os.path.exists(outdir):
        parser.error("output directory exists already")

    if options.fields:
        fields = options.fields.split(",")
        fields = None

    create_utfgrids(pgconn, minzoom, maxzoom, outdir, fields)


Once you’ve prepared, you can call it from the command line like this..

C:\xDev\utfgrids\ "PG:host= port=5432 user=postgres dbname=gis password=passw0rd tables=parcels_pmerc" 12 16 "C:/xGIS/tiles/utf" -f tms,owner_name

  • Hopefully the PostgreSQL connection string ("PG:host=..") makes sense.
  • 12 and 16 represent the minimum and maximum zoom levels to be rendered, respectively.
  • The directory “C:/xGIS/tiles/utf” is where your UTFGrid tiles will be saved.
  • And -f tms,owner_name,the_wkt represents a comma-separated list of data fields you want in your UTFGrid.


  • Both and require your geodata table to be in a Web Mercator projection (EPSG:3857)!
  • The script assumes a top-origin tile scheme, like OSM and others.
  • The script will only work with polygons.
  • While the OGR PostgreSQL connection string has a tables parameter, this implementation will only accept one table.
  • The script will create your target directory, in the example case, utf, and it will throw an error if you create this directory in advance.

[1] Many thanks to Matthew Perry, Klokan Petr Přidal, and Dane Springmeyer for their collective efforts and for sharing their work.

Written by elrobis

September 26th, 2014 at 1:26 pm

OGR VRT: Connect to PostGIS DataSource

without comments

I needed an OGR VRT for something and didn’t find a clear example on the web all in one place, so here goes.

Somewhere on your system, create a new file with a .ovf extension. Inside that file, add some XML like the following to define your PostgreSQL connection:

That name=”WKTGrid” is semantically unrelated here. I have been experimenting with including WKT geometry data in UtfGrid tiles, and that name is relative to my experiments. You can provide most any value for name. However, do note that the layer name is referenced in the ogrinfo command.

  <OGRVRTLayer name="WKTGrid">
    <SrcDataSource>PG:host= user=postgres dbname=gis password=l00per</SrcDataSource>
    <SrcSQL>SELECT tms, owner_name, the_wkt FROM parcels_cama_20140829_pmerc</SrcSQL>

  • OGRVRTLayer The layer name attribute you assign can be anything you want.
  • SrcDataSource The data source value defines your PostgreSQL connection parameters.
  • SrcLayer The source layer identifies the target table in your PostgreSQL instance.
  • SrcSQL [Optional] OGR SQL can be used to target specific fields, define field aliases, and even refine the data set using WHERE clauses, etc.

After you make a VRT, it’s smart to test it in ogrinfo before you use it for anything serious. It’s easy to test a VRT in ogrinfo, and if ogrinfo makes sense of it, then you know you’ve got a good VRT.

A command like this uses ogrinfo and OGR_SQL to open the VRT and isolate one feature, showing you its attributes.

ogrinfo C:\xGIS\Vector\parcels\parcels_cama_20140829_pmerc.ovf -sql " SELECT tms, owner_name, the_wkt FROM WKTGrid WHERE tms = 'R39200-02-31' "

In some cases, OGR may have trouble identifying your geometry field, or you may have multiple geometry fields and want to specify one field in particular. If so, note the following changes, specifically the modification to the SrcSQL node and the added GeometryField node.

  <OGRVRTLayer name="WKTGrid">
    <SrcDataSource>PG:host= user=postgres dbname=gis password=l00per</SrcDataSource>
    <SrcSQL>SELECT ST_AsBinary(wkb_geometry) as geomm, tms, owner_name FROM parcels_cama_20140829_pmerc</SrcSQL>
    <GeometryField encoding="WKB" field="geomm"></GeometryField>

And this is just scratching the surface. Make sure to check out the OGR VRT driver page for a complete list of options available to you.

Written by elrobis

September 24th, 2014 at 4:45 pm

MySQL Implementation of Google’s Encoded Polyline Algorithm

without comments

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!

Written by elrobis

September 22nd, 2014 at 4:26 pm

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

without comments

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();

Written by elrobis

June 6th, 2014 at 10:44 am

PostGREsql/PostGIS Implementation of Google’s Encoded Polyline Algorithm

with 6 comments

[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.
  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.
  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.
 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.
 gEncoded = '';
 ng = ST_NumGeometries(g1);
 g = 1;
     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;
 RETURN gEncoded;
$$ LANGUAGE plpgsql;


 * First of two methods. Pass in a geometry (LINESTRING only).
 * Returns ASCII-encoded point array for use in Google Maps.
  p INT; np INT;
  deltaX INT;
  deltaY INT;
  enX VARCHAR(255);
  enY VARCHAR(255);
  gEncoded TEXT;
  gEncoded = '';
  np = ST_NPoints(g);

  IF np > 3 THEN
    g = ST_SimplifyPreserveTopology(g, 0.00001);
    np = ST_NPoints(g);

  pt1 = ST_SetSRID(ST_MakePoint(0, 0),4326);

    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);
RETURN gEncoded;
$$ LANGUAGE plpgsql;


 * Second of two methods. Accepts a signed integer (LON or LAT
 * by 1e5) and returns an ASCII-encoded coordinate expression.
  e VARCHAR(255);
  s BIT(32);
  b BIT(6);
  n INT;
 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;
 e = e || chr(s::INT+63);

$$ 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]);
       var polygonObject = new google.maps.Polygon(
         paths: polyPaths,
         strokeColor: '#890000',
         strokeOpacity: 1.0,
         strokeWeight: 2


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.

Written by elrobis

January 27th, 2014 at 12:20 pm

PostGIS: query all multipolygon parcels with at least one hole

without comments

I was writing some code to iterate over Well Known Text expressions for polygon features, and I decided I needed to test the most complex edge-case I could think of–multipolygon geometries where at least one of the bound polygons has a hole (i.e. an interior ring).

I ended up with the following query. This seems like the kind of thing I’ll want to reuse later, so I’m noting it here. For good measure, I also use a rudimentary technique to sort the output with the most complicated geometries in the table at the top of the list. Basically, the more “text” it takes to describe the geometry using Well Known Text, the larger and more complex I figure it must be!

  SomePrimaryId,   /* your primary key, i.e. ogc_fid, etc. */
  SomeUniqueId,    /* your descriptive id, i.e. a parcel number */
  ST_NumGeometries(wkb_geometry) AS num_geoms,
  ST_NRings(wkb_geometry) AS num_rings,
  ST_AsText(ST_Centroid(wkb_geometry)) AS center,
  Char_Length(ST_AsText(wkb_geometry)) AS len,
  ST_AsText(wkb_geometry) AS wkt
FROM SomePolygonTable
  ST_NumGeometries(wkb_geometry) > 1
  ST_NRings(wkb_geometry) > ST_NumGeometries(wkb_geometry)
ORDER BY Char_Length(ST_AsText(wkb_geometry)) ASC ;


Just for the sake of promoting caution, I’m not certain this is a definitive approach for identifying the largest geometry in a table, as the length of the binary representation and the length of the readable text representation do not correspond one-to-one. Moreover, a feature could have more vertices that required less precision to express (fewer decimal position), than a geometry with fewer vertices that needed more precision, and then you have to ask, which is bigger, fewer vertices and more text, or more vertices that coincidentally did not require as much text? My conclusion is, the “most complicated geometry” is probably relative to the one asking the question. However for my purposes, this was close enough to put the most complicated stuff at the top of the list.

Written by elrobis

November 8th, 2013 at 10:26 am

Install httplib2 to your Preferred Python Runtime after ArcGIS Unceremoniously Hijacks the First Attempt

without comments

We’re going to use Google Maps Engine for some stuff, and so I thought I’d throw my first codes at it using Python. Well.. Google’s Python example requires the httplib2 library, so I needed to install that. When I did—following the library’s download/instructions—for some reason, the library went into service against ArcGIS’s embedded Python runtime, rather than my favored 2.7 instance, which I use for every single thing.


To fix this, I just ensured the new library’s script was called from my preferred Python runtime rather than let Windows, or some black magic, decide which instance it should integrate with httplib2.

Breaking it down:

1) cd into the unpacked httplib2 directory:

C:\>cd C:\Users\elijah\Downloads\httplib2-0.8

2) Next, you can qualify which Python runtime you want to execute the supplied install script—in my case, this worked to install httplib2 to that specific runtime. This was the exact command I used:

C:\Users\elijah\Downloads\httplib2-0.8>C:\Python27\python.exe install

Voila! Now when I pull up IDLE, relative to my preferred instance of Python 2.7, I can import httplib2 without issues.

It didn’t take too long to figure this out—but it seemed like an easy post should someone else be stumped and on the wrong foot.

Written by elrobis

August 5th, 2013 at 9:42 am

Posted in Python

Tagged with

osm2pgsql help and usage

without comments

I wanted a more convenient place to read the osm2pgsql help output, so this seemed as good a place as any for it. /E

osm2pgsql -h
osm2pgsql SVN version 0.80.0 (32bit id space)
        osm2pgsql [options] planet.osm
        osm2pgsql [options] planet.osm.{gz,bz2}
        osm2pgsql [options] file1.osm file2.osm file3.osm
This will import the data from the OSM file(s) into a PostgreSQL database
suitable for use by the Mapnik renderer
   -a|--append          Add the OSM file into the database without removing
                        existing data.
   -b|--bbox            Apply a bounding box filter on the imported data
                        Must be specified as: minlon,minlat,maxlon,maxlat
                        e.g. --bbox -0.5,51.25,0.5,51.75
   -c|--create          Remove existing data from the database. This is the
                        default if --append is not specified.
   -d|--database        The name of the PostgreSQL database to connect
                        to (default: gis).
   -i|--tablespace-index        The name of the PostgreSQL tablespace where
                        all indexes will be created.
                        The following options allow more fine-grained control:
      --tablespace-main-data    tablespace for main tables
      --tablespace-main-index   tablespace for main table indexes
      --tablespace-slim-data    tablespace for slim mode tables
      --tablespace-slim-index   tablespace for slim mode indexes
                        (if unset, use db's default; -i is equivalent to setting
                        --tablespace-main-index and --tablespace-slim-index)
   -l|--latlong         Store data in degrees of latitude & longitude.
   -m|--merc            Store data in proper spherical mercator (default)
   -M|--oldmerc         Store data in the legacy OSM mercator format
   -E|--proj num        Use projection EPSG:num
   -u|--utf8-sanitize   Repair bad UTF8 input data (present in planet
                        dumps prior to August 2007). Adds about 10% overhead.
   -p|--prefix          Prefix for table names (default planet_osm)
   -s|--slim            Store temporary data in the database. This greatly
                        reduces the RAM usage but is much slower. This switch is
                        required if you want to update with --append later.
                        This program was compiled on a 32bit system, so at most
                        3GB of RAM will be used. If you encounter problems
                        during import, you should try this switch.
      --drop            only with --slim: drop temporary tables after import (no updates).
   -S|--style           Location of the style file. Defaults to /usr/local/share/osm2pgsql/
   -C|--cache           Now required for slim and non-slim modes:
                        Use up to this many MB for caching nodes (default: 800)
   -U|--username        Postgresql user name
                        password can be given by prompt or PGPASS environment variable.
   -W|--password        Force password prompt.
   -H|--host            Database server hostname or socket location.
   -P|--port            Database server port.
   -e|--expire-tiles [min_zoom-]max_zoom        Create a tile expiry list.
   -o|--expire-output filename  Output filename for expired tiles list.
   -r|--input-reader    Input frontend.
                        libxml2   - Parse XML using libxml2. (default)
                        primitive - Primitive XML parsing.
                        pbf       - OSM binary format.
   -O|--output          Output backend.
                        pgsql - Output to a PostGIS database. (default)
                        gazetteer - Output to a PostGIS database suitable for gazetteer
                        null  - No output. Useful for testing.
                        Include attributes for each object in the database.
                        This includes the username, userid, timestamp and version.
                        Note: this option also requires additional entries in your style file.
   -k|--hstore          Add tags without column to an additional hstore (key/value) column to postgresql tables
      --hstore-match-only       Only keep objects that have a value in one of the columns
      -                         (normal action with --hstore is to keep all objects)
   -j|--hstore-all      Add all tags to an additional hstore (key/value) column in postgresql tables
   -z|--hstore-column   Add an additional hstore (key/value) column containing all tags
                        that start with the specified string, eg --hstore-column "name:" will
                        produce an extra hstore column that contains all name:xx tags
   -G|--multi-geometry  Generate multi-geometry features in postgresql tables.
   -K|--keep-coastlines Keep coastline data rather than filtering it out.
                        By default natural=coastline tagged data will be discarded based on the
                        assumption that post-processed Coastline Checker shapefiles will be used.
      --number-processes                Specifies the number of parallel processes used for certain operations
                        Default is 1
   -I|--disable-parallel-indexing       Disable indexing all tables concurrently.
      --unlogged        Use unlogged tables (lost on crash but faster). Requires PostgreSQL 9.1.
      --cache-strategy  Specifies the method used to cache nodes in ram.
                                Available options are:
                                dense: caching strategy optimised for full planet import
                                chunked: caching strategy optimised for non-contigouse memory allocation
                                sparse: caching strategy optimised for small extracts
                                optimized: automatically combines dense and sparse strategies for optimal storage efficiency.
                                           optimized may use twice as much virtual memory, but no more physical memory
                                   The default is "chunked"
   -h|--help            Help information.
   -v|--verbose         Verbose output.
Add -v to display supported projections.
Use -E to access any espg projections (usually in /usr/share/proj/epsg)

Written by elrobis

December 28th, 2012 at 9:21 am

Posted in osm2pgsql

Tagged with

Resident “lunatic” of LinkedIn

without comments

Just updated my LinkedIn profile.. This is just TOO FUNNY right now!

I totally had to save this for posterity. Unfortunately it will probably disappear before you can say “Jack Robinson”. 8)

Written by elrobis

December 21st, 2012 at 11:04 am

Posted in Uncategorized