# 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.

## Comments

### Re: XYZ

Author: Jeremy

Dojo has an implementation of zip/unzip in javascript.

http://servicesbeta.esri.com/dojox.functional/zip.html

Eugene has other goodies on his blog. http://lazutkin.com/blog/2008/jan/12/functional-fun-javascript-dojo/

### Re: XYZ

Author: Grant

I've recently begun using zip for staggered iteration over a list of coordinates, like so:

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)

Stumbled on the idea here: http://bit.ly/gOPACj