XYZ

While answering a fun question on Stack Overflow, I was reminded of the fantastic utility of Python's built in zip when working with coordinates in space. If you're not using zip every chance you get, you are not a true GIS owning Pythonista. What the function does is "zip" sequences together.

>>> zip(['x', 'y', 'z'], ['1', '2', '3'])
[('x', '1'), ('y', '2'), ('z', '3')]

M number of N length sequences become N number of M length tuples. This is handy for geospatial programming because we often go back and forth between sequences of x, y (ignoring z for the now) pairs (vertex semantics) and separate arrays of x and y values (computational efficiency). I'm going to show here how to use zip along with pyproj and Shapely to calculate the surface area of the State of Colorado to a good approximation.

let's say you have a GeoJSON representation of Colorado's boundary:

>>> co = {"type": "Polygon",
...  "coordinates": [[
...    [-102.05, 41.0],
...    [-102.05, 37.0],
...    [-109.05, 37.0],
...    [-109.05, 41.0]
...  ]]}

Note that Colorado is not exactly a rectangle in any coordinate system: the surveyed boundaries exhibit kinks on the borders with Utah and New Mexico. The method I'm going to show isn't invalidated by the boundary anomalies.

First, let's project the geographic coordinates onto a plane using a pyproj equal area projection centered on and bracketing the state.

>>> from pyproj import Proj
>>> projection = Proj(
...     "+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-105.55")

The projection is callable, but it wants unzipped sequences of longitudes and latitudes, not zipped GeoJSON coordinates. Handily, zip also unzips in conjunction with the * operator:

>>> lon, lat = zip(*co['coordinates'][0])
>>> lon
(-102.05, -102.05, -109.05, -109.05)
>>> lat
(41.0, 37.0, 37.0, 41.0)

Next, we project these values and zip them up again to make a projected GeoJSON representation.

>>> x, y = projection(lon, lat)
>>> co_projected = {"type": "Polygon", "coordinates": [zip(x, y)]}

Finally, make a Shapely geometric object from the projected feature and take its area.

>>> from shapely.geometry import shape
>>> shape(co_projected).area
268952044107.43542

This value is 0.3% smaller than the surveyed area of the state.

Use of zip eliminates needless, slower (the built-in is implemented in C) Python code using an idiom that isn't alien to GIS. Once you start to think of coordinate pairs as "zipped" and coordinate arrays as "unzipped", it makes good sense.

The itertools module also provides a iterable version of this function: izip.

import math
coords = [(0,0), (0,1), (2,6), (1,0)]
for cur,next in zip(coords, coords[1:]):
    delta = math.sqrt((next[0]-cur[0])**2 + (next[1]-cur[1])**2)