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


## ---- a1----------------------------------------------------------------------
if (!require("rspatial")) remotes::install_github('rspatial/rspatial')


## ---- a1b, error=TRUE---------------------------------------------------------
f <- system.file("WILDPOT.txt", package="rspatial")
basename(f)
d <- read.table(f, header=TRUE)


## ---- a2----------------------------------------------------------------------
# read all lines using UTF-8 encoding
d <- readLines(f, encoding='UTF-8')
# split each line into elements using the tabs
dd <- strsplit(d, '\t')
# show that the number of elements varies
table(sapply(dd, length))

# function to complete each line to 22 items
fun <- function(x) {
	r <- rep("", 22)
	r[1:length(x)] <- x
	r 
}

# apply function to each element of the list
ddd <- lapply(dd, fun)
# row bind all elements (into a matrix)
v <- do.call(rbind, ddd)
head(v)

#set the column names and remove them from the data
colnames(v) <- v[1,]
v <- v[-1,]

# coerce into a data.frame and change the type of some variables
# to numeric (instead of character)
v <- data.frame(v, stringsAsFactors=FALSE)


## ---- a3----------------------------------------------------------------------
# first coerce character values to numbers
for (i in c('LongD', 'LongM', 'LongS', 'LatD', 'LatM', 'LatS')) {
	v[, i] <- as.numeric(v[,i])
}
v$lon <- -1 * (v$LongD + v$LongM / 60 + v$LongS / 3600)
v$lat <- v$LatD + v$LatM / 60 + v$LatS / 3600

# Southern hemisphere gets a negative sign
v$lat[v$LatH == 'S'] <- -1 * v$lat[v$LatH == 'S']
head(v)


## ---- a4----------------------------------------------------------------------
library(raster)
library(rspatial)
cn <- sp_data('pt_countries')
cn


## ---- a5----------------------------------------------------------------------
plot(cn, xlim=c(-120, -40), ylim=c(-40,40), axes=TRUE)
points(v$lon, v$lat, cex=.5, col='red')


## ---- a6----------------------------------------------------------------------
sp <- v
coordinates(sp) <- ~lon + lat
proj4string(sp) <- CRS("+proj=longlat +datum=WGS84")


## ---- a7----------------------------------------------------------------------
sp <- SpatialPoints( v[, c('lon', 'lat')], 
             proj4string=CRS("+proj=longlat +datum=WGS84") )
sp <- SpatialPointsDataFrame(sp, v)


## ---- b1----------------------------------------------------------------------
table(v$COUNTRY)
# note Peru and PERU
v$COUNTRY <- toupper(v$COUNTRY)
table(v$COUNTRY)

# same fix for the SpatialPointsDataFrame
sp$COUNTRY <- toupper(sp$COUNTRY)


## ---- b2----------------------------------------------------------------------
ov <- over(sp, cn)
colnames(ov) <- 'name'
head(ov)
v <- cbind(v, ov)
table(v$COUNTRY)


## ---- b3----------------------------------------------------------------------
# some fixes first
# apparantly in the ocean (small island missing from polygon data)
v$name[is.na(v$name)] <- ''
# some spelling differenes
v$name[v$name=="UNITED STATES, THE"] <- "UNITED STATES"
v$name[v$name=="BRASIL"] <- "BRAZIL"

i <- which(toupper(v$name) != v$COUNTRY)
i
plot(cn, xlim=c(-120, -40), ylim=c(-40,40), axes=TRUE)
points(sp, cex=.25, pch='+', col='blue')
points(sp[i,], col='red', pch='x', cex=1.5)


## ---- b4----------------------------------------------------------------------
spc <- tapply(v$SPECIES, sp$COUNTRY, function(x)length(unique(x)) )
spc <- data.frame(COUNTRY=names(spc), nspp = spc)

