## ---- echo=FALSE, include=FALSE-----------------------------------------------
library(rgdal)


## ---- frac1-------------------------------------------------------------------
library(raster)
uk <- raster::getData('GADM', country='GBR', level=0)
par(mai=c(0,0,0,0))
plot(uk)


## -----------------------------------------------------------------------------
data.frame(uk)


## -----------------------------------------------------------------------------
prj <- "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=WGS84 +units=m"


## -----------------------------------------------------------------------------
library(rgdal)
guk <- spTransform(uk, CRS(prj)) 


## -----------------------------------------------------------------------------
duk <- disaggregate(guk)
head(duk)


## -----------------------------------------------------------------------------
a <- area(duk)
i <- which.max(a)
a[i] / 1000000
b <- duk[i,]


## ---- frac5-------------------------------------------------------------------
par(mai=rep(0,4))
plot(b)


## -----------------------------------------------------------------------------
measure_with_ruler <- function(pols, length, lonlat=FALSE) {
	# some sanity checking
	stopifnot(inherits(pols, 'SpatialPolygons'))
	stopifnot(length(pols) == 1)
	
	# get the coordinates of the polygon
	g <- geom(pols)[, c('x', 'y')]
	nr <- nrow(g)
	
	# we start at the first point
	pts <- 1
	newpt <- 1
	while(TRUE) {
	    # start here
		p <- newpt
		
		# order the points 
		j <- p:(p+nr-1)
		j[j > nr] <- j[j > nr] - nr
		gg <- g[j,]
		
		# compute distances
		pd <- pointDistance(gg[1,], gg, lonlat)
		
		# get the first point that is past the end of the ruler
		# this is precise enough for our high resolution coastline
		i <- which(pd > length)[1]
		if (is.na(i)) {
			stop('Ruler is longer than the maximum distance found')
		}
		
		# get the record number for new point in the original order 
		newpt <- i + p
		
		# stop if past the last point
		if (newpt >= nr) break
		
		pts <- c(pts, newpt)
	}
	# add the last (incomplete) stick.
	pts <- c(pts, 1)
	# return the locations
	g[pts, ]
}


## -----------------------------------------------------------------------------
y <- list()
rulers <- c(25,50,100,150,200,250) # km
for (i in 1:length(rulers)) {
	y[[i]] <- measure_with_ruler(b, rulers[i]*1000)
}


## ---- frac15, fig.width = 8, fig.height = 10----------------------------------
par(mfrow=c(2,3), mai=rep(0,4))
for (i in 1:length(y)) {
	plot(b, col='lightgray', lwd=2)
	p <- y[[i]]
	lines(p, col='red', lwd=3)
	points(p, pch=20, col='blue', cex=2)
	
    bar <- rbind(cbind(525000, 900000), cbind(525000, 900000-rulers[i]*1000))
	lines(bar, lwd=2)
	points(bar, pch=20, cex=1.5)
	text(525000, mean(bar[,2]), paste(rulers[i], '  km'), cex=1.5)
	text(525000, bar[2,2]-50000, paste0('(', nrow(p), ')'), cex=1.25)
}


## ---- frac20, fig.width=6, fig.height=6---------------------------------------
# number of times a ruler was used
n <- sapply(y, nrow)

# set up empty plot
plot(log(rulers), log(n), type='n', xlim=c(2,6), ylim=c(2,6), axes=FALSE,
       xaxs="i",yaxs="i", xlab='Ruler length (km)', ylab='Number of segments')

# axes
tics <- c(1,10,25,50,100,200,400)
axis(1, at=log(tics), labels=tics)
axis(2, at=log(tics), labels=tics, las=2)

# linear regression line
m <- lm(log(n)~log(rulers))
abline(m, lwd=3, col='lightblue')

# add observations
points(log(rulers), log(n), pch=20, cex=2, col='red')


## ----fracplot, fig.width=6, fig.height=6--------------------------------------
small_rulers <- c(0.000001, 0.00001, 0.0001, 0.001, 0.01)  # km
nprd <- exp(predict(m, data.frame(rulers=small_rulers)))
coast <- nprd * small_rulers
plot(small_rulers, coast, xlab='Length of ruler', ylab='Length of coast', pch=20, cex=2, col='red')


## -----------------------------------------------------------------------------
m


## -----------------------------------------------------------------------------
-1 * m$coefficients[2]