-
Notifications
You must be signed in to change notification settings - Fork 2
/
1_Record_dimension.R
90 lines (68 loc) · 4.3 KB
/
1_Record_dimension.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
library(sf)
library(terra)
library(spatstat)
library(geodata)
library(dplyr)
library(vroom)
# ---------------------------
# Data loading
# ---------------------------
# Definición del directorio de trabajo
wd <- '/GSI'
setwd(wd) # Establecer el directorio de trabajo
dir.create("1_Records/") # Crear un directorio para guardar resultados
# Cargar shapefile del área de estudio
##col <- geodata::gadm(country= "COL", level=0, path= ".") # función desactivada. Sirve para cargar shapefile de Colombia desde servidor remoto del paquete geodata
col <- read_sf('/Carpeta_shapefile', 'shapefile')
# Definición del sistema de referencia geográfico (GCS) y proyectado (CRS)
GRS.geo<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # Geographic Reference to AOI
CRS.proj<-'+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs' # Planar Reference to AOI
# Transformación del shapefile al sistema de coordenadas proyectadas
proj.col <- sf::st_transform(col, crs = CRS.proj)
# Convertir el shapefile a una ventana espacial compatible con spatstat
shape_zoneOwin <- as.owin(proj.col)
# Cargar registros de especies desde un archivo txt
Data <- vroom("Archivo_registros.txt", col_names = TRUE)
Data2 <- Data[, c('gbifID', 'decimalLatitude', 'decimalLongitude_x')] #Select ad first column the name of the id record. And for second and third column, latitude and longitude, respectively.
Data2 <- Data2[!is.na(Data2$decimalLatitude),] # Eliminar registros con latitud NA
Data2 <- Data2[!is.na(Data2$decimalLongitude),] # Eliminar registros con longitud NA
# Eliminar duplicados basados en coordenadas
Data2 <- Data2[!duplicated(Data2[c("decimalLongitude", "decimalLatitude")]), ]
# De acuerdo con García Márquez et al. (2012), esta dimensión se basa en la lista de localidades de colecta como puntos
# con los cuales se genera la "densidad de localidades" ("density of collection localities").
# Convertir a objeto sf y transformar las coordenadas al sistema proyectado
coords <- st_as_sf(Data2, coords = c("decimalLongitude", "decimalLatitude"), crs = GRS.geo)
coordinates.col <- st_transform(coords, crs = CRS.proj)
# Intersección espacial para verificar si las coordenadas están dentro del AOI (Area of Interest)
system.time(over.coords <- st_intersects(coordinates.col, proj.col, sparse = FALSE)); head(over.coords)
Data2 <- Data2[over.coords == TRUE,]
#system.time(over.coords <- st_intersection(coordinates.col, proj.col, sparse = FALSE)); head(over.coords) #Identify coordinates within the AOI
#coordinates.col <- over.coords #Este paso asume que el st_intersect solo genera en su output coordenadas que efextivamente haver overlap (no genera NAs).
# El comando con st_intersection genera el objetvo espacial de los overlap, pero toma más tiempo.
#system.time(over.coords <- st_within(coordinates.col, proj.col)) #otra alternativa, mirar cuál es más rápida
# Crear un patrón de puntos con la ventana espacial definida
coords_df <- st_coordinates(coordinates.col)
p <- ppp(coords_df[, "X"], coords_df[, "Y"], window = shape_zoneOwin, unitname = c("metre", "metres"))
summary(p)
plot(p)
# Estimación de la densidad del kernel
diggle <- bw.diggle(p) # Selección del valor de suavizamiento
plot(diggle)
plot(diggle, xlim= c(0,100), main="Smoothing bandwidth for the kernel estimation")
# Calcular la intensidad suavizada del patrón de puntos (1km de resolución)
system.time(diggle_den <- density.ppp(p, diggle, eps=10000)) # Kernel Smoothed Intensity of Point Pattern 1km - 17 min
# Graficar la densidad suavizada
plot(diggle_den,diggle, main='Gap density diggle 1km')
# Convertir la densidad a un objeto raster
densi_dig1km <- rast(diggle_den)
rescal_dig1km <- densi_dig1km/max(diggle_den) # Reescalar usando el valor máximo
# Asignar el sistema de referencia proyectado al raster
crs(rescal_dig1km) <- CRS.proj
# Proyectar el raster reescalado al sistema de referencia geográfico (WGS84)
rescal_dig1km_wgs84 <- terra::project(rescal_dig1km, crs(GRS.geo))
# Graficar el raster proyectado
plot(rescal_dig1km_wgs84)
# Guardar el raster resultante
writeRaster(rescal_dig1km_wgs84, filename="1_Records/Vac_dens_rescal_1km.tif", overwrite=TRUE)
# Guardar el entorno de trabajo en un archivo RData
save.image(paste0('1_Records/Records_R_object.RData'))