.adjFormat <- function(x) {
x <- x[, 1:3]
names(x) <- c("species", "decimallongitude", "decimallatitude")
return(x)
}
.BarChartPoly <- function(x) {
liste <- names(x$spec_table[-1])
leng <- length(liste)
if (leng == 0) {
cat("No point fell in any polygon")
} else {
for (i in 1:leng) {
subs <- subset(x$spec_table, x$spec_table[, i] > 0)
datsubs <- subs[order(subs[, i]), ]
if (nrow(subs) == 0) {
plot(1:10, 1:10, type = "n", xlab = "", ylab = "Number of occurences")
text(3, 6, labels = "No species occurred in this polygon.",
adj = 0)
title(liste[i])
} else {
out <- ggplot() + ggplot2::geom_bar(data = datsubs, aes_string(x = "rownames(datsubs)",
y = "datsubs[,i]"), stat = "identity") + ggplot2::theme_bw() +
ggplot2::xlab("Species") + ggplot2::ylab("Number of occurrences") +
ggplot2::ggtitle(liste[i])
print(out)
}
}
}
}
.BarChartSpec <- function(x, mode = c("percent", "total")) {
match.arg(mode)
switch(mode, percent = {
dat.plo <- x$spec_table/rowSums(x$spec_table) * 100
leng <- length(rownames(dat.plo))
for (i in 1:leng) {
dat.sub <- as.numeric(as.vector(dat.plo[i, ]))
out <- ggplot2::ggplot() + ggplot2::geom_bar(data = data.frame(dat.sub),
aes_string(x = "names(x$spec_table)", y = "dat.sub"), stat = "identity") +
ggplot2::theme_bw() + ggplot2::ylab("Percent of occurrences") +
ggtitle(rownames(dat.plo[i, ])) + theme(axis.title.x = element_blank())
print(out)
}
}, total = {
dat.plo <- x$spec_table
leng <- length(rownames(dat.plo))
for (i in 1:leng) {
dat.sub <- as.numeric(as.vector(dat.plo[i, ]))
out <- ggplot2::ggplot() + ggplot2::geom_bar(data = data.frame(dat.sub),
aes_string(x = "names(x$spec_table)", y = "dat.sub"), stat = "identity") +
ggplot2::theme_bw() + ggplot2::ylab("Number of occurrences") +
ggtitle(rownames(dat.plo[i, ])) + theme(axis.title.x = element_blank())
print(out)
}
})
}
.ConvHull <- function(x, type) {
if (type == "euclidean") {
conv.hull <- grDevices::chull(x$decimallongitude, x$decimallatitude)
dat2 <- x[conv.hull, ]
out <- rbind(dat2[, c(2, 3)], dat2[1, c(2, 3)])
}
if (type == "pseudospherical") {
cv <- grDevices::chull(x$decimallongitude, x$decimallatitude)
dat2 <- x[cv, ]
dat2 <- rbind(dat2[, c(2, 3)], dat2[1, c(2, 3)])
out <- geosphere::makePoly(dat2)
}
poly <- sp::Polygons(list(Polygon(out)), ID = paste(x[1, 1], "_convhull",
sep = ""))
poly <- sp::SpatialPolygons(list(poly), proj4string = CRS("+proj=longlat +datum=WGS84"))
return(poly)
}
.ConvArea <- function(x, reps, repfrac, repsize, terrestrial, biome, type, cropper) {
if (!is.null(repfrac)) {
if (repfrac > 1) {
stop("'repfrac' must be between 0 and 1")
}
repsize <- round(nrow(x[!duplicated(x), ]) * repfrac)
if (repsize < 5 & reps > 1) {
repsize <- 5
warning("'repfrac' to small, repsize set to 5")
}
}
if (nrow(x[!duplicated(x), ]) > repsize) {
samp <- as.list(data.frame((replicate(reps, sample(1:nrow(x), size = repsize,
replace = FALSE)))))
samp <- lapply(samp, function(k) {
x[k, ]
})
pols <- lapply(samp, ".ConvHull", type = type)
nam <- names(pols)
names(pols) <- NULL
if (length(pols) > 1) {
pols <- do.call(raster::bind, pols)
} else {
slot(slot(pols[[1]], "polygons")[[1]], "ID") <- "1"
pols <- pols[[1]]
}
pols <- sp::SpatialPolygonsDataFrame(pols, data = data.frame(species = nam))
if (terrestrial) {
pols <- rgeos::gIntersection(pols, cropper, byid = TRUE)
}
if (!is.null(biome)) {
if (!"BIOME" %in% names(biome)) {
stop("'BIOME' not found in 'biome'")
}
pts <- sp::SpatialPoints(x[, c("decimallongitude", "decimallatitude")])
biome <- raster::crop(biome, extent(pts))
test <- sp::over(pts, biome)
overl <- biome[biome$BIOME %in% test$BIOME, ]
overl <- aggregate(overl)
pols <- rgeos::gIntersection(pols, overl, byid = TRUE)
pols <- gBuffer(pols, byid = TRUE, width = 0)
}
pol.area <- geosphere::areaPolygon(pols)
pol.area <- median(pol.area)
} else {
pols <- .ConvHull(x, type = type)
if (terrestrial) {
pols <- rgeos::gIntersection(pols, cropper, byid = TRUE)
pols <- rgeos::gBuffer(pols, byid = TRUE, width = 0)
}
if (!is.null(biome)) {
if (!"BIOME" %in% names(biome)) {
stop("'BIOME' not found in 'biome'")
}
pts <- sp::SpatialPoints(x[, c("decimallongitude", "decimallatitude")])
biome <- raster::crop(biome, extent(pts))
test <- sp::over(pts, biome)
overl <- biome[biome$BIOME %in% test$BIOME, ]
overl <- aggregate(overl)
pols <- rgeos::gIntersection(pols, overl, byid = TRUE)
pols <- rgeos::gBuffer(pols, byid = TRUE, width = 0)
}
pol.area <- geosphere::areaPolygon(pols)
}
pol.area <- round(pol.area/(1000 * 1000), 0)
return(pol.area)
}
.ConvertPoly <- function(x) {
x <- read.table(x, sep = "\t")
out2 <- vector()
for (j in 1:dim(x)[1]) {
aa <- as.character(x[j, ])
ff <- t(aa)
bb <- unlist(strsplit(ff[1], split = ":"))
bb <- c(bb[1], unlist(strsplit(bb[2], split = " ")))
out <- c(1, 1, 1)
for (i in 2:length(bb)) {
dd <- c(bb[1], unlist(strsplit(as.character(bb[i]), split = ",")))
out <- rbind(out, dd)
}
out2 <- rbind(out2, out[-c(1, 2), ])
}
colnames(out2) <- c("species", "decimallongitude", "decimallatitude")
rownames(out2) <- NULL
out2 <- as.data.frame(out2)
return(out2)
}
.Cord2Polygon <- function(x) {
if (is.character(x)) {
tt <- read.table(x, sep = "\t")
} else {
tt <- x
}
dat <- split(tt, f = tt[, 1])
col <- lapply(dat, function(k) {
pp <- Polygon(k[, 2:3])
po <- Polygons(list(pp), ID = unique(k[, 1]))
})
polys <- SpatialPolygons(col)
pol.data <- data.frame(unique(tt[, 1]), row.names = unique(tt[, 1]))
names(pol.data) <- names(x)[1]
polys <- SpatialPolygonsDataFrame(polys, data = pol.data)
return(polys)
}
.getEle <- function(x) {
ele <- try(getData("SRTM", lon = round(as.numeric(x[2]), 2), lat = round(as.numeric(x[3]),
2)))
if (class(ele) == "try-error") {
elevation <- "NA"
} else {
if (!is.na(extract(ele, SpatialPoints(data.frame(round(as.numeric(x[2]),
2), round(as.numeric(x[3]), 2)))))) {
elevation <- extract(ele, data.frame(round(as.numeric(x[2]), 2),
round(as.numeric(x[3]), 2)))
} else {
elevation <- "NA"
}
}
return(elevation)
}
.GetElevation <- function(x) {
if (is.data.frame(x)) {
inp <- x
}
if (class(x) == "spgeoOUT") {
inp <- data.frame(species = x$species, decimallongitude = x$species_coordinates[,
1], decimallatitude = x$species_coordinates[, 2])
}
if (is.character(x) & length(grep(".txt", x)) == 0) {
if (!requireNamespace("rgbif", quietly = TRUE)) {
stop("rgbif needed for species name option. Please install it.",
call. = FALSE)
}
coords <- rgbif::occ_search(scientificName = x, return = "data", limit = 2e+05,
hasCoordinate = TRUE, spatialIssues = FALSE, fields = c("species", "decimalLongitude",
"decimalLatitude"))
coords <- do.call("rbind", coords)
names(coords) <- c("species", "decimallongitude", "decimallatitude")
coords <- data.frame(coords[complete.cases(coords), ])
warning(paste(dim(inp)[1], "geo-referenced records found in GBIF; no data cleaning was performed",
sep = " "))
}
if (is.character(x) & length(grep(".txt", x)) > 0) {
inp <- read.table(x, sep = "\t", header = TRUE)
names(inp) <- c("species", "decimallongitude", "decimallatitude")
}
tt <- list()
for (i in 1:dim(inp)[1]) {
tt[[i]] <- .getEle(inp[i, ])
}
ele.vector <- suppressWarnings(as.numeric(unlist(tt)))
if (is.character(x) & length(grep(".txt", x)) == 0) {
ele.vector <- cbind(inp, ele.vector)
return(ele.vector)
} else {
return(ele.vector)
}
}
.MapAll <- function(x, buffer = 1) {
dat <- x$samples
dat$homepolygon <- as.character(dat$homepolygon)
dat$homepolygon[dat$homepolygon != "not_classified"] <- "classified"
# prepare background
e <- raster::extent(SpatialPoints(x$samples[, 2:3])) + buffer
bgmap <- speciesgeocodeR::landmass
bgmap <- raster::crop(bgmap, e)
bgmap <- ggplot2::fortify(bgmap)
pols <- ggplot2::fortify(x$polygons)
# plot results
plo <- ggplot2::ggplot() + ggplot2::geom_polygon(data = bgmap, aes_string(x = "long",
y = "lat", group = "group"), fill = "grey60") + ggplot2::geom_polygon(data = pols,
aes_string(x = "long", y = "lat", group = "group"), fill = rgb(0, 100,
0, 100, maxColorValue = 255)) + ggplot2::geom_point(data = dat,
aes_string(x = "decimallongitude", y = "decimallatitude", color = "homepolygon")) +
ggplot2::scale_colour_manual(values = c("blue", "red")) + ggplot2::coord_fixed() +
ggplot2::theme_bw() + theme(legend.title = element_blank())
print(plo)
}
.MapPerPoly <- function(x, buffer = 1) {
# background plot data
e <- raster::extent(SpatialPoints(x$samples[, 2:3])) + buffer
areanames <- x$areanam
bgmap <- speciesgeocodeR::landmass
bgmap <- raster::crop(bgmap, e)
bgmap <- ggplot2::fortify(bgmap)
# backgroundplot
plo <- ggplot2::ggplot() + ggplot2::geom_polygon(data = bgmap, aes_string(x = "long",
y = "lat", group = "group"), fill = "grey60") + ggplot2::coord_fixed() +
ggplot2::theme_bw()
# per polygon plots
liste <- unique(as.character(x$polygons@data[, areanames]))
outp.li <- list()
for (i in liste) {
# select polygon
pols <- x$polygons[x$polygons[[areanames]] == i, ]
pols <- ggplot2::fortify(pols)
# select points
pts <- x$samples[x$samples$homepolygon == i, ]
# plot
plo2 <- plo + ggplot2::geom_polygon(data = pols, aes_string(x = "long",
y = "lat", group = "group"), fill = "grey90") + ggplot2::geom_point(data = pts,
aes_string(x = "decimallongitude", y = "decimallatitude", colour = "species")) +
ggplot2::ggtitle(i) + theme(legend.position = "right")
outp.li[[i]] <- plo2
}
lapply(outp.li, "print")
}
.MapPerSpecies <- function(x, buffer = 1) {
# create background plot
e <- raster::extent(SpatialPoints(x$samples[, 2:3])) + buffer
bgmap <- speciesgeocodeR::landmass
bgmap <- raster::crop(bgmap, e)
bgmap <- ggplot2::fortify(bgmap)
pols <- ggplot2::fortify(x$polygons)
# plot results
plo <- ggplot2::ggplot() + ggplot2::geom_polygon(data = bgmap, aes_string(x = "long",
y = "lat", group = "group"), fill = "grey60") + ggplot2::geom_polygon(data = pols,
aes_string(x = "long", y = "lat", group = "group"), fill = rgb(0, 100,
0, 100, maxColorValue = 255)) + ggplot2::coord_fixed() + ggplot2::theme_bw() +
theme(legend.title = element_blank())
# create per species plots
inp <- x$samples
inp$homepolygon <- as.character(inp$homepolygon)
inp$homepolygon[inp$homepolygon != "not_classified"] <- "classified"
liste <- sort(unique(inp$species))
for (i in liste) {
dat <- inp[inp$species == i, ]
plo2 <- plo + ggplot2::geom_point(data = dat, aes_string(x = "decimallongitude",
y = "decimallatitude", color = "homepolygon")) + ggplot2::scale_colour_manual(values = c("blue",
"red")) + ggplot2::ggtitle(i)
print(plo2)
}
}
.MapUnclassified <- function(x, buffer = 1) {
# pick unclassified species
dat <- x$samples[x$samples$homepolygon == "not_classified", ]
if (nrow(dat) == 0) {
plot(c(1:20), c(1:20), type = "n", axes = FALSE, xlab = "", ylab = "")
text(10, 10, labels = paste("All points fell into the polygons and were classified.\n",
"No unclassified points", sep = ""))
} else {
# prepare background
e <- raster::extent(SpatialPoints(x$samples[, 2:3])) + buffer
bgmap <- speciesgeocodeR::landmass
bgmap <- raster::crop(bgmap, e)
bgmap <- ggplot2::fortify(bgmap)
pols <- ggplot2::fortify(x$polygons)
# plot results
plo <- ggplot2::ggplot() + ggplot2::geom_polygon(data = bgmap, aes_string(x = "long",
y = "lat", group = "group"), fill = "grey60") + ggplot2::geom_polygon(data = pols,
aes_string(x = "long", y = "lat", group = "group"), fill = rgb(0,
100, 0, 100, maxColorValue = 255)) + ggplot2::geom_point(data = dat,
aes_string(x = "decimallongitude", y = "decimallatitude", color = "species")) +
ggplot2::coord_fixed() + ggplot2::theme_bw()
print(plo)
}
}
.NexusOut <- function(dat, verbose = FALSE) {
if (!is.spgeoOUT(dat)) {
tablist <- lapply(dat, function(x) x$spec_table)
for (i in 1:length(tablist)) {
names(tablist[[i]])[-1] <- paste(names(tablist)[i], names(tablist[[i]][-1]),
sep = "_")
}
speciestab <- Reduce(function(x, y) merge(x, y, all = TRUE), tablist)
} else {
speciestab <- dat$spec_table
}
if (!verbose) {
sink("species_classification.nex")
}
if (verbose) {
sink("species_classification_verbose.nex")
}
cat("#NEXUS \n")
cat("\n")
cat("begin data; \n")
cat(paste("\tdimensions ntax=", dim(speciestab)[1], " nchar=", dim(speciestab)[2] -
1, ";", sep = ""))
cat("\n")
cat("\tformat datatype=standard symbols=\"01\" gap=-;")
cat("\n")
cat("\tCHARSTATELABELS")
cat("\n")
if (length(speciestab) == 0) {
cat("No point fell in any of the polygons specified")
sink(NULL)
} else {
aa <- gsub(" ", "_", names(speciestab))
aa <- gsub("&", "_", aa)
aa <- gsub("__", "_", aa)
aa <- gsub("__", "_", aa)
bb <- seq(1, length(aa))
cat(paste("\t", bb[-length(bb)], " ", aa[-length(aa)], ",\n", sep = ""))
cat("\t", paste(bb[length(bb)], " ", aa[length(aa)], ";\n", sep = ""))
cat("\n")
cat("\tmatrix\n")
dd <- as.matrix(speciestab)
dd[dd > 0] <- 1
if (ncol(dd) > 1) {
dd <- data.frame(dd)
dd$x <- apply(dd[, names(dd)], 1, paste, collapse = "")
} else {
dd <- data.frame(dd, x = dd)
}
ff <- gsub(" ", "_", rownames(speciestab))
if (!verbose) {
ee <- paste("\t\t", ff, "\t", dd$x, "\n", sep = "")
cat(ee)
}
if (verbose) {
gg <- vector()
jj <- speciestab
for (i in 1:ncol(jj)) {
hh <- paste(dd[, i], "[", jj[, i], "]", sep = "")
gg <- data.frame(cbind(gg, hh))
}
gg$x <- apply(gg[, names(gg)], 1, paste, collapse = "")
ee <- paste("\t\t", ff, "\t", gg$x, "\n", sep = "")
cat(ee)
}
cat("\t;\n")
cat("end;")
sink(NULL)
}
}
.OutBarChartPoly <- function(x, path, prefix, verbose = FALSE) {
if (verbose) {
cat("Creating barchart per polygon: barchart_per_polygon.pdf. \n")
}
if (missing(path)) {
path <- getwd()
}
pdf(file = file.path(path, paste(prefix, "barchart_per_polygon.pdf", sep = "")),
paper = "special", width = 10.7, height = 7.2, onefile = TRUE)
suppressMessages(.BarChartPoly(x))
dev.off()
}
.OutBarChartSpec <- function(x, path, prefix, verbose = FALSE, ...) {
if (verbose) {
cat("Creating barchart per species: barchart_per_species.pdf. \n")
}
if (missing(path)) {
path <- getwd()
}
pdf(file = file.path(path, paste(prefix, "barchart_per_species.pdf", sep = "")),
paper = "special", width = 10.7, height = 7.2, onefile = TRUE)
suppressMessages(.BarChartSpec(x, mode = "percent"))
dev.off()
}
.OutBayArea <- function(x, prefix) {
# taxon to area classification
dat <- x[["spec_table"]]
# exclude not classified column
dat <- dat[!names(dat) == "not_classified"]
# set from occurrence number to present/absent
dat[dat > 0] <- 1
# dat[, 1] <- rownames(dat)
dat <- cbind(sp = rownames(dat), dat)
dat$sp <- as.character(dat$sp)
out <- rbind(rep("", ncol(dat)), dat)
out[, 1] <- paste(out[, 1], "\t", sep = "")
out[1, 1] <- paste(nrow(dat), sum(names(out) != "sp"))
write.table(out, paste(prefix, "bayarea_data.txt", sep = "_"), sep = "",
na = "", col.names = FALSE, row.names = FALSE, quote = FALSE)
# area coordinates
coor <- x[["polygons"]]
coor <- data.frame(sp::coordinates(coor))
coor <- coor[, c(2, 1)]
coor <- rbind(c("# 0.0", ""), coor)
write.table(coor, paste(prefix, "bayarea_areas.txt", sep = "_"), sep = "\t",
na = "", col.names = FALSE, row.names = FALSE, quote = FALSE)
}
.OutMapAll <- function(x, path, prefix, areanames = "", verbose = FALSE) {
if (verbose) {
cat("Creating overview map: map_samples_overview.pdf. \n")
}
if (missing(path)) {
path <- getwd()
}
pdf(file = file.path(path, paste(prefix, "map_samples_overview.pdf", sep = "")),
paper = "special", width = 10.7, height = 7.2, onefile = TRUE)
suppressMessages(.MapAll(x))
suppressMessages(.MapUnclassified(x))
dev.off()
}
.OutMapPerPoly <- function(x, path, prefix, verbose = FALSE) {
if (verbose) {
cat("Creating map per polygon: map_samples_per_polygon.pdf. \n")
}
if (missing(path)) {
path <- getwd()
}
pdf(file.path(path, file = paste(prefix, "map_samples_per_polygon.pdf",
sep = "")), paper = "special", width = 10.7, height = 7.2, onefile = TRUE)
suppressMessages(.MapPerPoly(x))
dev.off()
}
.OutMapPerSpecies <- function(x, path, prefix, verbose = FALSE) {
if (verbose) {
cat("Creating map per species: map_samples_per_species.pdf. \n")
}
if (missing(path)) {
path <- getwd()
}
pdf(file = file.path(path, paste(prefix, "map_samples_per_species.pdf",
sep = "")), paper = "special", width = 10.7, height = 7.2, onefile = TRUE)
.MapPerSpecies(x)
dev.off()
}
.OutPlotSpPoly <- function(x, path, prefix, verbose = FALSE) {
if (verbose) {
cat("Creating species per polygon barchart: number_of_species_per_polygon.pdf. \n")
}
if (missing(path)) {
path <- getwd()
}
pdf(file = file.path(path, paste(prefix, "number_of_species_per_polygon.pdf",
sep = "")), paper = "special", width = 10.7, height = 7.2, onefile = TRUE)
suppressMessages(.PlotSpPoly(x))
dev.off()
}
.PipSamp <- function(x, columnname, verbose = FALSE) {
if (class(x$polygons) == "SpatialPolygonsDataFrame") {
if (any(is.na(x$polygons@data[, columnname]))) {
stop("area names contain missing data (#N/A). Renamed to unnamed. This migh cause problems")
}
occ <- SpatialPoints(x$species_coordinates[, c(1, 2)], proj4string = sp::CRS(sp::proj4string(x$polygons)))
outp <- as.character(over(occ, x$polygon)[, columnname])
outp[is.na(outp)] <- "not_classified"
return(outp)
} else {
if (verbose) {
cat("Performing point in polygon test \n")
}
pip <- over(occ, x$polygons)
if (verbose) {
cat("Done \n")
}
for (i in 1:length(names(x$polygons))) {
pip$homepolygon[pip$homepolygon == i] <- names(x$polygons)[i]
}
return(pip)
}
}
.PlotSpPoly <- function(x) {
dat.plo <- data.frame(x$polygon_table)
out <- ggplot() + ggplot2::geom_bar(data = dat.plo, aes_string(x = "rownames(dat.plo)",
y = "dat.plo[,1]"), stat = "identity") + ggplot2::theme_bw() + ggplot2::ylab("Species") +
ggplot2::theme(axis.title.x = element_blank(), axis.text.x = element_text(angle = 90,
hjust = 1))
print(out)
}
.rasterSum <- function(x, ras, type) {
po <- sp::SpatialPoints(x[, 2:3])
ras_sub <- raster::rasterize(po, ras, fun = "count")
if (type == "div") {
ras_sub[ras_sub >= 1] <- 1
}
ras_sub[is.na(ras_sub)] <- 0
return(ras_sub)
}
.ReadPoints <- function(x, y, areanames = NA, verbose = FALSE) {
res <- list()
if (all(!is.character(x), !is.data.frame(x))) {
stop(sprintf("function not defined for class %s", dQuote(class(x))))
}
if (is.character(x) & length(grep(".txt", x)) == 0) {
if (!requireNamespace("rgbif", quietly = TRUE)) {
stop("rgbif needed for species name option. Please install it.",
call. = FALSE)
}
coords <- rgbif::occ_search(scientificName = x, return = "data", limit = 2e+05,
hasCoordinate = TRUE, hasGeospatialIssue = FALSE, fields = c("species",
"decimalLongitude", "decimalLatitude"))
if (length(x) == 1) {
coords <- do.call("cbind.data.frame", coords)
} else {
coords <- do.call("rbind.data.frame", coords)
}
coords <- data.frame(coords[complete.cases(coords), ])
}
if (is.character(x) & length(grep(".txt", x)) > 0) {
coords <- read.table(x, sep = "\t", header = TRUE, row.names = NULL)
}
if (is.data.frame(x)) {
names(x) <- tolower(names(x))
if (ncol(x) > 3 & all(c("species", "decimallatitude", "decimallongitude") %in%
names(x))) {
coords <- x[, c("species", "decimallongitude", "decimallatitude")]
rownames(coords) <- 1:nrow(coords)
} else {
coords <- x
rownames(coords) <- 1:nrow(coords)
}
if (ncol(coords) < 3) {
stop(paste("wrong input format: \n", "Inputfile for coordinates must have at least three columns",
sep = ""))
}
}
if (is.character(y) | is.data.frame(y)) {
if (is.character(y) & length(grep(".shp", y)) > 0) {
poly <- maptools::readShapeSpatial(y)
} else {
if (is.character(y)) {
polycord <- read.table(y, sep = "\t", header = TRUE)
}
if (is.data.frame(y)) {
polycord <- y
}
if (ncol(polycord) != 3) {
stop("Wrong input format;\ninputfile for polygons must be a tab-delimited text file with three columns")
}
if (!is.numeric(polycord[, "decimallongitude"]) || !is.numeric(polycord[,
"decimallatitude"])) {
stop("wrong input format:\nInput polygon coordinates (columns 2 and 3) must be numeric.")
}
if (!is.character(polycord[, "species"]) && !is.factor(polycord[,
"species"])) {
warning("polygon identifier (column 1) should be a string or a factor")
}
if (max(polycord[, "decimallongitude"]) > 180) {
warning(sprintf("check polygon input coordinates; file contains longitude values outside possible range in row: \n %s\n Coordinates set to maximum: 180.\n",
rownames(polycord[polycord[, 2] > 180, ])))
polycord[polycord[, "decimallongitude"] > 180, "decimallongitude"] <- 180
}
if (min(polycord[, "decimallongitude"]) < -180) {
warning(paste("check polygon input coordinates. File contains longitude values outside possible range in row: ",
rownames(polycord[polycord[, "decimallongitude"] < -180, ]),
"\n", "Coordinates set to minimum: -180", sep = ""))
polycord[polycord[, "decimallongitude"] < -180, ] <- -180
}
if (max(polycord[, "decimallatitude"]) > 90) {
warning(paste("check polygon input coordinates. File contains latitude values outside possible range in row:",
rownames(polycord[polycord[, "decimallatitude"] > 90, ]),
"\n", "Coordinates set to maximum: 90", sep = ""))
polycord[polycord[, "decimallatitude"] > 90, "decimallatitude"] <- 90
}
if (min(polycord[, "decimallatitude"]) < -90) {
warning(paste("check polygon input coordinates. File contains latitude values outside possible range in row:",
rownames(polycord[polycord[, "decimallatitude"] < -90, ]),
"\n", "Coordinates set to minimum: -90", sep = ""))
polycord[polycord[, "decimallatitude"] < -90, "decimallatitude"] <- -90
}
poly <- .Cord2Polygon(polycord)
}
}
if (class(y) == "SpatialPolygonsDataFrame") {
poly <- y
}
if (!is.numeric(coords[, "decimallongitude"]) || !is.numeric(coords[, "decimallatitude"])) {
stop(paste("wrong input format: \n", "Input point coordinates (columns 2 and 3) must be numeric",
sep = ""))
}
if (max(coords[, "decimallongitude"]) > 180) {
warning(paste("longitude values outside possible range in row:", rownames(coords[coords[,
"decimallongitude"] > 180, ]), ". ", "Row deleted", sep = ""))
coords <- coords[!coords[, "decimallongitude"] > 180, ]
}
if (min(coords[, "decimallongitude"]) < -180) {
warning(paste("longitude values outside possible range in row: ", rownames(coords[coords[,
"decimallongitude"] < -180, ]), ". ", "Row deleted", sep = ""))
coords <- coords[!coords[, "decimallongitude"] < -180, ]
}
if (max(coords[, "decimallatitude"]) > 90) {
warning(paste("latitude values outside possible range in row:", rownames(coords[coords[,
"decimallatitude"] > 90, ]), ". ", "Row deleted", sep = ""))
coords <- coords[!coords[, "decimallatitude"] > 90, ]
}
if (min(coords[, "decimallatitude"]) < -90) {
warning(paste("latitude values outside possible range in row:", rownames(coords[coords[,
"decimallatitude"] < -90, ]), ". ", "Row deleted", sep = ""))
coords <- coords[!coords[, "decimallatitude"] < -90, ]
}
if (!is.character(coords[, "species"]) && !is.factor(coords[, "species"])) {
warning("species name (column 1) should be a string or a factor")
}
coords[, 1] <- as.factor(coords[, "species"])
coordi <- coords[, c(2, 3)]
names(coordi) <- c("decimallongitude", "decimallatitude")
areanam <- areanames
res <- list(species = coords[, 1], species_coordinates = coordi, polygons = poly,
areanam = areanam)
class(res) <- "spgeoIN"
return(res)
}
.SpGeoCodH <- function(x, areanames = NULL, occ.thresh = 0) {
if (class(x$polygons) == "SpatialPolygonsDataFrame") {
if (!areanames %in% names(x$polygons@data)) {
stop(sprintf("column '%s' not found", areanames))
}
nam.test <- as.vector(unlist(x$polygons@data[, areanames]))
if ("NA" %in% nam.test) {
warning("the polygondata contain a polygon named NA. this can cause problems. Please rename")
}
}
# point in polygon test
kkk <- .PipSamp(x, columnname = areanames)
# number of records per species per polygon, everything belo occ.thresh is
# discarded
spsum <- .SpSumH(kkk, y = x$species, occ.thresh = occ.thresh)
if (length(spsum) == 0) {
namco <- c("species", names(x$polygons))
fill <- matrix(0, nrow = length(unique(kkk)), ncol = length(names(x$polygons)))
fill <- data.frame(fill)
spsum <- data.frame(cbind(as.character(unique(kkk)), fill))
names(spsum) <- namco
}
# species number per polygon
sppol <- spsum
sppol[sppol > 1] <- 1
sppol <- colSums(sppol)
# SpatialPolygonsDataFrame with species number per polygon based on areanam
pol.df <- as(x$polygons, "data.frame")
nums <- data.frame(sppol)
pol.df.m <- merge(pol.df, nums, sort = FALSE, by.x = areanames, by.y = "row.names",
all.x = TRUE)
pol.df.m <- pol.df.m[match(pol.df[, areanames], pol.df.m[, areanames]),
]
rownames(pol.df.m) <- rownames(pol.df)
pol.df.m[is.na(pol.df.m$sppol), "sppol"] <- 0
pol <- SpatialPolygonsDataFrame(as(x$polygons, "SpatialPolygons"), data = pol.df.m)
# create output
out <- list(samples = data.frame(species = x$species, decimallongitude = x$species_coordinates[,
1], decimallatitude = x$species_coordinates[, 2], homepolygon = kkk),
polygons = pol, spec_table = spsum, polygon_table = sppol, areanam = areanames)
class(out) <- "spgeoOUT"
return(out)
}
.SpSumH <- function(x, y, occ.thresh = occ.thresh) {
liste <- as.character(na.omit(unique(x)))
if (length(liste) == 0) {
spec_sum <- NULL
} else {
dat <- data.frame(x, y)
spec_sum <- as.data.frame.matrix(t(table(dat)))
if (occ.thresh > 0) {
filtperc <- apply(spec_sum, 2, function(k) {
k/rowSums(spec_sum) * 100
})
spec_sum[filtperc < occ.thresh] <- 0
}
}
return(spec_sum)
}
.WriteTablesSpGeo <- function(x, path, prefix = "", verbose = FALSE) {
if (missing(path)) {
path <- getwd()
}
write.table(x$samples, file = file.path(path, paste(prefix, "sample_classification_to_polygon.txt",
sep = "")), row.names = FALSE, sep = "\t")
write.table(x$spec_table, file = file.path(path, paste(prefix, "species_occurences_per_polygon.txt",
sep = "")), row.names = TRUE, sep = "\t")
write.table(x$polygon_table, file = file.path(path, paste(prefix, "speciesnumber_per_polygon.txt",
sep = "")), row.names = TRUE, sep = "\t", col.names = FALSE)
write.table(x$samples[x$samples$homepolygon == "not_classified", ], file = file.path(path,
paste(prefix, "unclassified samples.txt", sep = "")), row.names = FALSE,
sep = "\t")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.