.Rhistory

# Plot the transformed probabilities
plot(multipl.years$date, multipl.years$smoothed_prob_presence, type = "l", xlab = "Date", ylab = "Probability of Presence", main = "Transformed Probability of Presence")
pa.ts <- ts(multipl.years$smoothed_prob_presence, start = 1, frequency = 52) ### for computing the sine/cosine waves is the same as zoo above - start = 1 is ok even if the first starting time/day is at
pa.ts
aaarrr2 <- data.frame(forecast::fourier(pa.ts, K=2)) ### on average enough
aaarrr2
aaarrr$abb<- multipl.years$abundance
aaarrr$pa<-multipl.years$pa
aaarrr$pa<-multipl.years$smoothed_prob_presence
aaarrr2 <- data.frame(forecast::fourier(pa.ts, K=2)) ### on average enough
aaarrr2$abb<- multipl.years$abundance
aaarrr2$pa<-multipl.years$smoothed_prob_presence
#---- 4. Derive probability of presence of the species over time ----
m1 <- glm(abb~., family="poisson", data=aaarrr)
summary(m1)
plot(m1)
plot(multipl.years$date, m1$fitted.values)
m2 <- glm(pa~., family="binomial", data=aaarrr)
summary(m2)
plot(m2)
plot(multipl.years$date, m2$fitted.values)
m1.1 <- glm(abb~., family="poisson", data=aaarrr2)
summary(m1.1)
plot(m1.1)
plot(multipl.years$date, m1.1$fitted.values)
m2.2 <- glm(pa~., family="binomial", data=aaarrr2)
summary(m2.2)
plot(multipl.years$date, m2.2$fitted.values)
# multiple location and compute photoperiod as proxy of seasonality
set.seed(123)
# Define the number of points to simulate
num_points <- 12
# Simulate longitude and latitude coordinates for Europe
# For demonstration purposes, we'll generate coordinates within a bounding box roughly covering Europe
# You may want to adjust this bounding box based on your specific area of interest
min_lon <- -10  # Minimum longitude (Western Europe)
max_lon <- 30   # Maximum longitude (Eastern Europe)
min_lat <- 35   # Minimum latitude (Southern Europe)
max_lat <- 60   # Maximum latitude (Northern Europe)
# Generate random longitude and latitude coordinates within the specified bounding box
lon <- runif(num_points, min = min_lon, max = max_lon)
lat <- runif(num_points, min = min_lat, max = max_lat)
# Create a data frame to store the coordinates
points_data <- data.frame(lon, lat)
# Print the simulated spatial points
print(points_data)
st_as_sf(points_data, coord=c("lon", "lat"), crs=4326)
library(sf)
st_as_sf(points_data, coord=c("lon", "lat"), crs=4326)
points_data
?st_as_sf
st_as_sf(x=points_data, coord=c("lon", "lat"), crs=4326)
st_as_sf(x=points_data, coords=c("lon", "lat"), crs=4326)
sp.pts <- st_as_sf(x=points_data, coords=c("lon", "lat"), crs=4326)
mapview(sp.pts)
library(mapview)
mapview(sp.pts)
# multiple location and compute photoperiod as proxy of seasonality
library(geosphere)
?daylength
# Set seed for reproducibility
set.seed(123)
# Define the start and end date of the time series
start_date <- as.Date("2024-01-01")
end_date <- as.Date("2024-12-31")
# Create a sequence of dates from start_date to end_date with a weekly frequency
dates <- seq(start_date, end_date, by = "week")
length(dates)
dates
daylength(points_data$lat[1], doy=dates)
apply(points_data, 1, function(x){daylength(points_data$lat, doy=dates)} )
sapply(points_data, function(x){daylength(points_data$lat, doy=dates)} )
sapply(points_data$lat, function(x){daylength(x, doy=dates)} )
# Print the simulated spatial points
points_data$id<-1:n()
# Print the simulated spatial points
points_data$id<-1:nrow(points_data)
print(points_data)
photo<-sapply(points_data$lat, function(x){daylength(x, doy=dates)} )
names(photo)<-points_data$id
photo
photo<-data.frame(sapply(points_data$lat, function(x){daylength(x, doy=dates)} ))
names(photo)<-points_data$id
photo
t(photo)
photo<-t(data.frame(sapply(points_data$lat, function(x){daylength(x, doy=dates)} )))
names(photo)<-dates
photo
photo<-t(data.frame(sapply(points_data$lat, function(x){daylength(x, doy=dates)} )))
photo
photo<-data.frame(t(sapply(points_data$lat, function(x){daylength(x, doy=dates)} )))
names(photo)<-dates
photo
library(dplyr)
library(tidyr)
photo<-cbind(points_data, photo)
photo
photo<-data.frame(t(sapply(points_data$lat, function(x){daylength(x, doy=dates)} )))
names(photo)<-dates
photo<-bind_cols(points_data, photo) %>%
pivot_longer(!c(ID, lon, lat))
photo<-bind_cols(points_data, photo) %>%
pivot_longer(!c(id, lon, lat))
photo
library(ggplot2)
photo %>%
mutate(name=as.Date(name))
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name)
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value)
ggplot(aes(x=date, y=photo, col=id))+
geom_line()+
theme_classic()
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value) %>%
ggplot(aes(x=date, y=photo, col=id))+
geom_line()+
theme_classic()
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value) %>%
ggplot(aes(x=date, y=photo, col=id, group=id))+
geom_line()+
theme_classic()
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value) %>%
ggplot(aes(x=date, y=photo, col=id, group=id))+
geom_line()+
labs(y="Daylength", x="Date")
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value) %>%
ggplot(aes(x=date, y=photo, col=id, group=id))+
geom_line()+
labs(y="Daylength", x="Date")+
theme_classic()
photo %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value) %>%
ggplot(aes(x=date, y=photo, col=lat, group=id))+
geom_line()+
labs(y="Daylength", x="Date")+
theme_classic()
sp.pts
sp.pts <- st_as_sf(x=points_data, coords=c("lon", "lat"), crs=4326)
sp.pts
mapview(sp.pts, zcol="id")
# Set seed for reproducibility
set.seed(123)
# Define the start and end date of the time series
start_date <- as.Date("2024-01-01")
end_date <- as.Date("2024-12-31")
# Create a sequence of dates from start_date to end_date with a weekly frequency
dates <- seq(start_date, end_date, by = "week")
length(dates)
# Define the seasonal pattern with peak at the end of July and abundance 0 around end of October (example with the Tiger mosquito in mind)
seasonal_pattern <- c(rep(0, 12),        # January to end of March: 0
seq(5, 20, by = 2),        # April to May: Low abundance (5)
seq(25, 50, by = 2),  # June-July to August: Increasing abundance to the peak (10 to 50)
seq(45, 30, length.out = 7), # August to September: decreasing abundance (45 to 30)
seq(25, 0, length.out = 5),  # September to October: Decreasing abundance (25 to 0)
rep(0, 8))         # November to December: 0
length(seasonal_pattern)
# Simulate counts from a Poisson process
counts <- rpois(length(dates), lambda = seasonal_pattern)
# Add jitter to the counts; repeat for three years to get enough variability
jitter_factor <- c(0.2, 0.1, 0.15)  # Adjust the amount of jitter as needed
multipl.years <- do.call(c, lapply(jitter_factor, function(x){jit <- jitter(ifelse(counts == 0, NA, counts), amount = max(counts) * x);
jit <- ifelse(is.na(jit), 0, jit); counts_jittered <- counts + jit; return(counts_jittered)})
)
multipl.years<- data.frame("date"= seq.Date(start_date, length.out = length(multipl.years) , by = "week"),
"abundance" = abs(multipl.years))
multipl.years$week <-lubridate::week(multipl.years$date)
# jit <- jitter(ifelse(counts == 0, NA, counts), amount = max(counts) * jitter_factor)
# jit <- ifelse(is.na(jit), 0, jit)
# counts_jittered <- counts + jit
# Plot the simulated counts with jitter
plot(multipl.years$date, multipl.years$abundance, type = "l", xlab = "Date", ylab = "Count", main = "Simulated Poisson Process with Seasonality and Jitter")
multipl.years
# assign latitute and longitute of one point to infer the seasonlity of the others
myObs<-multipl.years
myObs
photo
subset(photo, id==7)
photo<-bind_cols(points_data, photo) %>%
pivot_longer(!c(id, lon, lat)) %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value)
photo<-data.frame(t(sapply(points_data$lat, function(x){daylength(x, doy=dates)} )))
names(photo)<-dates
photo<-bind_cols(points_data, photo) %>%
pivot_longer(!c(id, lon, lat)) %>%
mutate(name=as.Date(name)) %>%
rename(date=name, photo=value)
photo
photo%>%
ggplot(aes(x=date, y=photo, col=lat, group=id))+
geom_line()+
labs(y="Daylength", x="Date")+
theme_classic()
subset(photo, id==7, select = c("id", "lon", "lat", "photo"))
# assign latitute and longitute of one point to infer the seasonlity of the others
myObs<-cbind(multipl.years,
subset(photo, id==7, select = c("id", "lon", "lat", "photo"))
)
myObs
myObs
PredMod<-glm(abundance~photo, data=myObs, family="poisson")
summary(PredMod)
plot(PredMod$fitted.values, PredMod$data$abundance)
photo
predict(PredMod, photo, type = "response")
photo$pred_ab<-predict(PredMod, photo, type = "response")
photo
photo%>%
ggplot(aes(x=date, y=pred_ab, col=lat, group=id))+
geom_line()+
labs(y="Daylength", x="Date")+
theme_classic()
install.packages("automap")
library(automap)
loadMeuse()
library(automap)
data(meuse)
data("meuse")
library(sp)
library(sf)
library(stars)
data(meuse)
coordinates(meuse) =~ x+y
data(meuse.grid)
gridded(meuse.grid) =~ x+y
# Ordinary kriging, no new_data object
kriging_result = autoKrige(zinc~1, meuse)
plot(kriging_result)
# Ordinary kriging
kriging_result = autoKrige(zinc~1, meuse, meuse.grid)
plot(kriging_result)
# Fixing the nugget to 0.2
kriging_result = autoKrige(zinc~1, meuse,
meuse.grid, fix.values = c(0.2,NA,NA))
plot(kriging_result)
# Universal kriging
kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid)
plot(kriging_result)
# Block kriging
kriging_result_block = autoKrige(zinc~soil+ffreq+dist,
meuse, meuse.grid, block = c(400,400))
plot(kriging_result_block)
plot(kriging_result_block)
# Dealing with duplicate observations
data(meuse)
meuse.dup = rbind(meuse, meuse[1,]) # Create duplicate
coordinates(meuse.dup) = ~x+y
kr = autoKrige(zinc~dist, meuse.dup, meuse.grid)
# Extracting parts from the autoKrige object
prediction_spdf = kr$krige_output
sample_variogram = kr$exp_var
variogram_model = kr$var_model
coordinates(meuse) = ~x + y
meuse = st_as_sf(meuse)
meuse.grid = st_as_stars(meuse.grid)
kriging_result = autoKrige(zinc~1, meuse,
meuse.grid, fix.values = c(0.2,NA,NA))
plot(kriging_result)
library(sp)
library(spacetime)
sumMetricVgm <- vgmST("sumMetric",
space = vgm( 4.4, "Lin", 196.6,  3),
time  = vgm( 2.2, "Lin",   1.1,  2),
joint = vgm(34.6, "Exp", 136.6, 12),
stAni = 51.7)
library(sp)
library(spacetime)
?vgmST
library(gstat)
sumMetricVgm <- vgmST("sumMetric",
space = vgm( 4.4, "Lin", 196.6,  3),
time  = vgm( 2.2, "Lin",   1.1,  2),
joint = vgm(34.6, "Exp", 136.6, 12),
stAni = 51.7)
sumMetricVgm
data(air)
air
suppressWarnings(proj4string(stations) <- CRS(proj4string(stations)))
stations
data(air)
stations
suppressWarnings(proj4string(stations) <- CRS(proj4string(stations)))
plot(stations)
rural = STFDF(stations, dates, data.frame(PM10 = as.vector(air)))
rural
rr <- rural[,"2005-06-01/2005-06-03"]
rr <- as(rr,"STSDF")
rr
x1 <- seq(from=6,to=15,by=1)
x2 <- seq(from=48,to=55,by=1)
DE_gridded <- SpatialPoints(cbind(rep(x1,length(x2)), rep(x2,each=length(x1))),
proj4string=CRS(proj4string(rr@sp)))
DE_gridded
gridded(DE_gridded) <- TRUE
DE_pred <- STF(sp=as(DE_gridded,"SpatialPoints"), time=rr@time)
DE_kriged <- krigeST(PM10~1, data=rr, newdata=DE_pred,
modelList=sumMetricVgm)
gridded(DE_kriged@sp) <- TRUE
stplot(DE_kriged)
stplot(rural)
dim(rural)
dim(rural)
rural
rural@data
head(rural@data)
head(rural@data, n=25)
head(rural@data, n=50)
rr <- rural[,"2005-06-01/2005-06-03"]
rr <- as(rr,"STSDF")
rr
stplot(rr)
DE_gridded <- SpatialPoints(cbind(rep(x1,length(x2)), rep(x2,each=length(x1))),
proj4string=CRS(proj4string(rr@sp)))
DE_gridded
gridded(DE_gridded) <- TRUE
plot(DE_gridded)
var <- variogramST(PM10~1,data=rr,tunit="days",assumeRegular=F,na.omit=T)
var <- variogramST(PM10~1,data=rural,tunit="days",assumeRegular=F,na.omit=T)
library(USE)
library(USE)
library(terra)
# library(raster)
library(virtualspecies)
library(sf)
library(ggplot2)
envData <- USE::Worldclim_tmp
envData <- terra::rast(envData,  type="xyz")
set.seed(123)
random.sp <- virtualspecies::generateRandomSp(raster::stack(envData),
convert.to.PA = FALSE,
species.type = "additive",
realistic.sp = TRUE,
plot = FALSE)
#reclassify suitability raster using a probability conversion rule
new.pres <- virtualspecies::convertToPA(x=random.sp,
beta=0.55,
alpha = -0.05, plot = FALSE)
#Sample true occurrences
presence.points <- virtualspecies::sampleOccurrences(new.pres,
n = 300, # The number of points to sample
type = "presence-absence",
sample.prevalence = 0.6,
detection.probability = 1,
correct.by.suitability = TRUE,
plot = TRUE)
myPres <- presence.points$sample.points[which(presence.points$sample.points$Observed==1), c( "x", "y",  "Observed")]
myPres <- st_as_sf(myPres, coords=c("x", "y"), crs=4326)
#PCA
rpc <- rastPCA(envData, stand = TRUE)
dt <- na.omit(as.data.frame(rpc$PCs[[c("PC1", "PC2")]], xy = TRUE))
dt <- sf::st_as_sf(dt, coords = c("PC1", "PC2"))
myRes <- USE::optimRes(sdf=dt,
grid.res=10,
perc.thr = 20,
showOpt = TRUE,
cr=5)
myRes
myRes <- USE::optimRes(sdf=dt,
grid.res=c(1, 10),
perc.thr = 20,
showOpt = TRUE,
cr=5)
myRes
uniformSampling
# dwelve into uniformSampling
sdf=dt
grid.res=myRes$Opt_res
n.tr = 5
sub.ts = TRUE
n.ts = 2
grid.res= 10
grid <- sf::st_make_grid(sdf, n = grid.res)
grid
plot(grid)
sdf
sdf$ID <- row.names(sdf)
plot(sdf, add=T)
plot(grid, add=T)
Pts_to_sample <- sdf[grid, ]
Pts_to_sample
if(nrow(Pts_to_sample) != 0) Pts_to_sample$ID <- i
seq_along(grid)
grid <- sf::st_make_grid(sdf, n = grid.res)
sdf$ID <- row.names(sdf)
Pseudo_unif <- do.call(rbind, lapply(seq_along(grid), function(i) {
Grid_to_sample <- grid[i]
sdf$ID <- row.names(sdf)
Pts_to_sample <- sdf[Grid_to_sample, ]
if(nrow(Pts_to_sample) != 0) Pts_to_sample$Cell_id <- i
return(Pts_to_sample)
}))
sum(duplicated(Pseudo_unif[[c('ID')]])) #101
plot(st_geometry(sdf), col = 'lightgrey')
plot(grid, add = T)
plot(st_geometry(Pseudo_unif[which(duplicated(Pseudo_unif[[c('ID')]])), ]), col = 'red', add = T)
uniformSampling
Pseudo_unif
plot(st_geometry(Pseudo_unif[which(duplicated(Pseudo_unif[[c('ID')]])), ]), col = 'red', add = T)
which(duplicated(Pseudo_unif[[c('ID')]])), ]
which(duplicated(Pseudo_unif[[c('ID')]]))
#---- debugging ----
debugonce(uniformSampling)
uniformSampling(sdf=dt, grid.res= 10, n.tr = 5,  sub.ts = TRUE,  n.ts = 2)
isTRUE(sub.ts)
abs_val <- sdf[!(sdf$ID %in% res$ID), ]
res_val <- do.call(rbind, lapply(rev(seq_len(length(grid))),
function(i) {
if (isTRUE(verbose))
message(paste("Processing tile", i, sep = " "))
if (isTRUE(plot_proc))
plot(grid[i], col = "red", add = TRUE)
subs <- abs_val[grid[i], ]
if (nrow(subs) <= n.ts) {
return(subs)
}
else {
subs <- subs[sample(nrow(subs), n.ts, replace = FALSE),
]
return(subs)
}
}))
res
duplicated(res$ID)
which(duplicated(res$ID))
res
res[2,"ID"]<-9634
res
which(duplicated(res$ID))
res[which(duplicated(res$ID)), ]
dim(res)
dim(res[distint(res$ID), ])
dim(res[distinct(res$ID), ])
dim(res[duplyr::distinct(res$ID), ])
dim(res[dplyr::distinct(res$ID), ])
dim(res[dplyr::distinct(res$ID), ])
str(res)
dim(res[dplyr::distinct(as.numeric(res$ID)), ])
dim(res[dplyr::distinct(as.numeric(res$ID)), ])
as.numeric(res$ID)
dplyr::distinct(as.numeric(res$ID))
dplyr::distinct(as.numeric(res$ID))
dplyr::distinct(res$ID)
?distinct
?distinct.sf
dplyr::distinct(res$ID)
res[!duplicated(res$ID)]
res[!duplicated(res$ID),]
dim(res)
dim(res[!duplicated(res$ID),])
res
if(duplicated(res$ID)){
print("bea!")
}
duplicated(res$ID)
TRUE %in% duplicated(res$ID)
if(TRUE %in% duplicated(res$ID)){
print("bea!")
}
?warning
if(TRUE %in% duplicated(res$ID)){
# res[!duplicated(res$ID),]
warning("Warning: pseudo-replicate sampled and removed it")
}
c(FALSE, FALSE, TRUE, FALSE)
any(c(FALSE, FALSE, TRUE, FALSE))
#install package
remove.packages("USE")
#------USE-----------------
setwd("/home/daniele/GitHub/USE/")
# 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
# token
library(devtools)
library(roxygen2)
# library(usethat)
# devtools::create("yourPkg")
# devtools::install_github("mattmar/dynamAedes")
devtools::load_all(".") # Working directory should be in the package directory
#update documentation
devtools::document()
#check everything is ok for the CRAN
devtools::check()#vignettes=FALSE
devtools::install(".")
vignette( package="USE")
# if you change the readme description .Rmd, then run
# use_readme_rmd(open = rlang::is_interactive())
# devtools::build_readme()
# pkgdown::clean_site()
build_news(pkg = ".", override = list(), preview = NA)
#website
# devtools::install_github('r-lib/pkgdown')
library(pkgdown)
# if you change the readme description .Rmd, then run
# use_readme_rmd(open = rlang::is_interactive())
# devtools::build_readme()
# pkgdown::clean_site()
build_news(pkg = ".", override = list(), preview = NA)
build_news(pkg = ".", override = list(), preview = NA)
# Run to build the website
pkgdown::build_site("/home/daniele/GitHub/USE/", install = FALSE,
examples = FALSE)