Title:  Spherical Trigonometry 

Description:  Spherical trigonometry for geographic applications. That is, compute distances and related measures for angular (longitude/latitude) locations. 
Authors:  Robert J. Hijmans [cre, aut], Charles Karney [ctb] (GeographicLib), Ed Williams [ctb], Chris Vennes [ctb] 
Maintainer:  Robert J. Hijmans <[email protected]> 
License:  GPL (>=3) 
Version:  1.520 
Built:  20241102 03:12:12 UTC 
Source:  https://github.com/rspatial/geosphere 
This package implements functions that compute various aspects of distance, direction, area, etc. for geographic (geodetic) coordinates. Some of the functions are based on an ellipsoid (spheroid) model of the world, other functions use a (simpler, but less accuarate) spherical model. Functions using an ellipsoid can be recognized by having arguments to specify the ellipsoid's radius and flattening (a
and f
). By setting the value for f
to zero, the ellipsoid becomes a sphere.
There are also functions to compute intersections of of rhumb lines. There are also functions to compute the distance between points and polylines, and to characterize spherical polygons; for random sampling on a sphere, and to compute daylength. See the vignette vignette('geosphere')
for examples.
Geographic locations must be specified in latitude and longitude in degrees (NOT radians). Degrees are (obviously) in decimal notation. Thus 12 degrees, 30 minutes, 10 seconds = 12 + 30/60 + 10/3600 = 12.50278 degrees. The Southern and Western hemispheres have a negative sign.
The default unit of distance is meter; but this can be adjusted by supplying a different radius r
to functions.
Directions are expressed in degrees (North = 0 and 360, East = 90, South = 180, and West = 270 degrees).
David Purdy, Bill Monahan and others for suggestions to improve the package.
Robert Hijmans, using code by C.F.F. Karney and Chris Veness; formulas by Ed Williams; and with contributions from George Wang, Elias Pipping and others. Maintainer: Robert J. Hijmans <[email protected]>
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
https://www.edwilliams.org/avform147.htm
https://www.movabletype.co.uk/scripts/latlong.html
https://en.wikipedia.org/wiki/Great_circle_distance
https://mathworld.wolfram.com/SphericalTrigonometry.html
The "along track distance" is the distance from the start point (p1) to the closest point on the path to a third point (p3),
following a great circle path defined by points p1 and p2. See dist2gc
for the "cross track distance"
alongTrackDistance(p1, p2, p3, r=6378137)
alongTrackDistance(p1, p2, p3, r=6378137)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above 
p3 
as above 
r 
radius of the earth; default = 6378137m 
A distance in units of r (default is meters)
Ed Williams and Robert Hijmans
alongTrackDistance(c(0,0),c(60,60),c(50,40))
alongTrackDistance(c(0,0),c(60,60),c(50,40))
Compute an antipode, or check whether two points are antipodes. Antipodes are places on Earth that are diametrically opposite to one another; and could be connected by a straight line through the centre of the Earth.
Antipodal points are connected by an infinite number of great circles (e.g. the meridians connecting the poles), and can therefore not be used in some great circle based computations.
antipode(p) antipodal(p1, p2, tol=1e9)
antipode(p) antipodal(p1, p2, tol=1e9)
p 
Longitude/latitude of a single point, in degrees; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p1 
as above 
p2 
as above 
tol 
tolerance for equality 
antipodal points or a logical value (TRUE
if antipodal)
Robert Hijmans
https://en.wikipedia.org/wiki/Antipodes
antipode(rbind(c(5,52), c(120,37), c(60,0), c(0,70))) antipodal(c(0,0), c(180,0))
antipode(rbind(c(5,52), c(120,37), c(60,0), c(0,70))) antipodal(c(0,0), c(180,0))
Compute the area of a polygon in angular coordinates (longitude/latitude) on an ellipsoid.
## S4 method for signature 'matrix' areaPolygon(x, a=6378137, f=1/298.257223563, ...) ## S4 method for signature 'SpatialPolygons' areaPolygon(x, a=6378137, f=1/298.257223563, ...)
## S4 method for signature 'matrix' areaPolygon(x, a=6378137, f=1/298.257223563, ...) ## S4 method for signature 'SpatialPolygons' areaPolygon(x, a=6378137, f=1/298.257223563, ...)
x 
longitude/latitude of the points forming a polygon; Must be a matrix or data.frame of 2 columns (first one is longitude, second is latitude) or a SpatialPolygons* object 
a 
major (equatorial) radius of the ellipsoid 
f 
ellipsoid flattening. The default value is for WGS84 
... 
Additional arguments. None implemented 
area in square meters
Use raster::area for polygons that have a planar (projected) coordinate reference system.
This function calls GeographicLib code by C.F.F. Karney
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
p < rbind(c(180,20), c(140,55), c(10, 0), c(140,60), c(180,20)) areaPolygon(p) # Be careful with very large polygons, as they may not be what they seem! # For example, if you wanted a polygon to compute the area equal to about 1/4 of the ellipsoid # this won't work: b < matrix(c(180, 0, 90, 90, 0, 0, 180, 0), ncol=2, byrow=TRUE) areaPolygon(b) # Becausee the shortest path between (180,0) and (0,0) is # over one of the poles, not along the equator! # Inserting a point along the equator fixes that b < matrix(c(180, 0, 0, 0, 90,0, 180, 0), ncol=2, byrow=TRUE) areaPolygon(b)
p < rbind(c(180,20), c(140,55), c(10, 0), c(140,60), c(180,20)) areaPolygon(p) # Be careful with very large polygons, as they may not be what they seem! # For example, if you wanted a polygon to compute the area equal to about 1/4 of the ellipsoid # this won't work: b < matrix(c(180, 0, 90, 90, 0, 0, 180, 0), ncol=2, byrow=TRUE) areaPolygon(b) # Becausee the shortest path between (180,0) and (0,0) is # over one of the poles, not along the equator! # Inserting a point along the equator fixes that b < matrix(c(180, 0, 0, 0, 90,0, 180, 0), ncol=2, byrow=TRUE) areaPolygon(b)
Get the initial bearing (direction; azimuth) to go from point p1
to point p2
(in longitude/latitude) following the shortest path on an ellipsoid (geodetic). Note that the bearing of travel changes continuously while going along the path. A route with constant bearing is a rhumb line (see bearingRhumb
).
bearing(p1, p2, a=6378137, f=1/298.257223563)
bearing(p1, p2, a=6378137, f=1/298.257223563)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above. Can also be missing, in which case the bearing is computed going from the first point to the next and continuing along the following points 
a 
major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
ellipsoid flattening. The default value is for WGS84 
Bearing in degrees
use f=0
to get a bearing on a sphere (great circle)
Robert Hijmans
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
bearing(c(10,10),c(20,20))
bearing(c(10,10),c(20,20))
Bearing (direction of travel; true course) along a rhumb line (loxodrome) between two points.
bearingRhumb(p1, p2)
bearingRhumb(p1, p2)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above 
A direction (bearing) in degrees
Unlike most great circles, a rhumb line is a line of constant bearing (direction), i.e. tracks of constant true course. The meridians and the equator are both rhumb lines and great circles. Rhumb lines approaching a pole become a tightly wound spiral.
Chris Veness and Robert Hijmans, based on formulae by Ed Williams
https://www.edwilliams.org/avform147.htm#Rhumb
https://en.wikipedia.org/wiki/Rhumb_line
bearingRhumb(c(10,10),c(20,20))
bearingRhumb(c(10,10),c(20,20))
Compute the centroid of longitude/latitude polygons. Unlike other functions in this package, there is no spherical trigonometry involved in the implementation of this function. Instead, the function projects the polygon to the (conformal) Mercator coordinate reference system, computes the centroid, and then inversely projects it to longitude and latitude. This approach fails for polygons that include one of the poles (and is rather biased for anything close to the poles). The function should work for polygons that cross the 180/180 meridian (date line).
centroid(x, ...)
centroid(x, ...)
x 
SpatialPolygons* object, or a 2column matrix or data.frame reprenting a single polgyon (longitude/latitude) 
... 
Additional arguments. None implemented 
A matrix (longitude/latitude)
For multipart polygons, the centroid of the largest part is returned.
Robert J. Hijmans
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) centroid(pol)
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) centroid(pol)
Compute daylength (photoperiod) for a latitude and date.
daylength(lat, doy)
daylength(lat, doy)
lat 
latitude, in degrees. I.e. between 90.0 and 90.0 
doy 
integer, day of the year (1..365) for common (nonleap) years; or an object of class Date; or a character that can be coerced into a date, using 'yyyymmdd' format, e.g. '19821123' 
Daylength in hours
Robert J. Hijmans
Forsythe, William C., Edward J. Rykiel Jr., Randal S. Stahl, Hsini Wu and Robert M. Schoolfield, 1995. A model comparison for daylength as a function of latitude and day of the year. Ecological Modeling 80:8795.
daylength(25, '20101010') daylength(45, 1:365) # average monthly daylength dl < daylength(45, 1:365) tapply(dl, rep(1:12, c(31,28,31,30,31,30,31,31,30,31,30,31)), mean)
daylength(25, '20101010') daylength(45, 1:365) # average monthly daylength dl < daylength(45, 1:365) tapply(dl, rep(1:12, c(31,28,31,30,31,30,31,31,30,31,30,31)), mean)
Given a start point, initial bearing (direction), and distance, this function computes the destination point travelling along a the shortest path on an ellipsoid (the geodesic).
destPoint(p, b, d, a=6378137, f=1/298.257223563, ...)
destPoint(p, b, d, a=6378137, f=1/298.257223563, ...)
p 
Longitude and Latitude of point(s), in degrees. Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
b 
numeric. Bearing (direction) in degrees 
d 
numeric. Distance in meters 
a 
major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
ellipsoid flattening. The default value is for WGS84 
... 
additional arguments. If an argument 'r' is supplied, this is taken as the radius of the earth (e.g. 6378137 m) and computations are for a sphere (great circle) instead of an ellipsoid (geodetic). This is for backwards compatibility only 
A pair of coordinates (longitude/latitude)
Direction changes continuously when travelling along a geodesic. Therefore, the final direction is not the same as the initial direction. You can compute the final direction with finalBearing
(see examples, below)
This function calls GeographicLib code by C.F.F. Karney
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
p < cbind(5,52) d < destPoint(p,30,10000) d #final direction, when arriving at endpoint: finalBearing(d, p)
p < cbind(5,52) d < destPoint(p,30,10000) d #final direction, when arriving at endpoint: finalBearing(d, p)
Calculate the destination point when travelling along a 'rhumb line' (loxodrome), given a start point, direction, and distance.
destPointRhumb(p, b, d, r = 6378137)
destPointRhumb(p, b, d, r = 6378137)
p 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
b 
bearing (direction) in degrees 
d 
distance; in the same unit as 
r 
radius of the earth; default = 6378137 m 
Coordinates (longitude/latitude) of a point
Chris Veness; ported to R by Robert Hijmans
https://www.edwilliams.org/avform147.htm#Rhumb
https://www.movabletype.co.uk/scripts/latlong.html
https://en.wikipedia.org/wiki/Rhumb_line
destPointRhumb(c(0,0), 30, 100000, r = 6378137)
destPointRhumb(c(0,0), 30, 100000, r = 6378137)
Compute the distance of a point to a greatcircle path (also referred to as the cross track distance or cross track error). The great circle is defined by p1
and p2
, while p3
is the point away from the path.
dist2gc(p1, p2, p3, r=6378137, sign=FALSE)
dist2gc(p1, p2, p3, r=6378137, sign=FALSE)
p1 
Start of great circle path. longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
End of great circle path. As above 
p3 
Point away from the great cricle path. As for p2 
r 
radius of the earth; default = 6378137 
sign 
logical. If 
A distance in units of r
(default is meters)
If sign=TRUE
, the sign indicates which side of the path p3 is on. Positive means right of the course from p1 to p2, negative means left.
Ed Williams and Robert Hijmans
https://www.movabletype.co.uk/scripts/latlong.html
https://www.edwilliams.org/ftp/avsig/avform.txt
dist2gc(c(0,0),c(90,90),c(80,80))
dist2gc(c(0,0),c(90,90),c(80,80))
The shortest distance between points and polylines or polygons.
dist2Line(p, line, distfun=distGeo)
dist2Line(p, line, distfun=distGeo)
p 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
line 
longitude/latitude of line as a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialLines* or SpatialPolygons* object 
distfun 
A distance function, such as distGeo 
matrix with distance and lon/lat of the nearest point on the line. Distance is in the same unit as r
in the distfun
(default is meters). If line
is a Spatial*
object, the ID (index) of (one of) the nearest objects is also returned. Thus if the objects are polygons and the point is inside a polygon the function may return the ID of a neighboring polygon that shares the nearest border. You can use the intersect
function in packages terra
.
George Wang and Robert Hijmans
line < rbind(c(180,20), c(150,10), c(140,55), c(10, 0), c(140,60)) pnts < rbind(c(170,0), c(75,0), c(70,10), c(80,20), c(100,50), c(100,60), c(100,40), c(100,20), c(100,10), c(100,0)) d = dist2Line(pnts, line) plot( makeLine(line), type='l') points(line) points(pnts, col='blue', pch=20) points(d[,2], d[,3], col='red', pch='x') for (i in 1:nrow(d)) lines(gcIntermediate(pnts[i,], d[i,2:3], 10), lwd=2)
line < rbind(c(180,20), c(150,10), c(140,55), c(10, 0), c(140,60)) pnts < rbind(c(170,0), c(75,0), c(70,10), c(80,20), c(100,50), c(100,60), c(100,40), c(100,20), c(100,10), c(100,0)) d = dist2Line(pnts, line) plot( makeLine(line), type='l') points(line) points(pnts, col='blue', pch=20) points(d[,2], d[,3], col='red', pch='x') for (i in 1:nrow(d)) lines(gcIntermediate(pnts[i,], d[i,2:3], 10), lwd=2)
The shortest distance between two points (i.e., the 'greatcircledistance' or 'as the crow flies'), according to the 'law of the cosines'. This method assumes a spherical earth, ignoring ellipsoidal effects.
distCosine(p1, p2, r=6378137)
distCosine(p1, p2, r=6378137)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above 
r 
radius of the earth; default = 6378137 m 
Vector of distances in the same unit as r
(default is meters)
Robert Hijmans
https://en.wikipedia.org/wiki/Great_circle_distance
distGeo, distHaversine, distVincentySphere, distVincentyEllipsoid, distMeeus
distCosine(c(0,0),c(90,90))
distCosine(c(0,0),c(90,90))
Highly accurate estimate of the shortest distance between two points on an ellipsoid (default is WGS84 ellipsoid). The shortest path between two points on an ellipsoid is called the geodesic.
distGeo(p1, p2, a=6378137, f=1/298.257223563)
distGeo(p1, p2, a=6378137, f=1/298.257223563)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first column is longitude, second column is latitude) or a SpatialPoints* object 
p2 
as above; or missing, in which case the sequential distance between the points in p1 is computed 
a 
numeric. Major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
numeric. Ellipsoid flattening. The default value is for WGS84 
Parameters from the WGS84 ellipsoid are used by default. It is the best available global ellipsoid, but for some areas other ellipsoids could be preferable, or even necessary if you work with a printed map that refers to that ellipsoid. Here are parameters for some commonly used ellipsoids. Also see the refEllipsoids
function.
ellipsoid 
a 
f 