# merge with country SpatialPolygonsDataFrame
cn <- merge(cn, spc, by='COUNTRY')
print(spplot(cn, 'nspp', col.regions=rev(terrain.colors(25))))


## ---- b5----------------------------------------------------------------------
tb <- table(v[ c('COUNTRY', 'SPECIES')])
# a big table
dim(tb)
# show two columns
tb[,2:3]


## ---- c1----------------------------------------------------------------------
library(rgdal)
# "proj.4" notation of CRS
projection(cn) <- "+proj=longlat +datum=WGS84"
# the CRS we want
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(cn, laea)
pts <- spTransform(sp, laea)
plot(clb, axes=TRUE)
points(pts, col='red', cex=.5)


## ---- d1----------------------------------------------------------------------
r <- raster(clb)
# 200 km = 200000 m
res(r) <- 200000


## ---- d2----------------------------------------------------------------------
rich <- rasterize(pts, r, 'SPECIES', function(x, ...) length(unique(na.omit(x))))
plot(rich)
plot(clb, add=TRUE)


## ---- d3----------------------------------------------------------------------
obs <- rasterize(pts, r, field='SPECIES', fun=function(x, ...)length((na.omit(x))) )
plot(obs)
plot(clb, add=TRUE)


## ---- d3b---------------------------------------------------------------------
plot(obs, rich, cex=1, xlab='Observations', ylab='Richness')


## ---- d4----------------------------------------------------------------------
d <- v[, c('lat', 'SPECIES')]
d$lat <- round(d$lat)
g <- tapply(d$SPECIES, d$lat, function(x) length(unique(na.omit(x))) )
plot(names(g), g)
# moving average
lines(names(g), movingFun(g, 3))


## ---- f1----------------------------------------------------------------------
# get the (Lambert AEA) coordinates from the SpatialPointsDataFrame
xy <- coordinates(pts)
# list of species
sp <- unique(pts$SPECIES)


## ---- f2----------------------------------------------------------------------
maxD <- vector(length=length(sp))
for (s in 1:length(sp)) {
	# get the coordinates for species 's'
	p <- xy[pts$SPECIES == sp[s], ]
	# distance matrix
	d <- as.matrix(dist(p))
	# ignore the distance of a point to itself
	diag(d) <- NA
	# get max value
	maxD[s] <- max(d, na.rm=TRUE)
}

# Note the typical J shape
plot(rev(sort(maxD))/1000, ylab='maxD (km)')


## ---- f3----------------------------------------------------------------------
library(dismo)
library(rgeos)
CA <- vector(length=length(sp))
for (s in 1:length(sp)) {
	p <- xy[pts$SPECIES == sp[s], ,drop=FALSE]
	# run "circles" model
	m <- try(circles(p, d=50000, lonlat=FALSE), silent=TRUE)
	if (!inherits(m, "try-error")) {
   	CA[s] <- area(polygons(m))
	}
}
# standardize to the size of one circle
CA <- CA / (pi * 50000^2)
plot(rev(sort(CA)), ylab='CA50')


## ---- f4----------------------------------------------------------------------
hull <- list()
for (s in 1:length(sp)) {
	p <- unique(xy[pts$SPECIES == sp[s], ,drop=FALSE])
	# need at least three points for hull
	if (nrow(p) > 3) {
		h <- convHull(p, lonlat=FALSE)
		pol <- polygons(h)
		hull[[s]] <- pol
	}
}


## ---- f4b---------------------------------------------------------------------
# which elements are NULL
i <- which(!sapply(hull, is.null))
h <- hull[i]
# combine them
hh <- do.call(bind, h)
plot(hh)


## ---- f4c---------------------------------------------------------------------
ahull <- sapply(hull, function(i) ifelse(is.null(i), 0, area(i)))
plot(rev(sort(ahull))/1000, ylab='Area of convex hull')


## ---- f5----------------------------------------------------------------------
  d <- cbind(maxD, CA, ahull)
  pairs(d)