-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathred.R
1743 lines (1622 loc) · 90.6 KB
/
red.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#####RED - IUCN Redlisting Tools
#####Version 1.5.0 (2020-05-04)
#####By Pedro Cardoso
#####Maintainer: pedro.cardoso@helsinki.fi
#####Reference: Cardoso, P.(2017) An R package to facilitate species red list assessments according to the IUCN criteria. Biodiversity Data Journal 5: e20530 doi: 10.3897/BDJ.5.e20530
#####Changed from v1.4.0:
#####added function rli.predict to interpolate and extrapolate linearly beyond the years assessed
#####added new options in functions rli and rli.multi on how to deal with DD species when bootstrapping
#####required packages
library("BAT")
library("dismo")
library("gdistance")
library("geosphere")
library("graphics")
library("grDevices")
library("jsonlite")
library("maptools")
library("methods")
library("raster")
library("rgdal")
library("rgeos")
library("sp")
library("stats")
library("utils")
#' @import gdistance
#' @import graphics
#' @import jsonlite
#' @import maptools
#' @import rgdal
#' @import rgeos
#' @import sp
#' @import stats
#' @import utils
#' @importFrom BAT contribution
#' @importFrom geosphere areaPolygon
#' @importFrom grDevices chull dev.copy dev.off pdf
#' @importFrom methods slot
#' @importFrom raster area cellStats clump crop extent extract getValues layerStats mask raster rasterize rasterToPoints rasterToPolygons reclassify res sampleRandom scalebar terrain trim writeRaster xmax xmin
raster::rasterOptions(maxmemory = 2e+09)
globalVariables(c("worldborders"))
###############################################################################
##############################AUX FUNCTIONS####################################
###############################################################################
longlat2utm <- function(longlat){
longlat = as.matrix(longlat)
minlong = min(longlat[,1])
zone = floor((minlong + 180) / 6) + 1
res = rgdal::project(longlat, paste("+proj=utm +zone=",zone," ellps=WGS84",sep=''))
return(res)
}
utm2longlat <- function(utm, zone){
if(class(utm) == "RasterLayer"){
if(!is.null(zone))
raster::crs(utm) <- paste("+proj=utm +zone=", zone, sep="")
res <- raster::projectRaster(utm, crs = "+proj=longlat +datum=WGS84", method='ngb')
} else {
utm <- SpatialPoints(utm, CRS(paste("+proj=utm +zone=", zone,sep="")))
res <- as.data.frame(spTransform(utm,CRS(paste("+proj=longlat"))))
}
return(res)
}
##warn if maxent.jar is not available
warnMaxent <- function(){
warning("RED could not find maxent.jar.
1. Download the latest version of maxent from:
https://biodiversityinformatics.amnh.org/open_source/maxent/
2. Move the file maxent.jar to the java directory inside dismo package
(there should be a file named dismo.jar already there)
3. Install the latest version of java runtime environment (JRE) with the same architecture (32 or 64 bits) as your version of R:
http://www.oracle.com/technetwork/java/javase/downloads/jre8-downloads-2133155.html")
}
##detect which layers are categorical by checking if all values are integers and if the max is less than 50 (may fail, just an attempt)
find.categorical <- function(layers){
categorical = c()
for(l in 1:(dim(layers)[3])){
lay <- raster::as.matrix(layers[[l]])
lay <- as.vector(lay)
lay <- lay[!is.na(lay)]
if(sum(floor(lay)) == sum(lay) && length(unique(lay)) < 50)
categorical = c(categorical, l)
}
return(categorical)
}
##basic function to calculate the rli of any group of species
rli.calc <- function(spData, tree = NULL, boot = FALSE, dd = FALSE, runs = 1000){
if(all(is.na(spData)))
return(NA)
spData <- rli.convert(spData) ##call function to convert spData to a 0-1 scale
if(is.null(tree)){ ##if not weighted by PD or FD
if(!boot){ ##if no bootstrap to be made
return (mean(spData, na.rm = TRUE))
} else {
run <- rep(NA, runs)
if(!dd){
for(i in 1:runs){
rnd <- sample(spData, replace = TRUE) ##bootstrap with all species
run[i] <- mean(rnd, na.rm = TRUE)
}
} else { ##bootstrap with only DD species
nDD = sum(is.na(spData)) ##number of DD species
rliBase = sum(spData, na.rm = TRUE)
for(i in 1:runs){
rnd <- sample(spData[!is.na(spData)], nDD, replace = TRUE)
run[i] <- (rliBase + sum(rnd)) / length(spData)
}
}
res <- matrix(quantile(run, c(0.025, 0.5, 0.975)), nrow = 1)
colnames(res) <- c("LowCL", "Median", "UpCL")
return(res)
}
} else { ##if weighted by PD or FD, still to work, not available at the moment!!!!!!!!!!!!!!!!!!!!!!!!!!!!
comm <- matrix(1, nrow = 2, ncol = length(spData))
contrib <- BAT::contribution(comm, tree, relative = TRUE)[1,]
contrib <- contrib/sum(contrib[!is.na(spData)]) #needed to standardize the contribution by the total contribution of species living in the community
if(!boot){ ##if no bootstrap to be made
return(sum(spData * contrib, na.rm = TRUE))
} else {
run <- rep(NA, runs)
for(i in 1:runs){
rndSpp <- sample(length(spData), replace = TRUE)
rndComm <- spData[rndSpp]
rndContrib <- contrib[rndSpp]/sum(contrib[rndSpp])
run[i] <- sum(rndComm * rndContrib, na.rm = TRUE)
}
res <- matrix(quantile(run, c(0.025, 0.5, 0.975)), nrow = 1)
colnames(res) <- c("LowCL", "Median", "UpCL")
return(res)
}
}
}
##function to convert strings to numbers in the RLI
rli.convert <- function(spData){
if(!is.numeric(spData)){ ##if letters are given, convert to [0,1]
spData <- replace(spData, which(spData == "EX" ), 0)
spData <- replace(spData, which(spData == "EW" ), 0)
spData <- replace(spData, which(spData == "RE" ), 0)
spData <- replace(spData, which(spData == "CR" ), 0.2)
spData <- replace(spData, which(spData == "CR(PE)" ), 0.2)
spData <- replace(spData, which(spData == "EN" ), 0.4)
spData <- replace(spData, which(spData == "VU" ), 0.6)
spData <- replace(spData, which(spData == "NT" ), 0.8)
spData <- replace(spData, which(spData == "LC" ), 1)
spData <- replace(spData, which(spData == "DD" ), NA)
spData <- as.numeric(spData)
} else if (all(spData == floor(spData))){ #if all integers, a scale [0,5] is given, convert to [0,1]
spData <- 1 - spData/5
}
return(spData)
}
##################################################################################
##################################MAIN FUNCTIONS##################################
##################################################################################
#' Setup GIS directory.
#' @description Setup directory where GIS files are stored.
#' @param gisPath Path to the directory where the gis files are stored.
#' @details Writes a txt file in the red directory allowing the package to always access the world GIS files directory.
#' @export
red.setDir <- function(gisPath = NULL){
if(is.null(gisPath))
gisPath <- readline("Input directory for storing world gis layers:")
gisPath <- paste(gisPath, "/", sep = "")
redFile <- paste(find.package("red"), "/red.txt", sep = "")
dput(gisPath, redFile)
}
#' Read GIS directory.
#' @description Read directory where GIS files are stored.
#' @details Reads a txt file pointing to where the world GIS files are stored.
#' @export
red.getDir <- function(){
redFile <- paste(find.package("red"), "/red.txt", sep = "")
if (file.exists(redFile)){ #if there is already a file read from it
dir <- dget(redFile)
} else {
warning(paste(redFile, "not found, please run red.setDir()"))
return()
}
return(dir)
}
#' Download and setup GIS files.
#' @description Setup red to work with species distribution modelling and layers available online.
#' @details Please check that you have at least 50Gb free in your disk (and a fast internet connection) to download all files. In the end of the process "only" 17.4Gb will be left though. This function will:
#' 1. Check if maxent.jar is available in the dismo package directory.
#' 2. Ask user input for GIS directory.
#' 3. Download global bioclim and elevation files (20) from http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_30s_bio.zip.
#' 4. Download landcover files (12) from http://data.earthenv.org/consensus_landcover/without_DISCover/.
#' 5. Unzip all files and delete the originals.
#' 6. Create a new layer (1) with the dominant land cover at each cell.
#' 7. Resample all files (33) to approximately 10x10km (for use with widespread species) grid cells.
#' Sit back and enjoy, this should take a while.
#' @export
red.setup <- function(){
##test if maxent.jar is in the right directory
if(!file.exists(paste(.libPaths()[[1]], "/dismo/java/maxent.jar", sep=""))){
warnMaxent()
return()
}
oldwd = getwd()
on.exit(expr = setwd(oldwd))
gisdir = red.setDir()
setwd(gisdir)
##basic setup
pb <- txtProgressBar(min = 0, max = 33, style = 3)
##download and process bioclim
download.file("http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_30s_bio.zip", "bioclim2.zip")
unzip(zipfile = "bioclim.zip")
file.remove("bioclim.zip")
for(i in 1:19){
setTxtProgressBar(pb, i)
if(i < 10)
rast <- raster(paste("wc2.0_bio_30s_0", i, ".tif", sep=""))
else
rast <- raster(paste("wc2.0_bio_30s_", i, ".tif", sep=""))
rast <- crop(rast, c(-180, 180, -56, 90))
writeRaster(rast, paste("red_1km_", i, ".tif", sep=""))
rast <- aggregate(rast, 10)
writeRaster(rast, paste("red_10km_", i, ".tif", sep=""))
if(i < 10)
file.remove(paste("wc2.0_bio_30s_0", i, ".tif", sep=""))
else
file.remove(paste("wc2.0_bio_30s_", i, ".tif", sep=""))
gc()
}
##download and process altitude
setTxtProgressBar(pb, 20)
download.file("http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/alt_30s_bil.zip", "alt_30s_bil.zip")
unzip(zipfile = "alt_30s_bil.zip")
file.remove("alt_30s_bil.zip")
rast <- raster("alt.bil")
rast <- crop(rast, c(-180, 180, -56, 90))
writeRaster(rast, "red_1km_20.tif")
rast <- aggregate(rast, 10)
writeRaster(rast, "red_10km_20.tif")
file.remove("alt.bil")
file.remove("alt.hdr")
gc()
##download and process land cover
altmask1 = raster("red_1km_20.tif")
altmask10 = raster("red_10km_20.tif")
for(i in 5:12){
setTxtProgressBar(pb, (i+20))
download.file(paste("http://data.earthenv.org/consensus_landcover/without_DISCover/Consensus_reduced_class_", i, ".tif", sep=""), destfile = paste("Consensus_reduced_class_", i, ".tif", sep=""), mode = "wb")
rast <- raster(paste("Consensus_reduced_class_", i, ".tif", sep=""))
rast <- mask(rast, altmask1)
writeRaster(rast, paste("red_1km_", (i+20), ".tif", sep=""))
rast <- aggregate(rast, 10)
#maskLayer <- sum(altmask, rast)
#maskLayer[!is.na(maskLayer)] <- 1
rast <- mask(rast, altmask10)
writeRaster(rast, paste("red_10km_", (i+20), ".tif", sep=""))
file.remove(paste("Consensus_reduced_class_", i, ".tif", sep=""))
gc()
}
remove(rast)
##create new rasters with most common landcover at each cell
setTxtProgressBar(pb, 33)
max1 <- raster()
max10 <- raster()
for(i in 21:32){
rast <- raster(paste("red_1km_", i, ".tif", sep=""))
max1 <- raster::stack(max1, rast)
rast <- raster(paste("red_10km_", i, ".tif", sep=""))
max10 <- raster::stack(max10, rast)
}
max1 <- which.max(max1)
writeRaster(max1, "red_1km_33.tif")
max10 <- which.max(max10)
writeRaster(max10, "red_10km_33.tif")
remove(max1, max10)
gc()
setwd(oldwd)
##Now the files should be named as:
##red_1km_1.tif
##...
##red_10km_33.tif
##Where 1 to 19 are the corresponding bioclim variables, 20 is altitude, 21 to 32 are landcover proportion and 33 is most common landcover per cell
#download country borders (not working Feb. 2017)
#download.file("http://biogeo.ucdavis.edu/data/gadm2.6/countries_gadm26.rds", destfile = paste("worldcountries.rds"), mode = "wb")
}
#' Download taxon records from GBIF.
#' @description Downloads species or higher taxon data from GBIF and outputs non-duplicate records with geographical coordinates.
#' @param taxon Taxon name.
#' @details As always when using data from multiple sources the user should be careful and check if records "make sense". This can be done by either ploting them in a map (e.g. using red::map.draw()) or using red::outliers().
#' @return A data.frame with longitude and latitude, plus species names if taxon is above species.
#' @examples records("Nephila senegalensis")
#' @export
records <- function(taxon){
taxon = unlist(strsplit(taxon, split = " ")[[1]])
dat <- dismo::gbif(taxon[1], paste(taxon[2], "*", sep = ""))
dat <- dat[c("species","lon","lat")] #filter columns
dat <- dat[!(is.na(dat$lon) | is.na(dat$lat)),] #filter rows
dat <- unique(dat) #delete duplicate rows
colnames(dat) <- c("Species", "long", "lat")
if (length(taxon) == 1){ #if genus
dat[which(is.na(dat[,1])),1] <- paste(taxon, "sp.")
} else { #if species
dat <- dat[,-1]
}
return(dat)
}
#' Move records to closest non-NA cell.
#' @description Identifies and moves presence records to cells with environmental values.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of species occurrence records.
#' @param layers Raster* object as defined by package raster.
#' @param buffer Maximum distance in map units that a record will move. If 0 all NA records will be changed.
#' @details Often records are in coastal or other areas for which no environmental data is available. This function moves such records to the closest cells with data so that no information is lost during modelling.
#' @return A matrix with new coordinate values.
#' @examples rast <- raster::raster(matrix(c(rep(NA,100), rep(1,100), rep(NA,100)), ncol = 15))
#' pts <- cbind(runif(100, 0, 0.55), runif(100, 0, 1))
#' raster::plot(rast)
#' points(pts)
#' pts <- move(pts, rast)
#' raster::plot(rast)
#' points(pts)
#' @export
move <- function(longlat, layers, buffer = 0){
layers <- layers[[1]]
values <- extract(layers, longlat) #get values of each record
suppressWarnings(
for(i in which(is.na(values))){ #if a value is NA, move it
distRaster = raster::distanceFromPoints(layers, longlat[i,])
distRaster = mask(distRaster, layers)
vmin = raster::minValue(distRaster)
if(buffer <= 0 || buffer > vmin){
vmin = rasterToPoints(distRaster, function(x) x == vmin)
longlat[i,] = vmin[1,1:2]
}
}
)
return(longlat)
}
#' Visual detection of outliers.
#' @description Draws plots of sites in geographical (longlat) and environmental (2-axis PCA) space.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of species occurrence records.
#' @param layers Raster* object as defined by package raster. It can be any set of environmental layers thought to allow the identification of environmental outliers.
#' @details Erroneous data sources or errors in transcriptions may introduce outliers that can be easily detected by looking at simple graphs of geographical or environmental space.
#' @return A data.frame with coordinate values and distance to centroid in pca is returned. Two plots are drawn for visual inspection. The environmental plot includes row numbers for easy identification of possible outliers.
#' @examples data(red.records)
#' data(red.layers)
#' outliers(red.records, red.layers[[1:3]])
#' @export
outliers <- function(longlat, layers){
if(dim(layers)[3] == 33) #if layers come from raster.read
pca <- raster.reduce(layers[[1:19]], n = 2)
else
pca <- raster.reduce(layers, n = 2)
##extract pca values from longlat
pca <- as.data.frame(raster::extract(pca, longlat))
goodRows <- which(!is.na(pca[,1]))
pca <- pca[goodRows,]
longlat <- longlat[goodRows,]
par(mfrow = c(1,2))
map.draw(longlat, layers[[1]], spName = "Geographical")
raster::plot(pca, main = "Environmental", type = "n")
centroid = colMeans(pca)
text(centroid[1], centroid[2], label = "X")
for(i in 1:nrow(pca)){
text(pca[i,1], pca[i,2], label = row.names(longlat)[i])
}
##build new matrix ordered by distance to centroid
dist2centroid = apply(pca, 1, function(x) dist(rbind(x, centroid)))
out = as.data.frame(cbind(longlat, dist2centroid))
out = out[order(-dist2centroid),]
return(out)
}
#' Spatial thinning of occurrence records.
#' @description Thinning of records with minimum distances either absolute or relative to the species range.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of species occurrence records.
#' @param distance Distance either in relative terms (proportion of maximum distance between any two records) or in raster units.
#' @param relative If TRUE, represents the proportion of maximum distance between any two records. If FALSE, is in raster units.
#' @param runs Number of runs
#' @details Clumped distribution records due to ease of accessibility of sites, emphasis of sampling on certain areas in the past, etc. may bias species distribution models.
#' The algorithm used here eliminates records closer than a given distance to any other record. The choice of records to eliminate is random, so a number of runs are made and the one keeping more of the original records is chosen.
#' @return A matrix of species occurrence records separated by at least the given distance.
#' @examples records <- matrix(sample(100), ncol = 2)
#' par(mfrow=c(1,2))
#' graphics::plot(records)
#' records <- thin(records, 0.1)
#' graphics::plot(records)
#' @export
thin <- function(longlat, distance = 0.01, relative = TRUE, runs = 100){
longlat = longlat[!duplicated(longlat),] #first, remove duplicate rows
nSites = nrow(longlat)
if(nSites < 4)
return(longlat)
##if relative, calculate maxDist between any two points
if(relative){
if(nSites < 40){ #if limited number of sites use all data
maxDist = 0
for(x in 1:(nSites-1)){
for(y in (x+1):nSites){
maxDist = max(maxDist,((longlat[x,1]-longlat[y,1])^2+(longlat[x,2]-longlat[y,2])^2)^.5)
}
}
} else { #if many sites use hypothenusa of square encompassing all of them
horiDist = max(longlat[,1]) - min(longlat[,1])
vertDist = max(longlat[,2]) - min(longlat[,2])
maxDist = (horiDist^2 + vertDist^2)^0.5
}
distance = maxDist*distance
}
listSites = matrix(longlat[1,], ncol=2, byrow = TRUE)
for (r in 1:runs){
longlat = longlat[sample(nSites),] ##shuffle rows (sites)
rndSites = longlat[1,] ##start with first random site
for(newSite in 2:nSites){
for(oldSite in 1:(newSite-1)){
addSite = TRUE
dist = ((longlat[newSite,1]-longlat[oldSite,1])^2+(longlat[newSite,2]-longlat[oldSite,2])^2)^.5
if(dist < distance){
addSite = FALSE
break
}
}
if(addSite)
rndSites = rbind(rndSites, longlat[newSite,])
}
if(nrow(rndSites) > nrow(listSites))
listSites = rndSites
}
return(as.matrix(listSites))
}
#' Read and buffer raster layers.
#' @description Read raster layers of environmental or other variables and crop them to a given extent around the known occurrences.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of species occurrence records.
#' @param layers Raster* object as defined by package raster.
#' @param ext Either extent of map or buffer around the known records used to crop layers. If buffer, it is relative to the maximum distance between any two records.
#' @details If layers are not given, the function will read either 30 arc-second (approx. 1km) or 5 arc-minutes (approx. 10km) resolution rasters from worldclim (Fick & Hijmans 2017) and landcover (Tuanmu & Jetz 2014) if red.setup() is run previously.
#' @return A RasterStack object (If no layers are given: Variables 1-19 = bioclim, 20 = elevation, 21-32 = proportion landcover, 33 = most common landcover).
#' @references Fick, S.E. & Hijmans, R.J. (2017) Worldclim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, in press.
#' @references Tuanmu, M.-N. & Jetz, W. (2014) A global 1-km consensus land-cover product for biodiversity and ecosystem modeling. Global Ecology and Biogeography, 23: 1031-1045.
#' @examples data(red.layers)
#' data(red.records)
#' par(mfrow=c(1,2))
#' raster::plot(red.layers[[1]])
#' points(red.records)
#' croppedLayers <- raster.read(red.records, red.layers, 0.1)
#' raster::plot(croppedLayers[[1]])
#' points(red.records)
#' @export
raster.read <- function(longlat, layers = NULL, ext = 1){
xmin = min(longlat[,1])
xmax = max(longlat[,1])
xlen = xmax - xmin
ymin = min(longlat[,2])
ymax = max(longlat[,2])
ylen = ymax - ymin
if(is.null(layers)){ ##if no layers are provided read the ones available
gisdir = red.getDir()
##calculate species range and buffer around it
if(eoo(longlat) < 200000){
layers <- raster::stack(raster::raster(paste(gisdir, "red_1km_1.tif", sep = "")))
for(i in 2:33)
layers <- raster::stack(layers, raster::raster(paste(gisdir, "red_1km_", i, ".tif", sep = "")))
} else {
layers <- raster::stack(raster::raster(paste(gisdir, "red_10km_1.tif", sep = "")))
for(i in 2:33)
layers <- raster::stack(layers, raster::raster(paste(gisdir, "red_10km_", i, ".tif", sep = "")))
}
##determine longitude limits of species to check if crop and paste are needed around longitude 180 for Pacific species
if(xmin < -90 && xmax > 90 && sum(longlat[longlat[,1] < 90 && longlat[,1] > -90,]) != 0){
##crop and merge layers
rightHalf = crop(layers, c(0,180,raster::extent(layers)@ymin,raster::extent(layers)@ymax))
raster::extent(rightHalf) <- c(-180,0,raster::extent(layers)@ymin,raster::extent(layers)@ymax)
leftHalf = crop(layers, c(-180,0,raster::extent(layers)@ymin,raster::extent(layers)@ymax))
raster::extent(leftHalf) <- c(0,180,raster::extent(layers)@ymin,raster::extent(layers)@ymax)
layers <- merge(rightHalf, leftHalf)
##modify longlat
for(i in 1:nrow(longlat))
if(longlat[i,1] > 0)
longlat[i,1] = longlat[i,1] - 180
else
longlat[i,1] = longlat[i,1] + 180
}
}
if(length(ext) == 4) ##if absolute extent is given crop and return, else calculate buffer
return(crop(layers, ext))
if(xlen == 0) ##in case some dimensions are inexistent consider equal to extent
xlen = ext
if(ylen == 0)
ylen = ext
##calculate new extent of layers and crop
ext = max(1, ((xlen + ylen) * ext))
xmin <- max(raster::extent(layers)@xmin, xmin-ext)
xmax <- min(raster::extent(layers)@xmax, xmax+ext)
ymin <- max(raster::extent(layers)@ymin, ymin-ext)
ymax <- min(raster::extent(layers)@ymax, ymax+ext)
layers <- crop(layers, c(xmin,xmax,ymin,ymax))
return(layers)
}
#' Uniformize raster layers.
#' @description Crop raster layers to minimum size possible and uniformize NA values across layers.
#' @param layers Raster* object as defined by package raster.
#' @details Excludes all marginal rows and columns with only NA values and change values to NA if they are NA in any of the layers.
#' @return A Raster* object, same class as layers.
#' @examples data(red.layers)
#' raster::plot(raster.clean(red.layers))
#' @export
raster.clean <- function(layers){
##apply mask to have NAs everywhere where any layer has NAs
maskLayer <- sum(layers)
maskLayer[!is.na(maskLayer)] <- 1
layers <- mask(layers, maskLayer)
##crop by excluding external rows and columns with NAs only
layers <- trim(layers)
return(layers)
}
#' Reduce dimensionality of raster layers.
#' @description Reduce the number of layers by either performing a PCA on them or by eliminating highly correlated ones.
#' @param layers Raster* object as defined by package raster.
#' @param method Either Principal Components Analysis ("pca", default) or Pearson's correlation ("cor").
#' @param n Number of layers to reduce to.
#' @param thres Value for pairwise Pearson's correlation above which one of the layers (randomly selected) is eliminated.
#' @details Using a large number of explanatory variables in models with few records may lead to overfitting. This function allows to avoid it as much as possible.
#' If both n and thres are given, n has priority. If method is not recognized and layers come from raster.read function, only landcover is reduced by using only the dominating landuse of each cell.
#' @return A RasterStack object.
#' @export
raster.reduce <- function(layers, method = "pca", n = NULL, thres = NULL){
##method = "pca, cor", if unrecognized method only reduce landcover but not climate
out <- raster::stack()
if(dim(layers)[3] == 33){ ##check if layers are obtained with raster.read
out <- raster::stack(layers[[33]])
layers = layers[[1:19]]
}
if(method == "cor"){ ##if correlation
if(is.null(n)){
if(is.null(thres))
thres = 0.7
for(i in 1:dim(layers)[3]){ ##delete layers until none are correlated above threshold
cor = as.matrix(as.dist(layerStats(layers, 'pearson', na.rm = TRUE)[[1]]))
if(max(cor) < thres)
break
corLayer = sample(which(cor == max(cor), arr.ind = TRUE)[,1],1)
layers = layers[[-corLayer]]
}
} else {
while (dim(layers)[3] > n){ ##delete layers until reaching n layers
cor = abs(as.matrix(as.dist(layerStats(layers, 'pearson', na.rm = TRUE)[[1]])))
corLayer = sample(which(cor == max(cor), arr.ind = TRUE)[,1],1)
layers = layers[[-corLayer]]
}
}
} else if(method == "pca"){ ##if pca
if(is.null(n))
n = 3
if(sum(!is.na(getValues(layers[[1]]))) > 2000)
sr <- sampleRandom(layers, 1000)
else
sr <- sampleRandom(layers, as.integer(sum(!is.na(getValues(layers[[1]])))/2))
pca <- prcomp(sr)
layers <- raster::predict(layers, pca, index = 1:n)
for(i in 1:n)
names(layers[[i]]) <- paste("pca",i)
}
out <- raster::stack(layers, out)
return(out)
}
#' Create distance layer.
#' @description Creates a layer depicting distances to records using the minimum, average, distance to the minimum convex polygon or distance taking into account a cost surface.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of species occurrence records.
#' @param layers Raster* object as defined by package raster to serve as model to create distance layer. Cost surface in case of param ="cost".
#' @param type text string indicating whether the output should be the "minimum", "average", "mcp" or "cost" distance to all records. "mcp" means the distance to the minimum convex polygon encompassing all records.
#' @details Using distance to records in models may help limiting the extrapolation of the predicted area much beyond known areas.
#' @return A RasterLayer object.
#' @examples data(red.layers)
#' alt = red.layers[[3]]
#' data(red.records)
#' par(mfrow=c(3,2))
#' raster::plot(alt)
#' points(red.records)
#' raster::plot(raster.distance(red.records, alt))
#' raster::plot(raster.distance(red.records, alt, type = "average"))
#' raster::plot(raster.distance(red.records, alt, type = "mcp"))
#' raster::plot(raster.distance(red.records, alt, type = "cost"))
#' @export
raster.distance <- function(longlat, layers, type = "minimum"){
if(dim(layers)[3] > 1)
layers <- layers[[1]]
layers[!is.na(layers)] <- 0
if(type == "average"){
for(d in 1:nrow(longlat)){
layers <- layers + raster::distanceFromPoints(layers, longlat[d,])
}
layers <- layers/nrow(longlat)
names(layers) <- "average distance"
} else if (type == "mcp"){
vertices <- chull(longlat)
vertices <- c(vertices, vertices[1])
vertices <- longlat[vertices,]
poly = Polygon(vertices)
poly = Polygons(list(poly),1)
poly = SpatialPolygons(list(poly)) ##minimum convex polygon
longlat = rasterToPoints(rasterize(poly, layers))[,1:2]
layers <- mask(raster::distanceFromPoints(layers, longlat), layers)
names(layers) <- "mcp distance"
} else if (type == "cost"){
layers <- transition(layers, function(x) 1/mean(x), 8)
layers <- geoCorrection(layers)
layers <- accCost(layers, as.matrix(longlat))
names(layers) <- "cost distance"
} else {
layers <- mask(raster::distanceFromPoints(layers, longlat), layers)
names(layers) <- "minimum distance"
}
return(layers)
}
#' Create longitude layer.
#' @description Create a layer depicting longitude based on any other.
#' @param layers Raster* object as defined by package raster.
#' @details Using longitude (and latitude) in models may help limiting the extrapolation of the predicted area much beyond known areas.
#' @return A RasterLayer object.
#' @examples data(red.layers)
#' raster::plot(raster.long(red.layers))
#' @export
raster.long <- function(layers){
if(dim(layers)[3] > 1)
layers <- layers[[3]]
x <- rasterToPoints(layers)[,1:2]
long <- rasterize(x, layers, x[,1])
long <- mask(long, layers)
names(long) <- "longitude"
return(long)
}
#' Create latitude layer.
#' @description Create a layer depicting latitude based on any other.
#' @param layers Raster* object as defined by package raster.
#' @details Using latitude (and longitude) in models may help limiting the extrapolation of the predicted area much beyond known areas.
#' @return A RasterLayer object.
#' @examples data(red.layers)
#' raster::plot(raster.lat(red.layers[[1]]))
#' @export
raster.lat <- function(layers){
if(dim(layers)[3] > 1)
layers <- layers[[3]]
x <- rasterToPoints(layers)[,1:2]
lat <- rasterize(x, layers, x[,2])
lat <- mask(lat, layers)
names(lat) <- "latitude"
return(lat)
}
#' Create eastness layer.
#' @description Create a layer depicting eastness based on an elevation layer.
#' @param dem RasterLayer object of elevation (a digital elevation model - DEM) as defined by package raster.
#' @details Using elevation, aspect can be calculated. Yet, it is a circular variable (0 = 360) and has to be converted to northness and eastness to be useful for modelling.
#' @return A RasterLayer object.
#' @examples data(red.layers)
#' raster::plot(raster.east(red.layers[[3]]))
#' @export
raster.east <- function(dem){
asp <- terrain(dem, opt = "aspect")
return(sin(asp))
}
#' Create northness layer.
#' @description Create a layer depicting northness based on an elevation layer.
#' @param dem RasterLayer object of elevation (a digital elevation model - DEM) as defined by package raster.
#' @details Using elevation, aspect can be calculated. Yet, it is a circular variable (0 = 360) and has to be converted to northness and eastness to be useful for modelling.
#' @return A RasterLayer object.
#' @examples data(red.layers)
#' raster::plot(raster.north(red.layers[[3]]))
#' @export
raster.north <- function(dem){
asp <- terrain(dem, opt = "aspect")
return(cos(asp))
}
#' Predict species distribution.
#' @description Prediction of potential species distributions using maximum entropy (maxent).
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of each occurrence record.
#' @param layers Predictor variables, a Raster* object as defined by package raster.
#' @param error Vector of spatial error in longlat (one element per row of longlat) in the same unit as longlat. Used to move any point randomly within the error radius.
#' @param year Vector of sampling years in longlat (one element per row of longlat). Used to exclude old records with a given probability proportional to time passed since sampling (never excluded only for current year).
#' @param idconf Vector of identification confidence in longlat (one element per row of longlat). Used to exclude uncertain records with a given probability. Can be on any scale where max values are certain (e.g. from 1 - very uncertain to 10 - holotype).
#' @param categorical Vector of layer indices of categorical (as opposed to quantitative) data. If NULL the package will try to find them automatically based on the data.
#' @param thres Threshold of logistic output used for conversion of probabilistic to binary (presence/absence) maps. If 0 this will be the value that maximizes the sum of sensitivity and specificity.
#' @param testpercentage Percentage of records used for testing only. If 0 all records will be used for both training and testing.
#' @param mcp Used for a precautionary approach. If TRUE, all areas predicted as present but outside the minimum convex hull polygon encompassing all occurrence records are converted to absence. Exceptions are cells connected to other areas inside the polygon.
#' @param points If TRUE, force map to include cells with presence records even if suitable habitat was not identified.
#' @param eval If TRUE, build a matrix with AUC, Kappa, TSS, EOO (from raw data), EOO (from model), AOO (from raw data) and AOO (from model).
#' @param runs If <= 0 no ensemble modelling is performed. If > 0, ensemble modelling with n runs is made. For each run, a new random sample of occurrence records (if testpercentage > 0), background points and predictive variables (if subset > 0) are chosen. In the ensemble model, each run is weighted as max(0, (runAUC - 0.5)) ^ 2.
#' @param subset Number of predictive variables to be randomly selected from layers for each run if runs > 0. If <= 0 all layers are used on all runs. Using a small number of layers is usually better than using many variables for rare species, with few occurrence records (Lomba et al. 2010, Breiner et al. 2015).
#' @details Builds maxent (maximum entropy) species distribution models (Phillips et al. 2004, 2006; Elith et al. 2011) using function maxent from R package dismo (Hijmans et al. 2017). Dismo requires the MaxEnt species distribution model software, a java program that can be downloaded from http://biodiversityinformatics.amnh.org/open_source/maxent. Copy the file 'maxent.jar' into the 'java' folder of the dismo package. That is the folder returned by system.file("java", package="dismo"). You need MaxEnt version 3.3.3b or higher. Please note that this program (maxent.jar) cannot be redistributed or used for commercial or for-profit purposes.
#' @return List with either one or two raster objects (depending if ensemble modelling is performed, in which case the second is a probabilistic map from all the runs) and, if eval = TRUE, a matrix with AUC, Kappa, TSS, EOO (from raw data), EOO (from model), AOO (from raw data) and AOO (from model). Aggregate values are taken from maps after transformation of probabilities to incidence, with presence predicted for cells with ensemble values > 0.5.
#' @references Breiner, F.T., Guisan, A., Bergamini, A., Nobis, M.P. (2015) Overcoming limitations of modelling rare species by using ensembles of small models. Methods in Ecology and Evolution, 6: 1210-1218.
#' @references Hijmans, R.J., Phillips, S., Leathwick, J., Elith, J. (2017) dismo: Species Distribution Modeling. R package version 1.1-4. https://CRAN.R-project.org/package=dismo
#' @references Lomba, A., Pellissier, L., Randin, C.F., Vicente, J., Moreira, F., Honrado, J., Guisan, A. (2010) Overcoming the rare species modelling paradox: a novel hierarchical framework applied to an Iberian endemic plant. Biological Conservation, 143: 2647-2657.
#' @references Phillips, S.J., Dudik, M., Schapire, R.E. (2004) A maximum entropy approach to species distribution modeling. Proceedings of the Twenty-First International Conference on Machine Learning. p. 655-662.
#' @references Phillips, S.J., Anderson, R.P., Schapire, R.E. (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190: 231-259.
#' @references Elith, J., Phillips, S.J., Hastie, T., Dudik, M., Chee, Y.E., Yates, C.J. (2011) A statistical explanation of MaxEnt for ecologists. Diversity and Distributions, 17: 43-57.
#' @export
map.sdm <- function(longlat, layers, error = NULL, year = NULL, idconf = NULL, categorical = NULL, thres = 0, testpercentage = 0, mcp = TRUE, points = FALSE, eval = TRUE, runs = 0, subset = 0){
raster::rasterOptions(maxmemory = 2e+09)
origLonglat = longlat
##if ensemble is to be done
if(runs > 0){
longlat = origLonglat
#if there is spatial error randomly move points within its radius
if(!is.null(error)){
for(i in 1:nrow(longlat)){
#move up to given error (angular movement converted to x and y)
rndAngle = sample(1:360, 1)
rndDist = runif(1, 0, error[i])
longlat[i,1] = longlat[i,1] + rndDist * cos(rndAngle)
longlat[i,2] = longlat[i,2] + rndDist * sin(rndAngle)
}
}
#if there is year
if(!is.null(year)){
for(i in 1:nrow(longlat)){
if(year[i] < sample(min(year):as.integer(substr(Sys.Date(), 1, 4)), 1))
longlat = longlat[-i,]
}
}
#if there is idconf
if(!is.null(idconf)){
for(i in 1:nrow(longlat)){
if(idconf[i] < sample(1:max(idconf), 1))
longlat = longlat[-i,]
}
}
if(eval)
runEval = matrix(NA, nrow = 1, ncol = 7)
runMap <- rasterize(longlat, layers[[1]], field = 0, background = 0)
pb <- txtProgressBar(min = 0, max = runs, style = 3)
totalAUC = 0
for(i in 1:runs){
if(subset > 0 && subset < dim(layers)[3]){
runLayers <- layers[[sample.int(dim(layers)[3], subset)]]
thisRun <- map.sdm(longlat, runLayers, error = NULL, year = NULL, idconf = NULL, categorical, thres, testpercentage, mcp, points, eval, runs = 0, subset = 0)
} else {
thisRun <- map.sdm(longlat, layers, error = NULL, year = NULL, idconf = NULL, categorical, thres, testpercentage, mcp, points, eval, runs = 0, subset = 0)
}
runAUC = 1
if(eval){
runAUC <- thisRun[[2]][1]
runAUC <- max(0, (runAUC - 0.5)) ^ 2 #weight the map by its AUC above 0.5 to the square
runEval <- rbind(runEval, thisRun[[2]])
thisRun <- thisRun[[1]]
}
totalAUC = totalAUC + runAUC
runMap <- runMap + (thisRun * runAUC)
setTxtProgressBar(pb, i)
}
runMap <- raster::calc(runMap, function(x) {x/totalAUC})
upMap <- reclassify(runMap, matrix(c(0,0.025,0,0.025,1,1), ncol = 3, byrow = TRUE))
consensusMap <- reclassify(runMap, matrix(c(0,0.499,0,0.499,1,1), ncol = 3, byrow = TRUE))
downMap <- reclassify(runMap, matrix(c(0,0.975,0,0.975,1,1), ncol = 3, byrow = TRUE))
if(mcp && aoo(consensusMap) >= 4)
consensusMap <- map.habitat(longlat, consensusMap, mcp = TRUE, eval = FALSE)
if(eval){
runEval <- runEval[-1,]
clEval <- matrix(NA, nrow = 3, ncol = 7)
colnames(clEval) <- c("AUC", "Kappa", "TSS", "EOO (raw)", "EOO (model)", "AOO (raw)", "AOO (model)")
rownames(clEval) <- c("UpCL", "Consensus", "LowCL")
clEval[1,] <- apply(runEval, 2, quantile, probs= 0.975, na.rm = TRUE)
clEval[2,] <- apply(runEval, 2, quantile, probs= 0.5, na.rm = TRUE)
clEval[3,] <- apply(runEval, 2, quantile, probs= 0.025, na.rm = TRUE)
clEval[1:3,4] <- eoo(longlat)
clEval[1:3,6] <- aoo(longlat)
clEval[1,5] <- eoo(upMap)
clEval[1,7] <- aoo(upMap)
clEval[2,5] <- eoo(consensusMap)
clEval[2,7] <- aoo(consensusMap)
clEval[3,5] <- eoo(downMap)
clEval[3,7] <- aoo(downMap)
return(list(consensusMap, runMap, clEval))
} else {
return (consensusMap)
}
}
longlat <- move(longlat, layers) #move all records falling on NAs
nPoints = min(1000, sum(!is.na(as.vector(layers[[1]])), na.rm = TRUE)/4)
bg <- dismo::randomPoints(layers, nPoints) ##extract background points
##if no categorical variables are given try to figure out which are
if(is.null(categorical))
categorical <- find.categorical(layers)
llTrain <- longlat
llTest <- longlat
if(testpercentage > 0){
testRecords <- sample(1:nrow(longlat), ceiling(nrow(longlat)*testpercentage/100))
llTrain <- longlat[-testRecords,]
llTest <- longlat[testRecords,]
}
mod <- dismo::maxent(layers, llTrain, a = bg, factors = categorical) ##build model
p <- raster::predict(mod, layers) ##do prediction
e <- dismo::evaluate(p = llTrain, a = bg, model = mod, x = layers) ##do evaluation of model
if(thres == 0)
thres <- dismo::threshold(e)$spec_sens ##extract threshold from evaluation
p <- reclassify(p, matrix(c(0,thres,0,thres,1,1), nrow=2, byrow = TRUE)) ##convert to presence/absence
if(mcp && aoo(p) >= 4)
p <- map.habitat(longlat, p, mcp = TRUE, eval = FALSE)
if(points)
p <- max(p, map.points(longlat, p, eval = FALSE))
if(eval){
e <- dismo::evaluate(p = llTest, a = bg, model = mod, x = layers, tr = thres) ##do evaluation of model with threshold
auc <- e@auc
kappa <- e@kappa
sensitivity <- as.numeric(e@TPR/(e@TPR+e@FNR))
specificity <- as.numeric(e@TNR/(e@TNR+e@FPR))
tss <- sensitivity + specificity - 1
eooRaw <- eoo(longlat)
aooRaw <- aoo(longlat)
aooModel <- aoo(p)
if(aooModel > 8)
eooModel <- eoo(p)
else
eooModel = aooModel
txtEval <- matrix(c(auc, kappa, tss, eooRaw, eooModel, aooRaw, aooModel), nrow = 1)
colnames(txtEval) <- c("AUC", "Kappa", "TSS", "EOO (raw)", "EOO (model)", "AOO (raw)", "AOO (model)")
return(list(p, txtEval))
} else {
return(p)
}
}
#' Map species distribution of habitat specialist.
#' @description Mapping of all habitat patches where the species is known to occur.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of each occurrence record.
#' @param layer RasterLayer object representing the presence/absence (1/0) of a single habitat type.
#' @param move If TRUE, identifies and moves presence records to closest cells with suitable habitat. Use when spatial error might put records outside the correct patch.
#' @param mcp If TRUE, all habitat patches inside the minimum convex hull polygon encompassing all occurrence records are converted to presence.
#' @param points If TRUE, force map to include cells with presence records even if suitable habitat was not identified.
#' @param eval If TRUE, build a matrix with EOO (from raw data), EOO (from model), AOO (from raw data) and AOO (from model).
#' @details In many cases a species has a very restricted habitat and we generally know where it occurs. In such cases using the distribution of the known habitat patches may be enough to map the species.
#' @return One raster object and, if eval = TRUE, a matrix with EOO (from raw data), EOO (from model), AOO (from raw data) and AOO (from model).
#' @export
map.habitat <- function(longlat, layer, move = TRUE, mcp = FALSE, points = FALSE, eval = TRUE){
if(points)
layer <- max(layer, map.points(longlat, layer, eval = FALSE))
if(move){
moveLayer <- layer
moveLayer[moveLayer == 0] <- NA
longlat <- move(longlat, moveLayer)
remove(moveLayer)
}
if(mcp){
vertices <- chull(longlat)
vertices <- c(vertices, vertices[1])
vertices <- longlat[vertices,]
poly = Polygon(vertices)
poly = Polygons(list(poly),1)
poly = SpatialPolygons(list(poly)) ##minimum convex polygon
patches <- raster::clump(layer, gaps=FALSE) ##individual patches, numbered
selPatches <- raster::unique(extract(patches, poly, df = TRUE, weights = TRUE)$clumps) ##which patches are inside polygon
} else {
patches <- raster::clump(layer, gaps=FALSE) ##individual patches, numbered
selPatches <- raster::unique(extract(patches, longlat, df = TRUE, weights = TRUE)$clumps) ##which patches have the species
}
selPatches <- selPatches[!is.na(selPatches)]
allPatches <- raster::unique(patches)
allPatches <- as.data.frame(cbind(allPatches, rep(0, length(allPatches))))
colnames(allPatches) <- c("patches", "selected")
allPatches[selPatches, 2] <- 1
patches <- raster::subs(patches, allPatches)
layer <- mask(layer, patches, maskvalue = 0, updatevalue = 0)
if(eval){
eooRaw <- eoo(longlat)
eooModel <- eoo(layer)
aooRaw <- aoo(longlat)
aooModel <- aoo(layer)
txtEval <- matrix(c(eooRaw, eooModel, aooRaw, aooModel), nrow = 1)
colnames(txtEval) <- c("EOO (raw)", "EOO (model)", "AOO (raw)", "AOO (model)")
return(list(layer, txtEval))
} else {
return(layer)
}
}
#' Map recorded distribution of species.
#' @description Mapping of all cells where the species is known to occur.
#' @param longlat Matrix of longitude and latitude or eastness and northness (two columns in this order) of each occurrence record.
#' @param layers Raster* object as defined by package raster. Any raster with the relevant extent and cell size can be used.
#' @param eval If TRUE, build a matrix with EOO and AOO calculated from occurrence records only.
#' @details To be used if either information on the species is very scarce (and it is not possible to model the species distribution) or, on the contrary, complete (and there is no need to model the distribution).
#' @return One raster object and, if EVAL = TRUE, a matrix with EOO and AOO.
#' @examples
#' data(red.records)
#' data(red.layers)
#' raster::plot(map.points(red.records, red.layers, eval = FALSE))
#' points(red.records)
#' @export
map.points <- function(longlat, layers, eval = TRUE){
p <- rasterize(longlat, layers[[1]], field = 1, background = 0)
maskLayer <- sum(layers)
maskLayer[!is.na(maskLayer)] <- 1
p <- mask(p, maskLayer)
if(eval){
eooRaw <- eoo(longlat)
aooRaw <- aoo(longlat)
txtEval <- matrix(c(eooRaw, aooRaw), nrow = 1)
colnames(txtEval) <- c("EOO", "AOO")
return(list(p, txtEval))
} else {
return(p)
}
}
#' Species distributions made easy (multiple species).
#' @description Single step for prediction of multiple species distributions. Output of maps (in pdf format), klms (for Google Earth) and relevant data (in csv format).
#' @param longlat data.frame of taxon names, longitude and latitude or eastness and northness (three columns in this order) of each occurrence record.
#' @param layers If NULL analyses are done with environmental layers read from data files of red.setup(). If a Raster* object as defined by package raster, analyses use these.
#' @param habitat Raster* object as defined by package raster. Habitat extent layer (0/1) used instead of layers if any species is an habitat specialist.
#' @param zone UTM zone if data is in metric units. Used only for correct placement of kmls and countries.
#' @param thin boolean defining if species data should be thinned before modeling (only for SDMs).
#' @param error Vector of spatial error in longlat (one element per row of longlat) in the same unit as longlat. Used to move any point randomly within the error radius.
#' @param move If TRUE, identifies and moves presence records to closest cells with environmental data. Use when spatial error might put records outside such data.
#' @param dem RasterLayer object. It should be a digital elevation model for calculation of elevation limits of the species. If NULL, dem from red.setup() is used if possible, otherwise it will be 0.
#' @param pca Number of pca axes for environmental data reduction. If 0 (default) no pca is made.
#' @param filename Name of output csv file with all results. If NULL it is named "Results_All.csv".
#' @param mapoption Vector of values within options: points, habitat and sdm; each value corresponding to the function to be used for each species (map.points, map.habitat, map.sdm). If a single value, all species will be modelled according to it. If NULL, the function will perform analyses using map.points. Species values must be in same order as latlong.
#' @param testpercentage Percentage of records used for testing only. If 0 all records will be used for both training and testing.
#' @param mintest Minimim number of total occurrence records of any species to set aside a test set. Only used if testpercentage > 0.
#' @param points If TRUE, force map to include cells with presence records even if suitable habitat was not identified.
#' @param runs If <= 0 no ensemble modelling is performed. If > 0, ensemble modelling with n runs is made. For each run, a new random sample of occurrence records (if testpercentage > 0), background points and predictive variables (if subset > 0) are chosen. In the ensemble model, each run is weighted as max(0, (runAUC - 0.5)) ^ 2.
#' @param subset Number of predictive variables to be randomly selected from layers for each run if runs > 0. If <= 0 all layers are used on all runs. Using a small number of layers is usually better than using many variables for rare species, with few occurrence records (Lomba et al. 2010, Breiner et al. 2015).
#' @return Outputs maps in asc, pdf and kml format, plus a file with EOO, AOO and a list of countries where the species is predicted to be present if possible to extract.
#' @references Breiner, F.T., Guisan, A., Bergamini, A., Nobis, M.P. (2015) Overcoming limitations of modelling rare species by using ensembles of small models. Methods in Ecology and Evolution, 6: 1210-1218.
#' @references Lomba, A., Pellissier, L., Randin, C.F., Vicente, J., Moreira, F., Honrado, J., Guisan, A. (2010) Overcoming the rare species modelling paradox: a novel hierarchical framework applied to an Iberian endemic plant. Biological Conservation, 143: 2647-2657.
#' @export
map.easy <- function(longlat, layers = NULL, habitat = NULL, zone = NULL, thin = TRUE, error = NULL, move = TRUE, dem = NULL, pca = 0, filename = NULL, mapoption = NULL, testpercentage = 0, mintest = 20, points = FALSE, runs = 0, subset = 0){
try(dev.off(), silent = TRUE)
spNames <- unique(longlat[,1])
nSp <- length(spNames)
if(is.null(mapoption))
mapoption = rep("points", nSp)
else if(length(mapoption) == 1)
mapoption = rep(mapoption, nSp)
else if(length(mapoption) != nSp)
return(warning("Number of species different from length of mapoption"))
if("sdm" %in% mapoption){
if(!file.exists(paste(.libPaths()[[1]], "/dismo/java/maxent.jar", sep=""))){
warnMaxent()