Source code for safe.common.geodesy

"""point.py - Represents a generic point on a sphere as a Python object.

   See documentation of class Point for details.
   Ole Nielsen, ANU 2002
"""


from math import cos, sin, pi
import numpy


[docs]def acos(c): """acos - Safe inverse cosine Input argument c is shrunk to admissible interval to avoid case where a small rounding error causes a math domain error. """ from math import acos as _acos if c > 1: c = 1 if c < -1: c = -1 return _acos(c)
[docs]class Point: """Definition of a generic point on the sphere. Defines a point in terms of latitude and longitude and computes distances to other points on the sphere. Initialise as Point(lat, lon), where lat and lon are in decimal degrees (dd.dddd) Public Methods: distance_to(P) bearing_to(P) dist(P) Author: Ole Nielsen, ANU 2002 """ # class constants R = 6372000 # Approximate radius of Earth (m) degrees2radians = pi / 180.0 def __init__(self, latitude=None, longitude=None): if latitude is None: msg = 'Argument latitude must be specified to Point constructor' raise Exception(msg) if longitude is None: msg = 'Argument longitude must be specified to Point constructor' raise Exception(msg) msg = 'Specified latitude %f was out of bounds' % latitude assert(latitude >= -90 and latitude <= 90.0), msg msg = 'Specified longitude %f was out of bounds' % longitude assert(longitude >= -180 and longitude <= 180.0), msg self.latitude = float(latitude) self.longitude = float(longitude) lat = latitude * self.degrees2radians # Converted to radians lon = longitude * self.degrees2radians # Converted to radians self.coslat = cos(lat) self.coslon = cos(lon) self.sinlat = sin(lat) self.sinlon = sin(lon) #--------------- # Public methods #---------------
[docs] def bearing_to(self, P): """Bearing (in degrees) to point P""" AZ = self.AZ(P) return int(round(AZ / self.degrees2radians))
[docs] def distance_to(self, P): """Distance to point P""" GCA = self.GCA(P) return self.R * GCA
[docs] def approximate_distance_to(self, P): """Very cheap and rough approximation to distance""" return max(abs(self.latitude - P.latitude), abs(self.longitude - P.longitude)) #----------------- # Internal methods #-----------------
def __repr__(self): """Readable representation of point """ d = 2 lat = round(self.latitude, d) lon = round(self.longitude, d) return ' (' + str(lat) + ', ' + str(lon) + ')'
[docs] def GCA(self, P): """Compute the Creat Circle Angle (GCA) between current point and P """ alpha = P.coslon * self.coslon + P.sinlon * self.sinlon # The original formula is alpha = cos(self.lon - P.lon) # but rewriting lets us make us of precomputed trigonometric values. x = alpha * self.coslat * P.coslat + self.sinlat * P.sinlat return acos(x)
[docs] def AZ(self, P): """Compute Azimuth bearing (AZ) from current point to P """ # Compute cosine(AZ), where AZ is the azimuth angle GCA = self.GCA(P) c = P.sinlat - self.sinlat * cos(GCA) c = c / self.coslat / sin(GCA) AZ = acos(c) # Reverse direction if bearing is westward, # i.e. sin(self.lon - P.lon) > 0 # Without this correction the bearing due west, say, will be 90 degrees # because the formulas work in the positive direction which is east. # # Precomputed trigonometric values are used to rewrite the formula: if self.sinlon * P.coslon - self.coslon * P.sinlon > 0: AZ = 2 * pi - AZ return AZ
[docs] def generate_circle(self, radius, resolution=1): """Make a circle about this point. Args: * radius: The desired cirle radius [m] * resolution (optional): Radial distance (degrees) between points on circle. Default is 1 making the circle consist of 360 points Returns: list of lon, lat coordinates defining the circle Note: The circle is defined in geographic coordinates so the distance in meters will be greater than the specified radius in the north south direction. """ # Find first point in circle to the east of the center by bisecting # until desired radius is achieved # Go north until distance is greater than desired step = 0.001 d = 0 while d < radius: stepmin = step step *= 2 p = Point(self.latitude + step, self.longitude) d = self.distance_to(p) # Then bisect until the correct distance is found in degrees stepmax = step while not numpy.allclose(d, radius, rtol=1.0e-6): step = (stepmax + stepmin) / 2 p = Point(self.latitude + step, self.longitude) d = self.distance_to(p) if d > radius: stepmax = step else: stepmin = step #print #print ('Found geographical distance %f degrees yielding radius %f m' # % (step, d)) r = step # Generate circle with geographic radius = step (degrees) P = [[p.longitude, p.latitude]] angle = 0 while angle < 360: theta = angle * self.degrees2radians #print angle, theta, self.bearing_to(p), self.distance_to(p) # Generate new point on circle lat = self.latitude + r * cos(theta) lon = self.longitude + r * sin(theta) p = Point(latitude=lat, longitude=lon) angle += resolution P.append([p.longitude, p.latitude]) # Close polygon P.append(P[0]) return numpy.array(P)