The method, which may be useful in determining coordinate systems for a mooring array or a ship transects, calculates (x,y) from distance calculations along geodesic curves. See “Caution”.
vector of longitude and latitude
numeric reference location. Poor results
will be returned if these values are not close to the locations described
an integer specifying whether debugging information is
to be printed during the processing. This is a general parameter that
is used by many
The calculation is as follows.
i-th point in the
vectors. To calculate
used is to find the distance along a
geodesic curve connecting (
latitude[i]). The resultant distance
is multiplied by -1 if
longitude[i]-longitudeRef is negative,
and the result is assigned to
A similar procedure is used for
geodXy returns a data frame of
geodesic distance components, measured in metres.
This scheme is without known precedent in the literature, and users should read the documentation carefully before deciding to use it.
On 2015-11-02, the names of the arguments were changed from
lon, etc., to
longitude, etc., to be in keeping with other oce functions.
On 2017-04-05, four changes were made.
1. Default values of
latitudeRef were removed,
since the old defaults were inappropriate to most work.
2. The argument called
rotate was eliminated, because it only made
sense if the mean resultant x and y were zero.
3. The example was made more useful.
4. Pointers were made to
lonlat2utm, which may be more useful.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## Not run: # Develop a transect-based axis system for final data(section) stations library(oce) data(section) lon <- tail(section[["longitude", "byStation"]], 26) lat <- tail(section[["latitude", "byStation"]], 26) lonR <- tail(lon, 1) latR <- tail(lat, 1) data(coastlineWorld) mapPlot(coastlineWorld, proj="+proj=merc", longitudelim=c(-75,-65), latitudelim=c(35,43), col="gray") mapPoints(lon, lat) XY <- geodXy(lon,lat,mean(lon), mean(lat)) angle <- 180/pi*atan(coef(lm(y~x, data=XY))) mapCoordinateSystem(lonR, latR, 500, angle, col=2) # Compare UTM calculation UTM <- lonlat2utm(lon, lat, zone=18) # we need to set the zone for this task! angleUTM <- 180/pi*atan(coef(lm(northing~easting, data=UTM))) mapCoordinateSystem(lonR, latR, 500, angleUTM, col=3) legend("topright", lwd=1, col=2:3, bg="white", title="Axis Rotation Angle", legend=c(sprintf("geod: %.1f deg", angle), sprintf("utm: %.1f deg",angleUTM))) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.