WGS84 
6378137 
1/298.257223563 

GRS80 
6378137 
1/298.257222101 

GRS67 
6378160 
1/298.25 

Airy 1830 
6377563.396 
1/299.3249646 

Bessel 1841 
6377397.155 
1/299.1528434 

Clarke 1880 
6378249.145 
1/293.465 

Clarke 1866 
6378206.4 
1/294.9786982 

International 1924 
6378388 
1/297 

Krasovsky 1940 
6378245 
1/298.2997381 

more info: https://en.wikipedia.org/wiki/Reference_ellipsoid
Vector of distances in meters
This function calls GeographicLib code by C.F.F. Karney
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
distCosine, distHaversine, distVincentySphere, distVincentyEllipsoid, distMeeus
distGeo(c(0,0),c(90,90))
distGeo(c(0,0),c(90,90))
The shortest distance between two points (i.e., the 'greatcircledistance' or 'as the crow flies'), according to the 'haversine method'. This method assumes a spherical earth, ignoring ellipsoidal effects.
distHaversine(p1, p2, r=6378137)
distHaversine(p1, p2, r=6378137)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above; or missing, in which case the sequential distance between the points in p1 is computed 
r 
radius of the earth; default = 6378137 m 
The Haversine ('halfversedsine') formula was published by R.W. Sinnott in 1984, although it has been known for much longer. At that time computational precision was lower than today (15 digits precision). With current precision, the spherical law of cosines formula appears to give equally good results down to very small distances. If you want greater accuracy, you could use the distVincentyEllipsoid
method.
Vector of distances in the same unit as r
(default is meters)
Chris Veness and Robert Hijmans
Sinnott, R.W, 1984. Virtues of the Haversine. Sky and Telescope 68(2): 159
https://www.movabletype.co.uk/scripts/latlong.html
https://en.wikipedia.org/wiki/Great_circle_distance
distGeo, distCosine, distVincentySphere, distVincentyEllipsoid, distMeeus
distHaversine(c(0,0),c(90,90))
distHaversine(c(0,0),c(90,90))
Distance matrix of a set of points, or between two sets of points
distm(x, y, fun=distGeo)
distm(x, y, fun=distGeo)
x 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
y 
Same as 
fun 
A function to compute distances (e.g., distCosine or distGeo) 
Matrix of distances
Robert Hijmans
https://en.wikipedia.org/wiki/Great_circle_distance
distGeo, distCosine, distHaversine, distVincentySphere, distVincentyEllipsoid
xy < rbind(c(0,0),c(90,90),c(10,10),c(120,45)) distm(xy) xy2 < rbind(c(0,0),c(10,10)) distm(xy, xy2)
xy < rbind(c(0,0),c(90,90),c(10,10),c(120,45)) distm(xy) xy2 < rbind(c(0,0),c(10,10)) distm(xy, xy2)
The shortest distance between two points on an ellipsoid (the 'geodetic'), according to the 'Meeus' method. distGeo
should be more accurate.
distMeeus(p1, p2, a=6378137, f=1/298.257223563)
distMeeus(p1, p2, a=6378137, f=1/298.257223563)
p1 
longitude/latitude of point(s), in degrees 1; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above; or missing, in which case the sequential distance between the points in p1 is computed 
a 
numeric. Major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
numeric. Ellipsoid flattening. The default value is for WGS84 
Parameters from the WGS84 ellipsoid are used by default. It is the best available global ellipsoid, but for some areas other ellipsoids could be preferable, or even necessary if you work with a printed map that refers to that ellipsoid. Here are parameters for some commonly used ellipsoids:
ellipsoid 
a 
f 

