Skip to content

Commit

Permalink
Update comments
Browse files Browse the repository at this point in the history
  • Loading branch information
mikeqfu committed Sep 10, 2019
1 parent b0e0bb3 commit d8a9b0d
Showing 1 changed file with 101 additions and 79 deletions.
180 changes: 101 additions & 79 deletions pyhelpers/geom.py
Expand Up @@ -8,6 +8,16 @@

# Convert British national grid (OSGB36) to latitude and longitude (WGS84)
def osgb36_to_wgs84(easting, northing):
"""
:param easting: [numbers.Number]
:param northing: [numbers.Number]
:return: [tuple] (numbers.Number, numbers.Number)
Convert British National grid coordinates (OSGB36 Easting, Northing) to WGS84 latitude and longitude.
Testing e.g.
easting, northing = 530034, 180381 # osgb36_to_wgs84(easting, northing) == (-0.12772404, 51.507407)
"""
import pyproj
osgb36 = pyproj.Proj(init='EPSG:27700') # UK Ordnance Survey, 1936 datum
wgs84 = pyproj.Proj(init='EPSG:4326') # LonLat with WGS84 datum used by GPS units and Google Earth
Expand All @@ -17,24 +27,31 @@ def osgb36_to_wgs84(easting, northing):

# Convert latitude and longitude (WGS84) to British national grid (OSGB36)
def wgs84_to_osgb36(longitude, latitude):
"""
:param longitude: [numbers.Number]
:param latitude: [numbers.Number]
:return: [tuple] ([numbers.Number], [numbers.Number])
Converts coordinates from WGS84 (latitude, longitude) to British National grid (OSGB36) (easting, northing).
Testing e.g.
longitude, latitude = -0.12772404, 51.507407 # wgs84_to_osgb36(longitude, latitude) == (530034, 180381)
"""
import pyproj
wgs84 = pyproj.Proj(init='EPSG:4326') # LonLat with WGS84 datum used by GPS units and Google Earth
osgb36 = pyproj.Proj(init='EPSG:27700') # UK Ordnance Survey, 1936 datum
easting, northing = pyproj.transform(wgs84, osgb36, longitude, latitude)
return easting, northing


# Convert british national grid (OSGB36) to latitude and longitude (WGS84)
# Convert british national grid (OSGB36) to latitude and longitude (WGS84) by calculation
def osgb36_to_wgs84_calc(easting, northing):
"""
:param easting: X
:param northing: Y
:return:
:param easting: [numbers.Number]
:param northing: [numbers.Number]
:return: [tuple] (numbers.Number, numbers.Number)
'easting' and 'northing' are the British national grid coordinates
Convert British National grid coordinates (OSGB36 Easting, Northing) to WGS84 latitude and longitude.
The code below was copied/adapted from Hannah Fry's blog; the original code and post can be found at:
The code below was copied/modified from:
http://www.hannahfry.co.uk/blog/2012/02/01/converting-british-national-grid-to-latitude-and-longitude-ii
"""
# The Airy 180 semi-major and semi-minor axes used for OSGB36 (m)
Expand Down Expand Up @@ -130,17 +147,14 @@ def osgb36_to_wgs84_calc(easting, northing):
return long, lat


# Convert latitude and longitude (WGS84) to british national grid (OSGB36)
# Convert latitude and longitude (WGS84) to British National Grid (OSGB36) by calculation
def wgs84_to_osgb36_calc(latitude, longitude):
"""
:param latitude:
:param longitude:
:return:
:param longitude: [numbers.Number]
:param latitude: [numbers.Number]
:return: [tuple] ([numbers.Number], [numbers.Number])
This function converts lat lon (WGS84) to british national grid (OSBG36)
Convert WGS84 (latitude and longitude) to British National grid coordinates (OSGB36 Easting, Northing).
The code below was copied/adapted from Hannah Fry's blog; the original code and post can be found at:
The code below was copied/modified from:
http://www.hannahfry.co.uk/blog/2012/02/01/converting-latitude-and-longitude-to-british-national-grid
"""
# First convert to radians. These are on the wrong ellipsoid currently: GRS80. (Denoted by _1)
Expand Down Expand Up @@ -228,62 +242,91 @@ def wgs84_to_osgb36_calc(latitude, longitude):
return x, y


