Mercurial > repos > ecology > obis_data

comparison obisindicators.r @ 0:ed3688d04c12 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/obisindicators commit 13ac67c0a21d742b29e6273cdff058560abad770
author ecology
date Tue, 05 Nov 2024 14:17:28 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:ed3688d04c12
1 #Rscript
2
3 ###########################################
4 ## Mapping alpha and beta diversity ##
5 ###########################################
6
7 #####Packages : obisindicators
8 # dplyr
9 # sf
10 # ggplot2
11 # rnaturalearth
12 # rnaturalearthdata
13 # viridis
14 # dggridr
15 library(magrittr)
16
17 ## remotes::install_github("r-barnes/dggridR")
18 #####Load arguments
19
20 args <- commandArgs(trailingOnly = TRUE)
21
22 # url for the S2 subset
23
24 if (length(args) < 4) {
25 stop("This tool needs at least 4 argument : longitude, latitude, species and number of records")
26 }else {
27 raster <- args[1]
28 hr <- args[2]
29 sep <- as.character(args[3])
30 longitude <- as.numeric(args[4])
31 latitude <- as.numeric(args[5])
32 spe <- as.numeric(args[6])
33 rec <- as.numeric(args[7])
34 crs <- as.numeric(args[8])
35 reso <- as.numeric(args[9])
36 source(args[10])
37 source(args[11])
38 source(args[12])
39 }
40
41 if (hr == "false") {
42 hr <- FALSE
43 }else {
44 hr <- TRUE
45 }
46
47 if (sep == "t") {
48 sep <- "\t"
49 }
50
51 if (crs == "0") {
52 crs <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
53 }
54 #####Import data
55 occ <- read.table(raster, sep = sep, dec = ".", header = hr, fill = TRUE, encoding = "UTF-8") # occ_1M OR occ_SAtlantic
56 occ <- na.omit(occ)
57 #Get biological occurrences
58 #Use the 1 million records subsampled from the full OBIS dataset
59 colnames(occ)[longitude] <- c("decimalLongitude")
60 colnames(occ)[latitude] <- c("decimalLatitude")
61 colnames(occ)[spe] <- c("species")
62 colnames(occ)[rec] <- c("records")
63
64 #Create a discrete global grid
65 #Create an ISEA discrete global grid of resolution 9 using the dggridR package:
66
67 dggs <- dggridR::dgconstruct(projection = "ISEA", topology = "HEXAGON", res = reso)
68
69 #Then assign cell numbers to the occurrence data
70 occ$cell <- dggridR::dgGEO_to_SEQNUM(dggs, occ$decimalLongitude, occ$decimalLatitude)[["seqnum"]]
71
72 #Calculate indicators
73 #The following function calculates the number of records, species richness, Simpson index, Shannon index, Hurlbert index (n = 50), and Hill numbers for each cell.
74
75 #Perform the calculation on species level data
76 idx <- calc_indicators(occ)
77 write.table(idx, file = "Index.csv", sep = ",", dec = ".", na = " ", col.names = TRUE, row.names = FALSE, quote = FALSE)
78
79 #add cell geometries to the indicators table (idx)
80 grid_idx <- sf::st_wrap_dateline(dggridR::dgcellstogrid(dggs, idx$cell))
81 colnames(grid_idx) <- c("cell", "geometry")
82
83 grid <- dplyr::left_join(grid_idx,
84 idx,
85 by = "cell")
86
87 #Plot maps of indicators
88 #Let’s look at the resulting indicators in map form.
89 #Indice ES(50)
90 es_50_map <- gmap_indicator(grid, "es", label = "ES(50)", crs = crs)
91 es_50 <- ggplot2::ggsave("ES_50.png", es_50_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE)
92
93 # Shannon index
94 shannon_map <- gmap_indicator(grid, "shannon", label = "Shannon index", crs = crs)
95 shannon <- ggplot2::ggsave("Shannon_index.png", shannon_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE)
96
97
98 # Number of records, log10 scale, Geographic projection
99 records_map <- gmap_indicator(grid, "n", label = "# of records", trans = "log10", crs = crs)
100 records <- ggplot2::ggsave("Records.png", records_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE)
101
102 # Simpson index
103 simpson_map <- gmap_indicator(grid, "simpson", label = "Simpson index", crs = crs)
104 simpson <- ggplot2::ggsave("Simpson_index.png", simpson_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE)
105
106 # maxp
107 maxp_map <- gmap_indicator(grid, "maxp", label = "maxp index", crs = crs)
108 maxp <- ggplot2::ggsave("Maxp.png", maxp_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE)
109
110 #Mapping
111 es_50
112 shannon
113 simpson
114 maxp
115 records