## ---- echo=FALSE, include=FALSE-----------------------------------------------
library(knitr)
opts_chunk$set(fig.width = 4, fig.height = 5, fig.cap = '', collapse = TRUE)
library(raster)
library(rgdal)


## ----getDataOSU1, echo=TRUE---------------------------------------------------
if (!require("devtools")) install.packages('devtools')
if (!require("rspatial")) remotes::install_github('rspatial/rspatial')


## ---- nzfrac1-----------------------------------------------------------------
library(rspatial)
coast <- sp_data("nz_coastline")
coast
par(mai=c(0,0,0,0))
plot(coast)


## ----crs----------------------------------------------------------------------
prj <- "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m"
library(rgdal)
mcoast <- spTransform(coast, CRS(prj)) 
mcoast


## ----yardfun------------------------------------------------------------------
stickpoints <- function(x, sticklength, lonlat) {
	# x is a matrix with two columns (x and y coordinates)
	nr <- nrow(x)
	pts <- 1
	pt <- 0
	while(TRUE) {
		pd <- pointDistance(x[1,], x, lonlat)
		# 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
g <- geom(mcoast)[, c('x', 'y')]
# reverse the order (to start at the top rather than at the bottom)
g <- g[nrow(g):1, ]

# three yardstick lengths
sticks <- c(10, 5, 2.5) # 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]*1000, FALSE)
}
# These last four lines are equivalent to:
# y <- lapply(sticks*1000, 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]]
	lines(g[stops, ], col="red", lwd=2)
	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]