Last active
August 29, 2015 14:14
-
-
Save seumasmorrison/24617c7fe6d99465cf89 to your computer and use it in GitHub Desktop.
Based upon Sean Gillies gist - http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# translate_multipoint | |
# Translation of Shapefile record geometries in a functional style. | |
from fiona import collection | |
from functools import partial | |
from itertools import imap | |
import logging | |
log = logging.getLogger() | |
# To begin: a few functions to translate geometry coordinates. They call each | |
# other semi-recursively. GeoJSON was designed to make this possible, BTW. | |
def translate_Point(coords, delta): | |
"""Returns translated Point coordinates, preserving order, number, | |
and dimensionality. | |
delta is a (delta_x, delta_y [, delta_y]) tuple.""" | |
print 'coords', coords | |
print 'delta', delta | |
result = [tuple(c + d for c, d in zip(coords[0], delta))] | |
print 'result', result | |
return result | |
def translate_LineString(coords, delta): | |
"""Returns translated LineString or Ring coordinates, preserving | |
order, number, and dimensionality. | |
delta is a (delta_x, delta_y [, delta_y]) tuple.""" | |
# Calls translate_Point. | |
return list(translate_Point(pt_coords, delta) for pt_coords in coords) | |
def translate(delta, rec): | |
print "translate" | |
"""Returns a record after translating its geometry's coordinates. | |
delta is a (delta_x, delta_y [, delta_y]) tuple.""" | |
# Use lexical dispatch on geometry type to get the proper translation | |
# function. | |
handlers = { | |
'MultiPoint': translate_Point} | |
try: | |
g = rec['geometry'] | |
print 'g', g | |
g['coordinates'] = handlers[g['type']](g['coordinates'], delta ) | |
rec['geometry'] = g | |
return rec | |
except Exception, e: | |
log.exception("Error processing record %s:", rec) | |
# And now the processing code. First, open a source file. | |
with collection("D:\\shawbost_friday_osgb_1936.shp", "r") as source: | |
# Create a sink file for processed features with the same format and | |
# coordinate reference system as the source. | |
with collection( | |
"translated_osgb_1936_dalmore2.shp", | |
"w", | |
driver=source.driver, | |
schema=source.schema, | |
crs=source.crs | |
) as sink: | |
# Example 2D translation, 1 unit eastward and northward. | |
results = imap(partial(translate, (-4153.8348, -2604.8443)), source) | |
sink.writerecords(results) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment