Month: October 2012

  • Making weird composite shapes with Shapely

    Let’s say you operate microwave radio links, and you want to see if you could build a link between, say, the huge comms tower at Clear Creek, ON and your tower just south of Aylmer. You know that Erie Shores Wind Farm is in the middle of the route. Do the 77m diameter turbines attenuate your signal to nothing?

    This is a real issue I have to deal with every day, except in reverse: when I’m planning wind farms, I don’t want to be blocking existing licensed links. I use Spectrum Direct to find these links (and precisely two years ago, I showed you how to write Spectrum Direct data as CSV), but up until now I worked out affected zones (called “consultation zones”) by hand.

    The Radio Advisory Board of Canada and the Canadian Wind Energy Association developed a joint protocol on sizing consultation zone for wind turbines, documented in the RABC CANWEA Guidelines. The one for a microwave link comprises:

    • A 1km buffer around the start and end points of the link;
    • A distance-, frequency- and wind-turbine diameter-specific path width between the points.

    Here I’m only going to consider the 2D geometry of the path. Lots of people in the radio industry get paid money to do line-of-sight terrain analysis, and I’m ignoring this for now. All I’m doing here is using Python’s Shapely library to make the dumb-bell shaped buffer around the path and endpoints. I learned of Shapely from Erik Westra’s book Python Geospatial Development. Here’s some code that outputs pipe-separated well-known text that imports nicely into QGIS:

    #!/usr/bin/python
    # -*- coding: utf-8 -*-
    # Calculate RABC/CanWEA Microwave Link Consultation Zone
    # scruss - 2012-10-29
    
    import pyproj
    import shapely.geometry
    from shapely.wkt import dumps
    
    # hard-coded details (sorry)
    startlat = 42.582892  # Clear Creek
    startlong = -80.603875
    endlat = 42.731776  # South Aylmer
    endlong = -80.98768
    frequency = 6e9  # 6 GHz
    rotor = 77.0  # wtg rotor diameter, m
    
    # calculate distance between points
    g = pyproj.Geod(ellps='WGS84')
    (az12, az21, dist) = g.inv(startlong, startlat, endlong, endlat)
    
    # Width of microwave link fresnel zone
    Lc = rotor + 52 * ((dist / 1000) / (frequency / 1e9)) ** 0.5
    
    # calculate line string along path with segments <= 1 km
    lonlats = g.npts(startlong, startlat, endlong, endlat,
                     1 + int(dist / 1000))
    
    # npts doesn't include start/end points, so prepend/append them
    lonlats.insert(0, (startlong, startlat))
    lonlats.append((endlong, endlat))
    
    # translate line string into projected coordinates
    utms = []
    srcp = pyproj.Proj(init='epsg:4326')  # WGS84
    destp = pyproj.Proj(init='epsg:32617')  # WGS84 UTM Zone 17N
    for (lon, lat) in lonlats:
        (utmx, utmy) = pyproj.transform(srcp, destp, lon, lat)
        utms.append((utmx, utmy))
    
    # turn start and end points into Shapely point objects
    startpt = shapely.geometry.Point(utms[0])
    endpt = shapely.geometry.Point(utms[-1])
    
    # draw a 1km buffer around start/end points
    startb = startpt.buffer(1000)
    endb = endpt.buffer(1000)
    
    # convert microwave path to a Shapely LineString
    path = shapely.geometry.LineString(utms)
    pathb = path.buffer(Lc / 2) # buffer it to fresnel zone width
    
    # join the buffers into one shape
    zone = startb.union(pathb).union(endb)
    
    # output the shape as well known text
    print 'id|wkt'
    print '|'.join(('1', dumps(zone)))
    

    That little snippet creates the buffer in projected coordinates from two geographic coordinates as inputs. Grabbing the Erie Shores coordinates from OpenStreetMap and the USGS WMS, let’s see if the beam will pass:

    Uh, no. That’s three turbines in the path in that little section alone. Back to the drawing board.

  • The Invisible City

    Now you see it …

    There’s not a whole lot north of North Bay. Highway 11 winds through some extensive geometry, but few habitations. Until something wonderful (and quite a bit wrong) happened at 46° 31′ 21″ N, 79° 33′ 9″ W on OpenStreetMap. A whole new town about 25 minutes out of North Bay sprung up overnight — and was as quickly deleted as mappers caught and reverted the vandalism.

    Now you don’t.

    The misguided mapper put quite a lot of work into this ephemeral town. It’s got urban rail, parklands, residential areas and more. You can see the detail on this large image of the town (650 KB PNG). If you want to play with the data, here’s a zipped OSM XML file of the area before it was reverted. But please, don’t re-upload it; OpenStreetMap should be for real features on the ground, and not for confusing map data consumers.

    Thanks to Bootprint for noticing this, and to rw__ for reverting the edits.