Skip to content

geodesics utils, including ellipsoid data and ellipsoidal distance

License

Notifications You must be signed in to change notification settings

lggomez/go-geodesy

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

23 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

geodesy

go-geodesy is a package of geodesy-related utilities made in golang, including ellipsoid constants for development of geographic and geospatial implementations

    import "github.com/lggomez/go-geodesy"

Usage

Point definitions

type Point

type Point [2]float64

Point represents a latitude-longitude pair in decimal degrees

Constants

const (
	LatLowerBound = float64(-90)
	LatUpperBound = float64(90)
	LonLowerBound = float64(-180)
	LonUpperBound = float64(180)
)

func (Point) Antipode

func (p Point) Antipode() Point

Antipode returns a new point representing the geographical antipode of p

func (Point) Equals

func (p Point) Equals(p2 Point) bool

Equals returns whether p is equal in latitude and longitude to p2

func (Point) IsAntipodeOf

func (p Point) IsAntipodeOf(p2 Point) bool

IsAntipodeOf returns whether p is the exact antipode of p2 or not

func (Point) Lat

func (p Point) Lat() float64

Lat returns point p's latitude

func (Point) LatRadians

func (p Point) LatRadians() float64

LatRadians returns point p's latitude in radians

func (Point) Lon

func (p Point) Lon() float64

Lon returns point p's longitude

func (Point) LonRadians

func (p Point) LonRadians() float64

LonRadians returns point p's longitude in radians

Calculating distances

    import "github.com/lggomez/go-geodesy/distance"

func Haversine

func Haversine(p1, p2 geodesy.Point) float64

Haversine calculates the ellipsoidal distance in meters between 2 points using the Haversine formula and the WGS-84 ellipsoid constants. If any of the points does not constitute a valid geographic coordinate, the returned distance will be math.NaN().

func VincentyInverse

func VincentyInverse(p1, p2 geodesy.Point, accuracy float64, calculateAzimuth bool) (float64, float64, float64)

VincentyInverse calculates the ellipsoidal distance in meters and azimuth in degrees between 2 points using the inverse Vincenty formulae and the WGS-84 ellipsoid constants. As it is an iterative operation it will converge to the defined accuracy, if accuracy < 0 it will use the default accuracy of 1e-12 (approximately 0.06 mm, magnitude should be no bigger than 1e-6). If calculateAzimuth is set to true, it will compute the forward and reverse azimuths (otherwise, these default to math.NaN()). If any of the points does not constitute a valid geographic coordinate, the returned distance will be math.NaN().

The following notations are used in the implementation:

* a 	length of semi-major axis of the ellipsoid (radius at equator)
* ƒ 	flattening of the ellipsoid
* b = (1 − ƒ) a 	length of semi-minor axis of the ellipsoid (radius at the poles)
* u1 = arctan( (1 − ƒ) tan lat1 ) 	reduced latitude for p1 (latitude on the auxiliary sphere);
* u2 = arctan( (1 − ƒ) tan lat2 ) 	reduced latitude for p2 (latitude on the auxiliary sphere);
* L1, L2 	longitude of the points;
* L = L2 − L1 	difference in longitude of two points;
* λ 	Difference in longitude of the points on the auxiliary sphere;
* α1, α2 	forward azimuths at the points;
* α 	forward azimuth of the geodesic at the equator, if it were extended that far;
* s 	ellipsoidal distance between the two points;
* σ 	angular separation between points;
* σ1 	angular separation between the point and the equator;
* σm 	angular separation between the midpoint of the line and the equator;

Ellipsoids

     import "github.com/lggomez/go-geodesy/ellipsoids"

GRS-80

const (
	// Geocentric gravitational constant GM, defined in (m^3)/(s^2)
	GRS80_GEOCENTRIC_GRAVITATIONAL_CONSTANT float64 = 3_986_005_000_000_000

	// Dynamical form factor J2; adimensional
	GRS80_DYNAMICAL_FORM_FACTOR float64 = 0.0108263

	// Dynamical form factor ω; defined in s^-1
	GRS80_ANGULAR_VELOCITY float64 = 0.0007292115
)

Defining physical constants

const (
	// Semi minor axis b, defined in meters (m)
	GRS80_SEMI_MINOR_AXIS float64 = 6_356_752.314140

	// Aspect ratio (b/a); adimensional
	GRS80_ASPECT_RATIO float64 = 0.996647189318816362

	// Mean radius R1 = (2a+b)/3, defined in meters (m)
	GRS80_MEAN_RADIUS = 6_371_008.7714
	// Mean radius R2, defined in meters (m)
	GRS80_AUTHALIC_MEAN_RADIUS = 6_371_007.1810
	// Radius of a sphere of the same volume R3 = ((a^2)*b)^(1/3); defined in meters (m)
	GRS80_SPHERE_RADIUS = 6_371_000.7900
	// Polar radius of curvature = (a^2)/b; defined in meters (m)
	GRS80_POLAR_CURVATURE_RADIUS = 6_399_593.6259
	// Equatorial radius of curvature for a meridian = (b^2)/a; defined in meters (m)
	GRS80_MERIDIAN_CURVATURE_EQUATORIAL_RADIUS = 6_335_439.3271

	// Meridian quadrant (meridian quarter); defined in meters (m)
	// See https://en.wikipedia.org/wiki/Meridian_arc#Quarter_meridian
	GRS80_MERIDIAN_QUADRANT = 10_001_965.7293

	// Linear eccentricity c = sqrt((a^2)-(b^2)); defined in meters (m)
	GRS80_LINEAR_ECCENTRICITY = 521_854.0097
	// Eccentricity of elliptical section through poles e = sqrt((a^2)-(b^2))/a; adimensional
	GRS80_LINEAR_ECCENTRICITY_POLES = 0.0818191910435

	// Flattening f; adimensional
	GRS80_FLATTENING float64 = 0.003352810681183637418
	// Flattening inverse (1/f); adimensional
	GRS80_FLATTENING_INVERSE float64 = 1 / 0.003352810681183637418
)

