# -*- 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