WGS84 
6378137 
1/298.257223563 

GRS80 
6378137 
1/298.257222101 

GRS67 
6378160 
1/298.25 

Airy 1830 
6377563.396 
1/299.3249646 

Bessel 1841 
6377397.155 
1/299.1528434 

Clarke 1880 
6378249.145 
1/293.465 

Clarke 1866 
6378206.4 
1/294.9786982 

International 1924 
6378388 
1/297 

Krasovsky 1940 
6378245 
1/298.2997381 

more info: https://en.wikipedia.org/wiki/Reference_ellipsoid
Distance value in the same units as parameter a
of the ellipsoid (default is meters)
This algorithm is also used in the spDists
function in the sp package
Robert Hijmans, based on a script by Stephen R. Schmitt
Meeus, J., 1999 (2nd edition). Astronomical algoritms. WillmanBell, 477p.
distGeo, distVincentyEllipsoid, distVincentySphere, distHaversine, distCosine
distMeeus(c(0,0),c(90,90)) # on a 'Clarke 1880' ellipsoid distMeeus(c(0,0),c(90,90), a=6378249.145, f=1/293.465)
distMeeus(c(0,0),c(90,90)) # on a 'Clarke 1880' ellipsoid distMeeus(c(0,0),c(90,90), a=6378249.145, f=1/293.465)
A rhumb line (loxodrome) is a path of constant bearing (direction), which crosses all meridians at the same angle.
distRhumb(p1, p2, r=6378137)
distRhumb(p1, p2, r=6378137)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above; or missing, in which case the sequential distance between the points in p1 is computed 
r 
radius of the earth; default = 6378137 m 
Rhumb (from the Spanish word for course, 'rumbo') lines are straight lines on a Mercator projection map. They were used in navigation because it is easier to follow a constant compass bearing than to continually adjust the bearing as is needed to follow a great circle, even though rhumb lines are normally longer than greatcircle (orthodrome) routes. Most rhumb lines will gradually spiral towards one of the poles.
distance in units of r (default=meters)
Robert Hijmans and Chris Veness
https://www.movabletype.co.uk/scripts/latlong.html
distCosine, distHaversine, distVincentySphere, distVincentyEllipsoid
distRhumb(c(10,10),c(20,20))
distRhumb(c(10,10),c(20,20))
The shortest distance between two points (i.e., the 'greatcircledistance' or 'as the crow flies'), according to the 'Vincenty (ellipsoid)' method. This method uses an ellipsoid and the results are very accurate. The method is computationally more intensive than the other greatcircled methods in this package.
distVincentyEllipsoid(p1, p2, a=6378137, b=6356752.3142, f=1/298.257223563)
distVincentyEllipsoid(p1, p2, a=6378137, b=6356752.3142, f=1/298.257223563)
p1 
longitude/latitude of point(s), in degrees 1; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above; or missing, in which case the sequential distance between the points in p1 is computed 
a 
Equatorial axis of ellipsoid 
b 
Polar axis of ellipsoid 
f 
Inverse flattening of ellipsoid 
The WGS84 ellipsoid is used by default. It is the best available global ellipsoid, but for some areas other ellipsoids could be preferable, or even necessary if you work with a printed map that refers to that ellipsoid. Here are parameters for some commonly used ellipsoids:
ellipsoid 
a 
b 
f 

