Mercurial > repos > ecology > obis_data
comparison obisindicators.r @ 0:1fcd81d65467 draft
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/obisindicators commit b377ff767e3051f301c2f02cfe3e1a17b285ede4
author | ecology |
---|---|
date | Thu, 18 Jan 2024 09:33:52 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:1fcd81d65467 |
---|---|
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 |