Making Shadowplay

README.md

shadowplay

Requirements:

• Python
• Shapely
• Fiona
• GNU Make or similar

Running it

make

Makefile

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

shadowplay.py

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