github.com/cockroachdb/cockroach@v20.2.0-alpha.1+incompatible/pkg/geo/geographiclib/geographiclib.go (about)

     1  // Copyright 2020 The Cockroach Authors.
     2  //
     3  // Use of this software is governed by the Business Source License
     4  // included in the file licenses/BSL.txt.
     5  //
     6  // As of the Change Date specified in that file, in accordance with
     7  // the Business Source License, use of this software will be governed
     8  // by the Apache License, Version 2.0, included in the file
     9  // licenses/APL.txt.
    10  
    11  // Package geographiclib is a wrapper around the GeographicLib library.
    12  package geographiclib
    13  
    14  // #cgo CXXFLAGS: -std=c++14
    15  // #cgo LDFLAGS: -lm
    16  //
    17  // #include "geodesic.h"
    18  // #include "geographiclib.h"
    19  import "C"
    20  
    21  import "github.com/golang/geo/s2"
    22  
    23  var (
    24  	// WGS84Spheroid represents the default WGS84 ellipsoid.
    25  	WGS84Spheroid = MakeSpheroid(6378137, 1/298.257223563)
    26  )
    27  
    28  // Spheroid is an object that can perform geodesic operations
    29  // on a given spheroid.
    30  type Spheroid struct {
    31  	cRepr        C.struct_geod_geodesic
    32  	Radius       float64
    33  	Flattening   float64
    34  	SphereRadius float64
    35  }
    36  
    37  // MakeSpheroid creates a spheroid from a radius and flattening.
    38  func MakeSpheroid(radius float64, flattening float64) Spheroid {
    39  	minorAxis := radius - radius*flattening
    40  	s := Spheroid{
    41  		Radius:       radius,
    42  		Flattening:   flattening,
    43  		SphereRadius: (radius*2 + minorAxis) / 3,
    44  	}
    45  	C.geod_init(&s.cRepr, C.double(radius), C.double(flattening))
    46  	return s
    47  }
    48  
    49  // Inverse solves the geodetic inverse problem on the given spheroid
    50  // (https://en.wikipedia.org/wiki/Geodesy#Geodetic_problems).
    51  // Returns s12 (distance in meters), az1 (azimuth to point 1) and az2 (azimuth to point 2).
    52  func (s *Spheroid) Inverse(a, b s2.LatLng) (s12, az1, az2 float64) {
    53  	var retS12, retAZ1, retAZ2 C.double
    54  	C.geod_inverse(
    55  		&s.cRepr,
    56  		C.double(a.Lat.Degrees()),
    57  		C.double(a.Lng.Degrees()),
    58  		C.double(b.Lat.Degrees()),
    59  		C.double(b.Lng.Degrees()),
    60  		&retS12,
    61  		&retAZ1,
    62  		&retAZ2,
    63  	)
    64  	return float64(retS12), float64(retAZ1), float64(retAZ2)
    65  }
    66  
    67  // InverseBatch computes the sum of the length of the lines represented
    68  // by the line of points.
    69  // This is intended for use for LineStrings. LinearRings/Polygons should use "AreaAndPerimeter".
    70  // Returns the sum of the s12 (distance in meters) units.
    71  func (s *Spheroid) InverseBatch(points []s2.Point) float64 {
    72  	lats := make([]C.double, len(points))
    73  	lngs := make([]C.double, len(points))
    74  	for i, p := range points {
    75  		latlng := s2.LatLngFromPoint(p)
    76  		lats[i] = C.double(latlng.Lat.Degrees())
    77  		lngs[i] = C.double(latlng.Lng.Degrees())
    78  	}
    79  	var result C.double
    80  	C.CR_GEOGRAPHICLIB_InverseBatch(
    81  		&s.cRepr,
    82  		&lats[0],
    83  		&lngs[0],
    84  		C.int(len(points)),
    85  		&result,
    86  	)
    87  	return float64(result)
    88  }
    89  
    90  // AreaAndPerimeter computes the area and perimeter of a polygon on a given spheroid.
    91  // The points must never be duplicated (i.e. do not include the "final" point of a Polygon LinearRing).
    92  // Area is in meter^2, Perimeter is in meters.
    93  func (s *Spheroid) AreaAndPerimeter(points []s2.Point) (area float64, perimeter float64) {
    94  	lats := make([]C.double, len(points))
    95  	lngs := make([]C.double, len(points))
    96  	for i, p := range points {
    97  		latlng := s2.LatLngFromPoint(p)
    98  		lats[i] = C.double(latlng.Lat.Degrees())
    99  		lngs[i] = C.double(latlng.Lng.Degrees())
   100  	}
   101  	var areaDouble, perimeterDouble C.double
   102  	C.geod_polygonarea(
   103  		&s.cRepr,
   104  		&lats[0],
   105  		&lngs[0],
   106  		C.int(len(points)),
   107  		&areaDouble,
   108  		&perimeterDouble,
   109  	)
   110  	return float64(areaDouble), float64(perimeterDouble)
   111  }