# Find the closest point of the given point to a list of points
def find_closest_point(pt, pts):
# Get the midpoint between <shapely.geometry.point.Point>s
def get_geometric_midpoint(pt_x, pt_y, as_geom=True):
"""
:param pt: [tuple] (lon, lat)
:param pts: [iterable] a sequence of reference points
:param pt_x: [shapely.geometry.point.Point]
:param pt_y: [shapely.geometry.point.Point]
:param as_geom: [bool] (default: True)
:return: [tuple]
"""
# Define a function calculating distance between two points
def distance(o, d):
"""
math.hypot(x, y) return the Euclidean norm, sqrt(x*x + y*y).
This is the length of the vector from the origin to point (x, y).
"""
return math.hypot(o[0] - d[0], o[1] - d[1])
# assert isinstance(x_pt, shapely.geometry.point.Point)
# assert isinstance(y_pt, shapely.geometry.point.Point)

# Find the min value using the distance function with coord parameter
return min(pts, key=functools.partial(distance, pt))
midpoint = (pt_x.x + pt_y.x) / 2, (pt_x.y + pt_y.y) / 2

if as_geom:
import shapely.geometry
midpoint = shapely.geometry.Point(midpoint)

return midpoint


# Get the midpoint between two points
def get_geometric_midpoint_along_earth_surface(pt_x, pt_y, as_geom=False):
"""
:param pt_x: [shapely.geometry.point.Point]
:param pt_y: [shapely.geometry.point.Point]
:param as_geom: [bool] (default: False)
:return: [tuple]
References:
http://code.activestate.com/recipes/577713-midpoint-of-two-gps-points/
http://www.movable-type.co.uk/scripts/latlong.html
"""
# Input values as degrees, convert them to radians
lon_1, lat_1 = math.radians(pt_x.x), math.radians(pt_x.y)
lon_2, lat_2 = math.radians(pt_y.x), math.radians(pt_y.y)

b_x, b_y = math.cos(lat_2) * math.cos(lon_2 - lon_1), math.cos(lat_2) * math.sin(lon_2 - lon_1)
lat_3 = math.atan2(
math.sin(lat_1) + math.sin(lat_2), math.sqrt((math.cos(lat_1) + b_x) * (math.cos(lat_1) + b_x) + b_y ** 2))
long_3 = lon_1 + math.atan2(b_y, math.cos(lat_1) + b_x)

midpoint = math.degrees(long_3), math.degrees(lat_3)

if as_geom:
import shapely.geometry
midpoint = shapely.geometry.Point(midpoint)

return midpoint


# Calculate distance between two points
def calc_distance_on_unit_sphere(x_pt, y_pt):
def calc_distance_on_unit_sphere(pt_x, pt_y):
"""
:param x_pt: [shapely.geometry.point.Point]
:param y_pt: [shapely.geometry.point.Point]
:param pt_x: [shapely.geometry.point.Point]
:param pt_y: [shapely.geometry.point.Point]
:return: [float]
The following function returns the distance between two locations based on each point’s longitude and latitude.
The distance returned is relative to Earth’s radius. To get the distance in miles, multiply by 3960. To get the
This function was copied/changed from http://www.johndcook.com/blog/python_longitude_latitude/.
It returns the distance between two locations based on each point's longitude and latitude.
The distance returned is relative to Earth's radius. To get the distance in miles, multiply by 3960. To get the
distance in kilometers, multiply by 6373.
Latitude is measured in degrees north of the equator; southern locations have negative latitude. Similarly,
longitude is measured in degrees east of the Prime Meridian. A location 10° west of the Prime Meridian,
for example, could be expressed as either 350° east or as -10° east.
The function assumes the earth is perfectly spherical. For a discussion of how accurate this assumption is,
see my blog post on http://www.johndcook.com/blog/2009/03/02/what-is-the-shape-of-the-earth/
see http://www.johndcook.com/blog/2009/03/02/what-is-the-shape-of-the-earth/
The algorithm used to calculate distances is described in detail at http://www.johndcook.com/lat_long_details.html
A web page to calculate the distance between to cities based on longitude and latitude is available at
A web page to calculate the distance between two cities based on longitude and latitude is available at
http://www.johndcook.com/lat_long_distance.html
This function is in the public domain. Do whatever you want with it, no strings attached.
Reference: http://www.johndcook.com/blog/python_longitude_latitude/
"""
# Convert latitude and longitude to spherical coordinates in radians.
degrees_to_radians = math.pi / 180.0