WGS84 
6378137 
6356752.3142 
1/298.257223563 

GRS80 
6378137 
6356752.3141 
1/298.257222101 

GRS67 
6378160 
6356774.719 
1/298.25 

Airy 1830 
6377563.396 
6356256.909 
1/299.3249646 

Bessel 1841 
6377397.155 
6356078.965 
1/299.1528434 

Clarke 1880 
6378249.145 
6356514.86955 
1/293.465 

Clarke 1866 
6378206.4 
6356583.8 
1/294.9786982 

International 1924 
6378388 
6356911.946 
1/297 

Krasovsky 1940 
6378245 
6356863 
1/298.2997381 

a
is the 'semimajor axis', and b
is the 'semiminor axis' of the ellipsoid. f
is the flattening.
Note that f = (ab)/a
more info: https://en.wikipedia.org/wiki/Reference_ellipsoid
Distance value in the same units as the ellipsoid (default is meters)
Chris Veness and Robert Hijmans
Vincenty, T. 1975. Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations. Survey Review Vol. 23, No. 176, pp8893. Available here:
https://www.movabletype.co.uk/scripts/latlongvincenty.html
https://en.wikipedia.org/wiki/Great_circle_distance
distGeo, distVincentySphere, distHaversine, distCosine, distMeeus
distVincentyEllipsoid(c(0,0),c(90,90)) # on a 'Clarke 1880' ellipsoid distVincentyEllipsoid(c(0,0),c(90,90), a=6378249.145, b=6356514.86955, f=1/293.465)
distVincentyEllipsoid(c(0,0),c(90,90)) # on a 'Clarke 1880' ellipsoid distVincentyEllipsoid(c(0,0),c(90,90), a=6378249.145, b=6356514.86955, f=1/293.465)
The shortest distance between two points (i.e., the 'greatcircledistance' or 'as the crow flies'), according to the 'Vincenty (sphere)' method.
This method assumes a spherical earth, ignoring ellipsoidal effects and it is less accurate then the distVicentyEllipsoid
method.
distVincentySphere(p1, p2, r=6378137)
distVincentySphere(p1, p2, r=6378137)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above; or missing, in which case the sequential distance between the points in p1 is computed 
r 
radius of the earth; default = 6378137 m 
Distance value in the same unit as r
(default is meters)
Robert Hijmans
https://en.wikipedia.org/wiki/Great_circle_distance
distGeo, distVincentyEllipsoid, distHaversine, distCosine, distMeeus
distVincentySphere(c(0,0),c(90,90))
distVincentySphere(c(0,0),c(90,90))
Get the final direction (bearing) when arriving at p2
after starting from p1
and following the shortest path on an ellipsoid (following a geodetic) or on a sphere (following a great circle).
finalBearing(p1, p2, a=6378137, f=1/298.257223563, sphere=FALSE)
finalBearing(p1, p2, a=6378137, f=1/298.257223563, sphere=FALSE)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first column is longitude, second column is latitude) or a SpatialPoints* object 
p2 
as above 
a 
major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
ellipsoid flattening. The default value is for WGS84 
sphere 
logical. If 
A vector of directions (bearings) in degrees
This function calls GeographicLib code by C.F.F. Karney
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
bearing(c(10,10),c(20,20)) finalBearing(c(10,10),c(20,20))
bearing(c(10,10),c(20,20)) finalBearing(c(10,10),c(20,20))
Get the two points where two great cricles cross each other. Great circles are defined by two points on it.
gcIntersect(p1, p2, p3, p4)
gcIntersect(p1, p2, p3, p4)
p1 
Longitude/latitude of a single point, in degrees; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
As above 
p3 
As above 
p4 
As above 
two points for each pair of great circles
Robert Hijmans, based on equations by Ed Williams (see reference)
https://www.edwilliams.org/intersect.htm
p1 < c(5,52); p2 < c(120,37); p3 < c(60,0); p4 < c(0,70) gcIntersect(p1,p2,p3,p4)
p1 < c(5,52); p2 < c(120,37); p3 < c(60,0); p4 < c(0,70) gcIntersect(p1,p2,p3,p4)
Get the two points where two great cricles cross each other. In this function, great circles are defined by a points and an initial bearing. In function gcIntersect
they are defined by two sets of points.
gcIntersectBearing(p1, brng1, p2, brng2)
gcIntersectBearing(p1, brng1, p2, brng2)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
brng1 
Bearing from p1 
p2 
As above. Should have same length as p1, or a single point (or vice versa when p1 is a single point 
brng2 
Bearing from p2 
a matrix with four columns (two points)
Chris Veness and Robert Hijmans based on code by Ed Williams
https://www.edwilliams.org/avform147.htm#Intersection
https://www.movabletype.co.uk/scripts/latlong.html
gcIntersectBearing(c(10,0), 10, c(10,0), 10)
gcIntersectBearing(c(10,0), 10, c(10,0), 10)
Latitude at which a great circle crosses a longitude
gcLat(p1, p2, lon)
gcLat(p1, p2, lon)
p1 
Longitude/latitude of a single point, in degrees; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
As above 
lon 
Longitude 
A numeric (latitude)
Robert Hijmans based on a formula by Ed Williams
https://www.edwilliams.org/avform147.htm#Int
gcLat(c(5,52), c(120,37), lon=120)
gcLat(c(5,52), c(120,37), lon=120)
Longitudes at which a great circle crosses a latitude (parallel)
gcLon(p1, p2, lat)
gcLon(p1, p2, lat)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above 
lat 
a latitude 
vector of two numbers (longitudes)
Robert Hijmans based on code by Ed Williams
https://www.edwilliams.org/avform147.htm#Intersection
gcLon(c(5,52), c(120,37), 40)
gcLon(c(5,52), c(120,37), 40)
What is northern most point that will be reached when following a great circle? Computed with Clairaut's formula. The southern most point is the antipode
of the northernmost point. This does not seem to be very precise; and you could use optimization instead to find this point (see examples)
gcMaxLat(p1, p2)
gcMaxLat(p1, p2)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above 
A matrix with coordinates (longitude/latitude)
Ed Williams, Chris Veness, Robert Hijmans
https://www.edwilliams.org/ftp/avsig/avform.txt
https://www.movabletype.co.uk/scripts/latlong.html
gcMaxLat(c(5,52), c(120,37)) # Another way to get there: f < function(lon){gcLat(c(5,52), c(120,37), lon)} optimize(f, interval=c(180, 180), maximum=TRUE)
gcMaxLat(c(5,52), c(120,37)) # Another way to get there: f < function(lon){gcLat(c(5,52), c(120,37), lon)} optimize(f, interval=c(180, 180), maximum=TRUE)
Highly accurate estimate of the 'geodesic problem' (find location and azimuth at arrival when departing from a location, given an direction (azimuth) at departure and distance) and the 'inverse geodesic problem' (find the distance between two points and the azimuth of departure and arrival for the shortest path. Computations are for an ellipsoid (default is WGS84 ellipsoid).
This is a direct implementation of the the GeographicLib code by C.F.F. Karney that is also used in several other functions in this package (for example, in distGeo
and areaPolygon
).
geodesic(p, azi, d, a=6378137, f=1/298.257223563, ...) geodesic_inverse(p1, p2, a=6378137, f=1/298.257223563, ...)
geodesic(p, azi, d, a=6378137, f=1/298.257223563, ...) geodesic_inverse(p1, p2, a=6378137, f=1/298.257223563, ...)
p 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first column is longitude, second column is latitude) or a SpatialPoints* object 
p1 
as above 
p2 
as above 
azi 
numeric. Azimuth of departure in degrees 
d 
numeric. Distance in meters 
a 
numeric. Major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
numeric. Ellipsoid flattening. The default value is for WGS84 
... 
additional arguments (none implemented) 
Parameters from the WGS84 ellipsoid are used by default. It is the best available global ellipsoid, but for some areas other ellipsoids could be preferable, or even necessary if you work with a printed map that refers to that ellipsoid. Here are parameters for some commonly used ellipsoids.
ellipsoid 
a 
f 

WGS84 
6378137 
1/298.257223563 

GRS80 
6378137 
1/298.257222101 

GRS67 
6378160 
1/298.25 

Airy 1830 
6377563.396 
1/299.3249646 

Bessel 1841 
6377397.155 
1/299.1528434 

Clarke 1880 
6378249.145 
1/293.465 

Clarke 1866 
6378206.4 
1/294.9786982 

International 1924 
6378388 
1/297 

Krasovsky 1940 
6378245 
1/298.2997381 

more info: https://en.wikipedia.org/wiki/Reference_ellipsoid
Three column matrix with columns 'longitude', 'latitude', 'azimuth' (geodesic); or 'distance' (in meters), 'azimuth1' (of departure), 'azimuth2' (of arrival) (geodesic_inverse)
This function calls GeographicLib code by C.F.F. Karney
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
geodesic(cbind(0,0), 30, 1000000) geodesic_inverse(cbind(0,0), cbind(90,90))
geodesic(cbind(0,0), 30, 1000000) geodesic_inverse(cbind(0,0), cbind(90,90))
mean location for spherical (longitude/latitude) coordinates that deals with the angularity. I.e., the mean of longitudes 179 and 178 is 179.5
geomean(xy, w)
geomean(xy, w)
xy 
matrix with two columns (longitude/latitude), or a SpatialPoints or SpatialPolygons object with a longitude/latitude CRS 
w 
weights (vector of numeric values, with a length that is equal to the number of spatial features in 
Ccoordinate pair (numeric)
Robert J. Hijmans
xy < cbind(x=c(179,179, 177), y=c(12,14,16)) xy geomean(xy)
xy < cbind(x=c(179,179, 177), y=c(12,14,16)) xy geomean(xy)
Get points on a great circle as defined by the shortest distance between two specified points
greatCircle(p1, p2, n=360, sp=FALSE)
greatCircle(p1, p2, n=360, sp=FALSE)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above 
n 
The requested number of points on the Great Circle 
sp 
Logical. Return a SpatialLines object? 
A matrix of points, or a list of such matrices (e.g., if multiple bearings are supplied)
Robert Hijmans, based on a formula provided by Ed Williams
https://www.edwilliams.org/avform147.htm#Int
greatCircle(c(5,52), c(120,37), n=36)
greatCircle(c(5,52), c(120,37), n=36)
Get points on a great circle as defined by a point and an initial bearing
greatCircleBearing(p, brng, n=360)
greatCircleBearing(p, brng, n=360)
p 
longitude/latitude of a single point. Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
brng 
bearing 
n 
The requested number of points on the great circle 
A matrix of points, or a list of matrices (e.g., if multiple bearings are supplied)
Robert Hijmans based on formulae by Ed Williams
https://www.edwilliams.org/avform147.htm#Int
greatCircleBearing(c(5,52), 45, n=12)
greatCircleBearing(c(5,52), 45, n=12)
Empirical function to compute the distance to the horizon from a given altitude. The earth is assumed to be smooth, i.e. mountains and other obstacles are ignored.
horizon(h, r=6378137)
horizon(h, r=6378137)
h 
altitude, numeric >= 0. Should have the same unit as r 
r 
radius of the earth; default value is 6378137 m 
Distance in units of h
(default is meters)
Robert J. Hijmans
https://www.edwilliams.org/avform147.htm#Horizon
Bowditch, 1995. American Practical Navigator. Table 12.
horizon(1.80) # me horizon(324) # Eiffel tower
horizon(1.80) # me horizon(324) # Eiffel tower
Get intermediate points (way points) between the two locations with longitude/latitude coordinates. gcIntermediate is based on a spherical model of the earth and internally uses distCosine
.
gcIntermediate(p1, p2, n=50, breakAtDateLine=FALSE, addStartEnd=FALSE, sp=FALSE, sepNA)
gcIntermediate(p1, p2, n=50, breakAtDateLine=FALSE, addStartEnd=FALSE, sp=FALSE, sepNA)
p1 
longitude/latitude of a single point, in degrees. This can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as for 
n 
integer. The desired number of intermediate points 
breakAtDateLine 
logical. Return two matrices if the dateline is crossed? 
addStartEnd 
logical. Add p1 and p2 to the result? 
sp 
logical. Return a SpatialLines object? 
sepNA 
logical. Rather than as a list, return the values as a two column matrix with lines seperated by a row of NA values? (for use in 'plot') 
matrix or list with intermediate points
Robert Hijmans based on code by Ed Williams (great circle)
https://www.edwilliams.org/avform147.htm#Intermediate
gcIntermediate(c(5,52), c(120,37), n=6, addStartEnd=TRUE)
gcIntermediate(c(5,52), c(120,37), n=6, addStartEnd=TRUE)
Compute the length of lines
lengthLine(line)
lengthLine(line)
line 
longitude/latitude of line as a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialLines* or SpatialPolygons* object 
length (in meters) for each line
For planar coordinates, see the terra or sf packages
line < rbind(c(180,20), c(150,10), c(140,55), c(10, 0), c(140,60)) d < lengthLine(line)
line < rbind(c(180,20), c(150,10), c(140,55), c(10, 0), c(140,60)) d < lengthLine(line)
Make a polygon or line by adding intermedate points (vertices) on the great circles inbetween the points supplied. This can be relevant when vertices are relatively far apart.
It can make the shape of the object to be accurate, when plotted on a plane. makePoly
will also close the polygon if needed.
makePoly(p, interval=10000, sp=FALSE, ...) makeLine(p, interval=10000, sp=FALSE, ...)
makePoly(p, interval=10000, sp=FALSE, ...) makeLine(p, interval=10000, sp=FALSE, ...)
p 
a 2column matrix (longitude/latitude) or a SpatialPolygons or SpatialLines object 
interval 
maximum interval of points, in units of r 
sp 
Logical. If 
... 
additional arguments passed to distGeo 
A matrix
Robert J. Hijmans
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) plot(pol) lines(pol, col='red', lwd=3) pol2 = makePoly(pol, interval=100000) lines(pol2, col='blue', lwd=2)
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) plot(pol) lines(pol, col='red', lwd=3) pol2 = makePoly(pol, interval=100000) lines(pol2, col='blue', lwd=2)
Transform longitude/latiude points to the Mercator projection. The main purpose of this function is to compute centroids, and to illustrate rhumb lines in the vignette.
mercator(p, inverse=FALSE, r=6378137)
mercator(p, inverse=FALSE, r=6378137)
p 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
inverse 
Logical. If 
r 
Numeric. Radius of the earth; default = 6378137 m 
matrix
Robert Hijmans
a = mercator(c(5,52)) a mercator(a, inverse=TRUE)
a = mercator(c(5,52)) a mercator(a, inverse=TRUE)
Find the point halfway between two points along an ellipsoid
midPoint(p1, p2, a=6378137, f = 1/298.257223563)
midPoint(p1, p2, a=6378137, f = 1/298.257223563)
p1 
longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
As above 
a 
major (equatorial) radius of the ellipsoid 
f 
ellipsoid flattening. The default value is for WGS84 
matrix with coordinate pairs
Elias Pipping and Robert Hijmans
midPoint(c(0,0),c(90,90)) midPoint(c(0,0),c(90,90), f=0)
midPoint(c(0,0),c(90,90)) midPoint(c(0,0),c(90,90), f=0)
Test if a point is on a great circle defined by two other points.
onGreatCircle(p1, p2, p3, tol=0.0001)
onGreatCircle(p1, p2, p3, tol=0.0001)
p1 
Longitude/latitude of the first point defining a great circle, in degrees; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object 
p2 
as above for the second point 
p3 
the point(s) to be tested if they are on the great circle or not 
tol 
numeric. maximum distance from the great circle (in degrees) that is tolerated to be considered on the circle 
logical
Robert Hijmans
onGreatCircle(c(0,0), c(30,30), rbind(c(10 11.33812), c(10,20)))
onGreatCircle(c(0,0), c(30,30), rbind(c(10 11.33812), c(10,20)))
Convert coordinates to the grid reference system used by the Ordnance Survey for Great Britain. Or do the inverse operation to get coordinates for a grid code.
OSGB(xy, precision, geo=FALSE, inverse=FALSE)
OSGB(xy, precision, geo=FALSE, inverse=FALSE)
xy 
x coordinate pairs (vector, matrix, data.frame 
; or grid codes if inverse=TRUE
.
precision 
character. One of "1m", "5m", "10m", "50m", "100m", "500m", "1km", "5km", "10km", "50km", "100km", "500km" 
geo 
If 
inverse 
If 
character
pnts < rbind(cbind(93555 , 256188), cbind(210637, 349798), cbind(696457, 481704)) g < OSGB(pnts, "1km", geo=FALSE) g OSGB(g, inverse=TRUE)
pnts < rbind(cbind(93555 , 256188), cbind(210637, 349798), cbind(696457, 481704)) g < OSGB(pnts, "1km", geo=FALSE) g OSGB(g, inverse=TRUE)
Compute the perimeter of a polygon (or the length of a line) with longitude/latitude coordinates, on an ellipsoid (WGS84 by default)
## S4 method for signature 'matrix' perimeter(x, a=6378137, f=1/298.257223563, ...) ## S4 method for signature 'SpatialPolygons' perimeter(x, a=6378137, f=1/298.257223563, ...) ## S4 method for signature 'SpatialLines' perimeter(x, a=6378137, f=1/298.257223563, ...)
## S4 method for signature 'matrix' perimeter(x, a=6378137, f=1/298.257223563, ...) ## S4 method for signature 'SpatialPolygons' perimeter(x, a=6378137, f=1/298.257223563, ...) ## S4 method for signature 'SpatialLines' perimeter(x, a=6378137, f=1/298.257223563, ...)
x 
Longitude/latitude of the points forming a polygon or line; Must be a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPolygons* or SpatialLines* object 
a 
major (equatorial) radius of the ellipsoid. The default value is for WGS84 
f 
ellipsoid flattening. The default value is for WGS84 
... 
Additional arguments. None implemented 
Numeric. The perimeter or length in m.
This function calls GeographicLib code by C.F.F. Karney
C.F.F. Karney, 2013. Algorithms for geodesics, J. Geodesy 87: 4355. doi:10.1007/s001900120578z. Addenda: https://geographiclib.sourceforge.io/geodaddenda.html. Also see https://geographiclib.sourceforge.io/
xy < rbind(c(180,20), c(140,55), c(10, 0), c(140,60), c(180,20)) perimeter(xy)
xy < rbind(c(180,20), c(140,55), c(10, 0), c(140,60), c(180,20)) perimeter(xy)
Plot polygons with arrow heads on each line segment, pointing towards the next vertex. This shows the direction of each line segment.
plotArrows(p, fraction=0.9, length=0.15, first='', add=FALSE, ...)
plotArrows(p, fraction=0.9, length=0.15, first='', add=FALSE, ...)
p 
Polygons (either a 2 column matrix or data.frame; or a SpatialPolygons* object 
fraction 
numeric between 0 and 1. When smaller then 1, interrupted lines are drawn 
length 
length of the edges of the arrow head (in inches) 
first 
Character to plot on first (and last) vertex 
add 
Logical. If 
... 
Additional arguments, see Details 
Based on an example in Software for Data Analysis by John Chambers (pp 250251) but adjusted such that the line segments follow great circles between vertices.
Robert J. Hijmans
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) plotArrows(pol)
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) plotArrows(pol)
randomCoordinates returns a 'uniform random sample' in the sense that the probability that a point is drawn from any region is equal to the area of that region divided by the area of the entire sphere. This would not happen if you took a random uniform sample of longitude and latitude, as the sample would be biased towards the poles.
regularCoordiaates returns a set of coordinates that are regularly distributed on the globe.
randomCoordinates(n) regularCoordinates(N)
randomCoordinates(n) regularCoordinates(N)
n 
Sample size (number of points (coordinate pairs)) 
N 
Number of 'parts' in which the earth is subdived ) 
Matrix of lon/lat coordiantes
Robert Hijmans, based on code by Nils Haeck (regularCoordinates) and on suggstions by Michael Orion (randomCoordinates)
randomCoordinates(3) regularCoordinates(1)
randomCoordinates(3) regularCoordinates(1)
This function returns a data.frame with parameters a
(semimajor axis) and 1/f
(inverse flattening) for a set of reference ellipsoids.
refEllipsoids()
refEllipsoids()
data.frame
To compute parameter b
you can do
Robert J. Hijmans
e < refEllipsoids() e[e$code=='WE', ] #to compute semiminor axis b: e$b < e$a  e$a / e$invf
e < refEllipsoids() e[e$code=='WE', ] #to compute semiminor axis b: e$b < e$a  e$a / e$invf
Compute the approximate surface span of polygons in longitude and latitude direction. Span is computed by rasterizing the polygons; and precision increases with the number of 'scan lines'. You can either use a fixed number of scan lines for each polygon, or a fixed bandwidth.
span(x, ...)
span(x, ...)
x 
a SpatialPolygons* object or a 2column matrix (longitude/latitude) 
... 
Additional arguments, see Details 
The following additional arguments can be passed, to replace default values for this function
nbands 
Character. Method to determine the number of bands to 'scan' the polygon. Either 'fixed' or 'variable'  
n 
Integer >= 1. If nbands='fixed' , how many bands should be used 

res 
Numeric. If nbands='variable' , what should the bandwidth be (in degrees)? 

fun 
Logical. A function such as mean or min. Mean computes the average span  
... 
further additional arguments passed to distGeo  
A list, or a matrix if a function fun
is specified. Values are in the units of r
(default is meter)
Robert J. Hijmans
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) plot(pol) lines(pol) # lon and lat span in m span(pol, fun=max) x < span(pol) max(x$latspan) mean(x$latspan) plot(x$longitude, x$lonspan)
pol < rbind(c(180,20), c(160,5), c(60, 0), c(160,60), c(180,20)) plot(pol) lines(pol) # lon and lat span in m span(pol, fun=max) x < span(pol) max(x$latspan) mean(x$latspan) plot(x$longitude, x$lonspan)
world coastline and country outlines in longitude/latitude (wrld) and in Mercator projection (merc).
data(wrld) data(merc)
data(wrld) data(merc)
Derived from the wrld_simpl data set in package maptools