# R/earth.poly.R In fossil: Palaeoecological and Palaeogeographical Analysis Tools

#### Defines functions `earth.poly`

````earth.poly` <-
function(lats) {
#tests if lats are a matrix or spatial points
if (inherits(lats, 'SpatialPoints')) lats<-coordinates(lats)
lats<-unique(lats)
rownames.lats<-rownames(lats)
nr <- nrow(lats)
#bearing matrix == bm == pairwise bearing matrix
bm <- matrix(, nr, nr)
for (i in 1:nr) {
for (j in 1:nr) {
if (i == j)
bm[i, j] <- -1
else bm[i, j] <- earth.bear(long1 = lats[i, 1], lat1 = lats[i, 2], long2 = lats[j, 1], lat2 = lats[j, 2])
}
}
#quadrant matrix == qm == which quadrants the bearings to other points fall in to
qm <- matrix(, nr, 4)
for (i in 1:nr) {
qm[i, 1] <- length(bm[1, ][bm[i, ] >= 0]) - length(bm[i,][bm[i, ] >= 90])
qm[i, 2] <- length(bm[1, ][bm[i, ] >= 90]) - length(bm[i,][bm[i, ] >= 180])
qm[i, 3] <- length(bm[1, ][bm[i, ] >= 180]) - length(bm[i,][bm[i, ] >= 270])
qm[i, 4] <- length(bm[1, ][bm[i, ] >= 270])
}
id <- numeric(nr)
#identifying which points have bearings in greater than 1 quadrant
#those that have bearings in less than 2 quadrants must be vertices (ie <180 degrees)
for (i in 1:nr) {
for (j in 1:4) {
if (length(qm[i, ][qm[i, ] > 0]) <= 2)
id[i] <- i
else {
if (length(qm[i, ][qm[i, ] > 0]) == 3) {
if (length(qm[i, 2:3][qm[i, 2:3] > 0]) == 2) {
if ((max(bm[i, ]) - min(bm[i, ][bm[i, ] >=
0])) < 180)
id[i] <- i
}
else {
if ((360 - min(bm[i, ][bm[i, ] > 180]) +
max(bm[i, ][bm[i, ] < 180])) < 180)
id[i] <- i
}
}
}
}
}
#mi == points on the perimeter
mi <- id[id > 0]
#the starting/primary point that will be used for the calculations
prpt <- mi[1]
#all points on the perimeter except the starting point
snpts <- mi[-1]
#the bearing for all the points from the primary point
pt <- bm[prpt, snpts]
#to get coordinates between -180 and 0 to between 180 and 360
if ((max(pt) - min(pt)) > 180) {
pt[pt < 180] <- pt[pt < 180] + 360
}
#ranked bearings, to give the points in the proper order so that the sequence of points is in order
rpt <- rank(pt)
ind <- numeric(length(pt))
for (i in 1:length(pt)) ind[rpt[i]] <- snpts[i]
ta <- numeric(1)
for (i in 1:(length(pt) - 1)) {
k <- i + 1
ta <- ta + earth.tri(long1 = lats[prpt, 1], lat1 = lats[prpt, 2], long2 = lats[ind[i], 1], lat2 = lats[ind[i], 2], long3 = lats[ind[k], 1], lat3 = lats[ind[k], 2])
}
verts <- c(prpt, ind)
return(list(area = ta, vertices = verts))
}
```

## Try the fossil package in your browser

Any scripts or data that you put into this service are public.

fossil documentation built on March 23, 2020, 5:06 p.m.