Derived geometrical constants (all rounded)

const (
	// Period of rotation (sidereal day) = 2π/ω; defined in seconds (s)
	GRS80_ROTATION_PERIOD float64 = 8_616.4100637
)

Derived physical constants (all rounded)

const (
	// Semi major axis a, defined in meters (m)
	GRS80_SEMI_MAJOR_AXIS float64 = 6_378_137
)

Derived geometrical constants (all rounded)

WGS-84

const (
	// WGS84_GEOCENTRIC_GRAVITATIONAL_CONSTANT Geocentric gravitational constant GM, defined in (m^3)/(s^2)
	WGS84_GEOCENTRIC_GRAVITATIONAL_CONSTANT float64 = 3_986_005_000_000_000

	// WGS84_DYNAMICAL_FORM_FACTOR Dynamical form factor J2; adimensional
	// See https://ahrs.readthedocs.io/en/latest/wgs84.html#ahrs.utils.wgs84.WGS.dynamical_form_factor
	WGS84_DYNAMICAL_FORM_FACTOR float64 = 0.0010826298213129219

	// WGS84_ANGULAR_VELOCITY Dynamical form factor ω; defined in s^-1
	WGS84_ANGULAR_VELOCITY float64 = 0.0007292115
)

Defining physical constants

const (
	// WGS84_SEMI_MINOR_AXIS Semi minor axis b, defined in meters (m)
	WGS84_SEMI_MINOR_AXIS float64 = 6_356_752.31424518

	// WGS84_ASPECT_RATIO Aspect ratio (b/a); adimensional
	WGS84_ASPECT_RATIO float64 = 0.9966471893352525

	// WGS84_MEAN_RADIUS Mean radius R1 = (2a+b)/3, defined in meters (m)
	WGS84_MEAN_RADIUS = 6_371_008.771415059
	// WGS84_AUTHALIC_MEAN_RADIUS Mean radius R2, defined in meters (m)
	WGS84_AUTHALIC_MEAN_RADIUS = 6_371_007.1809182055
	// WGS84_SPHERE_RADIUS Radius of a sphere of the same volume R3 = ((a^2)*b)^(1/3); defined in meters (m)
	WGS84_SPHERE_RADIUS = 6_371_000.79000916
	// WGS84_POLAR_CURVATURE_RADIUS Polar radius of curvature = (a^2)/b; defined in meters (m)
	WGS84_POLAR_CURVATURE_RADIUS = 6_399_593.625758493
	// WGS84_MERIDIAN_CURVATURE_EQUATORIAL_RADIUS Equatorial radius of curvature for a meridian = (b^2)/a; defined in meters (m)
	WGS84_MERIDIAN_CURVATURE_EQUATORIAL_RADIUS = 6_335_439.327292821

	// WGS84_MERIDIAN_QUADRANT Meridian quadrant (meridian quarter); defined in meters (m)
	// See https://en.wikipedia.org/wiki/Meridian_arc#Quarter_meridian
	WGS84_MERIDIAN_QUADRANT = 10_001_965.729

	// WGS84_LINEAR_ECCENTRICITY Linear eccentricity c = sqrt((a^2)-(b^2)); defined in meters (m)
	WGS84_LINEAR_ECCENTRICITY = 521_854.0084234
	// WGS84_LINEAR_ECCENTRICITY_POLES Eccentricity of elliptical section through poles e = sqrt((a^2)-(b^2))/a; adimensional
	WGS84_LINEAR_ECCENTRICITY_POLES = 0.0818191918426205

	// WGS84_FLATTENING Flattening f; adimensional
	WGS84_FLATTENING float64 = 1 / 298.257223563
	// WGS84_FLATTENING_INVERSE Flattening inverse (1/f); adimensional
	WGS84_FLATTENING_INVERSE float64 = 298.257223563
)

Defining geometrical constants

const (
	// WGS84_ROTATION_PERIOD Period of rotation (sidereal day) = 2π/ω; defined in seconds (s)
	WGS84_ROTATION_PERIOD float64 = 8_616.4100637
)

Derived physical constants (all rounded)

const (
	// WGS84_SEMI_MAJOR_AXIS Semi major axis a, defined in meters (m)
	WGS84_SEMI_MAJOR_AXIS float64 = 6_378_137
)

Defining geometrical constants