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.
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
y (ignoring z for the now) pairs (vertex semantics) and
separate arrays of
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:
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.
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
Next, we project these values and zip them up again to make a projected GeoJSON representation.
Finally, make a Shapely geometric object from the projected feature and take its area.
This value is 0.3% smaller than the surveyed area of the state.
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
itertools module also provides a iterable version of this function: