Source code for pycoalescence.spatial_algorithms

Simple spatial algorithms required for package functionality.

Algorithms include generation of Voronoi diagrams and spacing points on a landscape using Lloyd's algorithm.
from __future__ import absolute_import

import copy
import logging
import math

    from osgeo import osr, ogr
except ImportError as ie:  # pragma: no cover
    logging.warning("Could not import from osgeo: {}".format(ie))
    from scipy.spatial import Voronoi
except ImportError as ie:  # pragma: no cover
    logging.warning("Cannot import Voronoi from scipy.spatial: {}".format(ie))

[docs]def calculate_distance_between(x1, y1, x2, y2): """ Calculates the distance between the points (x1, y1) and (x2, y2) .. note:: Returns the absolute value :param x1: x coordinate of the first point :param y1: y coordinate of the first point :param x2: x coordinate of the second point :param y2: y coordinate of the second point :return: the absolute distance between the points """ return ((x1 - x2) ** 2 + (y1 - y2) ** 2) ** 0.5
[docs]def calculate_centre_of_mass(points_list): """ Calculates the centre of mass for the non-intersecting polygon defined by points_list. .. note:: the centre of mass will be incorrect for intersecting polygons. .. note:: it is assumed that points_list defines, in order, the vertices of the polygon. The last point is assumed to connect to the first point. :param points_list: a list of x, y points defining the non-intersecting polygon :return: the x,y centre of mass """ centre_x = 0 centre_y = 0 area = 0 size = len(points_list) points_list.append(points_list[0]) for i, each in enumerate(points_list): if i == size: break added_area = each[0] * points_list[i + 1][1] - points_list[i + 1][0] * each[1] area += added_area centre_x += (each[0] + points_list[i + 1][0]) * added_area centre_y += (each[1] + points_list[i + 1][1]) * added_area area *= 0.5 try: centre_x /= 6 * area centre_y /= 6 * area except ZeroDivisionError: raise ValueError("Points not ordered in polygon") return centre_x, centre_y
[docs]def lloyds_algorithm(points_list, maxima, n=7): """ Equally spaces the points in the given landscape defined by (0, x_max), (0, y_max) using Lloyd's algorithm. Algorthim is: - Reflect the points at x=0, x=x_max, y=0 and y=y_max to make boundaries of the Voronoi diagram on the original set of points have finite edges - Define the Voronoi diagram separating the points - Find the centres of the regions of the voronoi diagram for our original set of points - Move the our points to the centres of their voronoi regions - Repeat n times (for convergence) - Edits the points_list to contain the equally-spaced points .. note:: all points are assumed to be in the range x in (0, x_max) and y in (0, y_max) :param points_list: a list of points to be equally spaced in the landscape :param maxima: the maximum size of the landscape to space out within :param n: the number of iterations to perform Lloyd's algorthim for. :return list containing the new point centres. """ # Note that lloyd's algorithm is only tested at quite a high level for num in range(n): vor = Voronoi(reflect_dimensions(points_list, maximums=maxima)) include_regions = [] # Find our set of points for point_index, point in enumerate(vor.points): if 0 <= point[0] < maxima[0] and 0 <= point[1] < maxima[1]: include_regions.append(vor.point_region[point_index]) # Find which regions we want to use # Find the centres of our regions points_list = [] for region_index in include_regions: vertices_to_average = [] for vertex in vor.regions[region_index]: vertices_to_average.append([vor.vertices[vertex][0], vor.vertices[vertex][1]]) points_list.append(calculate_centre_of_mass(vertices_to_average)) return points_list
[docs]def reflect_dimensions(points, maximums): """ Reflects the provided points across x=0, y=0, x=x_max and y=y_max (essentially tiling the polygon 4 times, around the original polygon). :param list points: a list of 2-d points to reflect :param tuple maximums: tuple containing the x and y maximums :return: a list of reflected points """ max_x, max_y = maximums out_points = copy.copy(points) for reflection_y in [-0.0001, max_y + 0.0001]: for point in points: out_points.append([point[0], reflection_y + (reflection_y - point[1])]) for reflection_x in [-0.0001, max_x + 0.0001]: for point in points: out_points.append([reflection_x + (reflection_x - point[0]), point[1]]) return out_points
[docs]def archimedes_spiral(centre_x, centre_y, radius, theta): """ Gets the x, y coordinates on a spiral, given a radius and theta :param int centre_x: the x coordinate of the centre of the spiral :param int centre_y: the y coordinate of the centre of the spiral :param float radius: the distance from the centre of the spiral :param float theta: the angle of rotation :return: tuple of x and y coordinates :rtype: tuple """ return int(math.floor(radius * math.cos(theta) + centre_x)), int(math.floor(radius * math.sin(theta) + centre_y))
[docs]def convert_coordinates(x, y, input_srs, output_srs): """ Converts the coordinates from the input srs to the output srs. :param x: the x coordinate to transform :param y: the y coordinate to transform :param input_srs: the input srs to transform from :param output_srs: the output srs to transform to :rtype: list :return: transformed [x, y] coordinates """ coord_transform = osr.CoordinateTransformation(input_srs, output_srs) # create a geometry from coordinates point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(x, y) point.Transform(coord_transform) return point.GetX(), point.GetY()
[docs]def estimate_sigma_from_distance(distance, n): """ Estimates the sigma value from a rayleigh distribution (2-d normal) from a total distance travelled in n steps. :param float distance: the total distance travelled :param int n: the number of steps :return: an estimation of the sigma value required to generate the distance travelled in n steps """ return distance * (2.0 / (math.pi * n)) ** 0.5