- Author:
- dhke
- Posted:
- December 1, 2014
- Language:
- Python
- Version:
- 1.7
- Tags:
- Score:
- 0 (after 0 ratings)

GEOSGeometry.buffer() accepts a distance argument which is in the unit of the coordinate reference system of the geometry.

It is often necessary, however, to specify the buffer size in a more down-to-earth coordinate system, for example meters.

with_metric_buffer(geom, buf_size) implements this functionality.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 | ```
# -*- encoding=utf-8 -*-
import math
UTM_NORTHERN_BASE = 32600
UTM_SOUTHERN_BASE = 32700
DEFAULT_MAP_SRID = 3857
def get_utm_zone(point):
"""
Guess the UTM zone ID for the supplied point
"""
# hardcode the target reference system:
# this is WGS84 by spec
wgspoint = point.transform('WGS84', clone=True)
# remember, geographers are sideways
lon = wgspoint.x
return int(1 + (lon + 180.0) / 6.0)
def is_utm_northern(point):
"""
Determine if the supplied point is in the northern or southern part
for the UTM coordinate system.
"""
wgspoint = point.transform('WGS84', clone=True)
lat = wgspoint.y
return lat >= 0.0
def get_utm_srid(point):
"""
Given the input point, guess the UTM zone and hemissphere and
return the SRID for the UTM appropriate UTM zone.
Note that this does not do any range checking, so supplying bogus
points yields undefined results.
"""
utm_zone = get_utm_zone(point)
is_northern = is_utm_northern(point)
if is_northern:
return UTM_NORTHERN_BASE + utm_zone
else:
return UTM_SOUTHERN_BASE + utm_zone
def buffer_from_meters(geom, buffer_meters):
"""
Create a buffer around the supplied geometry
with the specified distance.
This is a wrapper around GEOMSGeometry.buffer()
but with the buffer distance specified in meters.
GEOM should be in the coordinate system the map will be drawn in,
this is usually "web mercator", i.e. EPSG:3857
"""
# The buffer calculation needs to happen in the source
# coordinate system, otherwise stretching occurs
# (e.g. circles become more and more egg-like further away from the equator)
#
# At the same time we want distances in meters to be
# as close to correct as possible given the local environment.
#
# The approach taken here is
# (1) use the centroid of the input geometry to determine a single reference
# point
# (2) transform this reference point into the appropriate UTM coordinate system
# This is done using the correct UTM zone for the reference put.
# to keep differences to a minimum.
# (3) shift the UTM version of the reference point to the north and east
# by buffer_meters / sqrt(2).
# (4) transform the shifted point back from UTM to the input coordinate system
# and calculate the distance between the shifted point and the original point
# in the input coordinate system.
# (5) Use this newly obtain distance value to create
# a buffered geometry from the original.
ref_point = geom.centroid.clone()
if not ref_point.srid:
# default to WGS84
ref_point.srid = 4326
utm_srid = get_utm_srid(ref_point)
utm_point = ref_point.transform(utm_srid, clone=True)
shift_distance = buffer_meters / math.sqrt(2.0)
shifted_point = utm_point.clone()
shifted_point.set_x(shifted_point.get_x() + shift_distance)
shifted_point.set_y(shifted_point.get_y() + shift_distance)
shifted_ref = shifted_point.transform(ref_point.srid, clone=True)
distance = shifted_ref.distance(ref_point)
return distance
def with_metric_buffer(geom, buffer_meters, map_srid=DEFAULT_MAP_SRID):
"""
Create a buffer around geom of the specified
size.
BUFFER_METERS is the size of the geometry buffer in meters.
This is calculated regardless of the input coordinate system.
MAP_SRID is the coordinate system for the map display.
This defaults to DEFAULT_MAP_SRID, i.e. web mercator, EPSG:3857.
The buffer is actually created in this coordinate system
so that geometric figures (circles) are preserved as
much as possible in the display.
For example, creating a buffer around a point object
results in a polygon that will largely resemble a circle
on the map even if it is clearly elliptical in real coordinates.
"""
if not geom.srid:
geom = geom.clone()
geom.srid = 4326
if map_srid:
map_geom = geom.transform(map_srid, clone=True)
else:
map_geom = geom
buf_size = buffer_from_meters(map_geom, buffer_meters)
buffered_geom = map_geom.buffer(buf_size)
buffered_geom.transform(geom.srid)
return buffered_geom
``` |

## Comments

Please login first before commenting.