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


## ----a1b, error=TRUE----------------------------------------------------------
f <- system.file("wildpot.csv", package="rspat")
basename(f)
v <- read.csv(f)


## ----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-----------------------------------------------------------------------
cn <- spat_data("pt_countries")
class(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 <- vect(v, crs="+proj=longlat +datum=WGS84")


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

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


## ----b2-----------------------------------------------------------------------
vv <- intersect(sp[, "COUNTRY"], cn)
names(vv)[1] <- "ptCountry"
head(vv)
table(vv$COUNTRY)


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

i <- which(toupper(vv$ptCountry) != vv$COUNTRY)
i
as.data.frame(vv[i,])
plot(cn, xlim=c(-120, -40), ylim=c(-40,40), axes=TRUE)
points(sp, cex=.25, pch='+', col='blue')

points(vv[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 SpatVector --- fix the names in the polygons this time
cn$COUNTRY[cn$COUNTRY=="UNITED STATES, THE"] <- "UNITED STATES"
cn$COUNTRY[cn$COUNTRY=="BRASIL"] <- "BRAZIL"

cns <- merge(cn, spc, by="COUNTRY", all.x=TRUE)
plot(cns, "nspp", col=rev(terrain.colors(25)), breaks=c(1,5,10,20,30,40,90))


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


## ----c1-----------------------------------------------------------------------
# the CRS we want
laea <-"+proj=laea  +lat_0=0 +lon_0=-80"
clb <- project(cn, laea)
pts <- project(sp, laea)
plot(clb)
points(pts, col='red', cex=.5)


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


## ----d2-----------------------------------------------------------------------
rich <- rasterize(pts, r, "SPECIES", function(x, ...) length(unique(na.omit(x))))
plot(rich)
lines(clb)


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


## ----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), raster::movingFun(g, 3))


## ----f2-----------------------------------------------------------------------
spp <- unique(pts$SPECIES)
maxD <- rep(NA, length(spp))
for (s in 1:length(spp)) {
	# get the coordinates for species 's'
	p <- pts[pts$SPECIES == spp[s], ]
	if (nrow(p) < 2) next
	# distance matrix
	d <- as.matrix(distance(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-----------------------------------------------------------------------
CA <- rep(NA, length(spp))
for (s in 1:length(spp)) {
	p <- pts[pts$SPECIES == spp[s], ]
	# run "circles" model
	m <- aggregate(buffer(p, 50000))
	CA[s] <- expanse(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(spp)) {
	p <- unique(pts[pts$SPECIES == spp[s], ])
	# need at least three (unique) points for hull
	if (nrow(p) > 3) {
		h <- convHull(p)
		if (geomtype(h) == "polygons") {
  		hull[[s]] <- h
		}
	}
}


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


## ----f4c----------------------------------------------------------------------
ahull <- expanse(hh)
plot(rev(sort(ahull))/1000, ylab="Area of convex hull")


## -----------------------------------------------------------------------------
cHull <- rep(NA, length(spp))
cHull[i] <- ahull


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