## ----remotes------------------------------------------------------------------
if (!require("remotes")) install.packages("remotes")


## ----getDataOSU1, echo=TRUE---------------------------------------------------
if (!require("rspat")) remotes::install_github("rspatial/rspat")


## ----nzfrac1------------------------------------------------------------------
library(terra)
library(rspat)
coast <- spat_data("nz_coastline")
coast
plot(coast)


## ----crs----------------------------------------------------------------------
prj <- "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +datum=WGS84 +units=m"
mcoast <- project(coast, prj) 
mcoast


## ----yardfun------------------------------------------------------------------
stickpoints <- function(x, sticklength, lonlat) {
	nr <- nrow(x)
	pts <- 1
	pt <- 0
	sticklength <- sticklength * 1000
	while(TRUE) {
		pd <- distance(x[1,], x)
		# i is the first point further than the yardstick
		i <- which(pd > sticklength)[1]
		
		# if we cannot find a point within yardsitck distance we 
		# break out of the loop
		if (is.na(i)) break
		
		# remove the all points we have passed
		x <- x[(i+1):nrow(x), ]
		pt <- pt + i
		pts <- c(pts, pt)
	}
	pts
}


## ----computen-----------------------------------------------------------------
# get the x and y coordinates of the nodes
g <- as.points(mcoast)
# reverse the order (to start at the top rather than at the bottom)
g <- rev(g)

# three yardstick lengths
sticks <- c(50, 25, 10) # km
# create an empty list for the results
y <- list()
# loop over the yardstick lengths
for (i in 1:length(sticks)) {
	y[[i]] <- stickpoints(g, sticks[i], FALSE)
}
# These last four lines are equivalent to:
# y <- lapply(sticks, function(s) stickpoints(g, s, FALSE))


## ----fracplot1, fig.width = 10, fig.height = 6--------------------------------
n <- sapply(y, length)
par(mfrow=c(1,3))
for (i in 1:length(y)) {
	plot(mcoast)
	stops <- y[[i]]
	points(g[stops, ], col="red", pch=20, cex=2)
	lines(g[stops, ], col="red", lwd=3)
	text(1715000, 5860000, paste(n[i], "x", sticks[i], "=", n[i] * sticks[i], "km"), cex=1.5) 
}


## ----fracplot2, fig.width=6, fig.height=6-------------------------------------
plot(sticks, n, log="xy", cex=3, pch=20, col="red", 
	xlab="stick length", ylab="number of measures", las=1)
m <- lm(log(n) ~ log(sticks))
lines(sticks, exp(predict(m)), lwd=2, col="blue")
cf <- round(coefficients(m) , 3)
txt <- paste("log N =", cf[2], "log L +", cf[1])
text(6, 222, txt)


## -----------------------------------------------------------------------------
-cf[2]