Mercurial > repos > ecology > obisindicators
comparison obisindicators.r @ 0:1015a0070bac draft
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/obisindicators commit 1e0b3b101d2380338030e607e6949331439c6dfc
author | ecology |
---|---|
date | Thu, 29 Dec 2022 13:05:04 +0000 |
parents | |
children | d7b6ff32d072 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:1015a0070bac |
---|---|
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 source(args[9]) | |
36 source(args[10]) | |
37 source(args[11]) | |
38 } | |
39 | |
40 if (hr == "false") { | |
41 hr <- FALSE | |
42 }else { | |
43 hr <- TRUE | |
44 } | |
45 | |
46 if (sep == "t") { | |
47 sep <- "\t" | |
48 } | |
49 | |
50 if (crs == "0") { | |
51 crs <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" | |
52 } | |
53 #####Import data | |
54 occ <- read.table(raster, sep = sep, dec = ".", header = hr, fill = TRUE, encoding = "UTF-8") # occ_1M OR occ_SAtlantic | |
55 occ <- na.omit(occ) | |
56 #Get biological occurrences | |
57 #Use the 1 million records subsampled from the full OBIS dataset | |
58 colnames(occ)[longitude] <- c("decimalLongitude") | |
59 colnames(occ)[latitude] <- c("decimalLatitude") | |
60 colnames(occ)[spe] <- c("species") | |
61 colnames(occ)[rec] <- c("records") | |
62 | |
63 #Create a discrete global grid | |
64 #Create an ISEA discrete global grid of resolution 9 using the dggridR package: | |
65 | |
66 dggs <- dggridR::dgconstruct(projection = "ISEA", topology = "HEXAGON", res = 9) | |
67 | |
68 #Then assign cell numbers to the occurrence data | |
69 occ$cell <- dggridR::dgGEO_to_SEQNUM(dggs, occ$decimalLongitude, occ$decimalLatitude)[["seqnum"]] | |
70 | |
71 #Calculate indicators | |
72 #The following function calculates the number of records, species richness, Simpson index, Shannon index, Hurlbert index (n = 50), and Hill numbers for each cell. | |
73 | |
74 #Perform the calculation on species level data | |
75 idx <- calc_indicators(occ) | |
76 write.table(idx, file = "Index.csv", sep = ",", dec = ".", na = " ", col.names = TRUE, row.names = FALSE, quote = FALSE) | |
77 | |
78 #add cell geometries to the indicators table (idx) | |
79 grid_idx <- sf::st_wrap_dateline(dggridR::dgcellstogrid(dggs, idx$cell)) | |
80 colnames(grid_idx) <- c("cell", "geometry") | |
81 | |
82 grid <- dplyr::left_join(grid_idx, | |
83 idx, | |
84 by = "cell") | |
85 | |
86 #Plot maps of indicators | |
87 #Let’s look at the resulting indicators in map form. | |
88 #Indice ES(50) | |
89 es_50_map <- gmap_indicator(grid, "es", label = "ES(50)", crs = crs) | |
90 es_50 <- ggplot2::ggsave("ES_50.png", es_50_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE) | |
91 | |
92 # Shannon index | |
93 shannon_map <- gmap_indicator(grid, "shannon", label = "Shannon index", crs = crs) | |
94 shannon <- ggplot2::ggsave("Shannon_index.png", shannon_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE) | |
95 | |
96 | |
97 # Number of records, log10 scale, Geographic projection | |
98 records_map <- gmap_indicator(grid, "n", label = "# of records", trans = "log10", crs = crs) | |
99 records <- ggplot2::ggsave("Records.png", records_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE) | |
100 | |
101 # Simpson index | |
102 simpson_map <- gmap_indicator(grid, "simpson", label = "Simpson index", crs = crs) | |
103 simpson <- ggplot2::ggsave("Simpson_index.png", simpson_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE) | |
104 | |
105 # maxp | |
106 maxp_map <- gmap_indicator(grid, "maxp", label = "maxp index", crs = crs) | |
107 maxp <- ggplot2::ggsave("Maxp.png", maxp_map, scale = 0.38, width = 12, height = 7, units = "in", dpi = 300, limitsize = TRUE) | |
108 | |
109 #Mapping | |
110 es_50 | |
111 shannon | |
112 simpson | |
113 maxp | |
114 records |