# phi = 90 - latitude
phi1 = (90.0 - x_pt.y) * degrees_to_radians
phi2 = (90.0 - y_pt.y) * degrees_to_radians
phi1 = (90.0 - pt_x.y) * degrees_to_radians
phi2 = (90.0 - pt_y.y) * degrees_to_radians

# theta = longitude
theta1 = x_pt.x * degrees_to_radians
theta2 = y_pt.x * degrees_to_radians
theta1 = pt_x.x * degrees_to_radians
theta2 = pt_y.x * degrees_to_radians

# Compute spherical distance from spherical coordinates.

Expand All @@ -295,45 +338,24 @@ def calc_distance_on_unit_sphere(x_pt, y_pt):
cosine = (math.sin(phi1) * math.sin(phi2) * math.cos(theta1 - theta2) + math.cos(phi1) * math.cos(phi2))
arc = math.acos(cosine) * 3960 # in miles

# Remember to multiply arc by the radius of the earth
# in your favorite set of units to get length.
# Remember to multiply arc by the radius of the earth in your favorite set of units to get length.
return arc


# Get the midpoint between two points
def get_midpoint_along_earth_surface(x_pt, y_pt):
"""
:param x_pt: [shapely.geometry.point.Point]
:param y_pt: [shapely.geometry.point.Point]
:return: [tuple]
References:
http://code.activestate.com/recipes/577713-midpoint-of-two-gps-points/
http://www.movable-type.co.uk/scripts/latlong.html
"""
# Input values as degrees, convert them to radians
lon_1, lat_1 = math.radians(x_pt.x), math.radians(x_pt.y)
lon_2, lat_2 = math.radians(y_pt.x), math.radians(y_pt.y)

b_x, b_y = math.cos(lat_2) * math.cos(lon_2 - lon_1), math.cos(lat_2) * math.sin(lon_2 - lon_1)
lat_3 = math.atan2(math.sin(lat_1) + math.sin(lat_2),
math.sqrt((math.cos(lat_1) + b_x) * (math.cos(lat_1) + b_x) + b_y ** 2))
long_3 = lon_1 + math.atan2(b_y, math.cos(lat_1) + b_x)

midpoint = math.degrees(long_3), math.degrees(lat_3)

return midpoint


# Get the midpoint between <shapely.geometry.point.Point>s
def get_geometric_midpoint(x_pt, y_pt):
# Find the closest point of the given point to a list of points
def find_closest_point(pt, pts):
"""
:param x_pt: [shapely.geometry.point.Point]
:param y_pt: [shapely.geometry.point.Point]
:param pt: [tuple] (lon, lat)
:param pts: [iterable] a sequence of reference points
:return: [tuple]
"""
# assert isinstance(x_pt, shapely.geometry.point.Point)
# assert isinstance(y_pt, shapely.geometry.point.Point)
# Define a function calculating distance between two points
def distance(o, d):
"""
math.hypot(x, y) return the Euclidean norm, sqrt(x*x + y*y).
This is the length of the vector from the origin to point (x, y).
"""
return math.hypot(o[0] - d[0], o[1] - d[1])

midpoint = (x_pt.x + y_pt.x) / 2, (x_pt.y + y_pt.y) / 2
return midpoint
# Find the min value using the distance function with coord parameter
return min(pts, key=functools.partial(distance, pt))

0 comments on commit d8a9b0d

Please sign in to comment.