Skip to content

Instantly share code, notes, and snippets.

@tmcw

tmcw/Makefile

Created Mar 19, 2013
Embed
What would you like to do?
Making Shadowplay
shadow.shp: ne_110m_land.shp shadowplay.py
python shadowplay.py ne_110m_land.shp shadow.shp
ne_110m_land.zip:
wget 'http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip'
ne_110m_land.shp:
unzip ne_110m_land.zip
touch ne_110m_land.shp
import sys
# shapely and fiona are nice bindings around OGR and GEOS.
# http://macwright.org/2012/10/31/gis-with-python-shapely-fiona.html
from shapely import speedups
from shapely.ops import cascaded_union
from shapely.prepared import prep
from shapely.geometry import LineString, Point, mapping, shape
from fiona import collection
import random
from math import sin
if speedups.available: speedups.enable()
with collection(sys.argv[1], "r") as input:
shapes = []
for f in input:
shapes.append(shape(f['geometry']))
# make one prepared, unified country so that it's simpler to
# do an intersection against
countries = prep(cascaded_union(shapes))
# the output is just linestrings with no properties
schema = { 'geometry': 'LineString', 'properties': {} }
with collection(sys.argv[2], "w", "ESRI Shapefile", schema) as output:
# we do smaller than normal steps - half-degrees - hence the 360 to
# 360 instead of 180 to 180
xs = [x * 0.5 for x in range(-360, 360)]
for y in range(-80, 80):
lsp = []
for x in xs:
p = Point((x, y))
# both lines that intersect countries and those that don't
# are noisy, but this cranks it up if there's an intersection
if countries.contains(p):
lsp.append((x, y + random.random() - 0.5))
else:
lsp.append((x, y + (sin(x * (0.05 + (abs(y) / 80.0) )))))
# construct an actual line from all of those points and write
# it to the file.
ls = LineString(lsp)
output.write({
'properties': { },
'geometry': mapping(ls)
})
print "created feature for %d" % y
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment