## Testeo y cacharreo
library(raster)
source('./inst/shiny-examples/respineapp/createLandscape.R')
source('./inst/shiny-examples/respineapp/initRichness.R')
source('./inst/shiny-examples/respineapp/dist2nf.R')
source('./inst/shiny-examples/respineapp/disper.R')
source('./inst/shiny-examples/respineapp/disper_time.R')
source('./inst/shiny-examples/respineapp/adjfactor.R')
library(landscapeR)
library(rasterVis)
library(rgeos)
### 1 ####
set.seed(333)
# m <- matrix(nrow=53, ncol=63, byrow = T)
# r <- raster(m)
# extent(r) <- matrix(c(0, 0, 630, 530), nrow=2)
# r[] <- 0
ancho <- 63 * 2
alto <- 53 * 2
m <- matrix(nrow=alto, ncol=ancho, byrow = T)
r <- raster(m)
extent(r) <- matrix(c(0, 0, ancho, alto), nrow=2)
r[] <- 0
## inputs
size_pp <- 1000
den_pp <- 900
size_nf <- 500
n_nf <- 5
pastUse <- 'Oak' # TODO switch past use
# Landscape parameteres
#size_pp <- input$size_pp # size pine plantations
# den_pp <- input$density_pp # density pine plantations
# size_nf <- input$size_nf # size natural forests
# n_nf <- input$n_nf # n patchs natural forests
### 2 ### CREA LANDSCAPE
myl <- createLandscape(r, size_pp = size_pp, size_nf = size_nf, n_nf = n_nf)
makeClass
pp <- makeClass(r, val=1, npatch = 1, rast=TRUE,
size = size_pp,
pts = matrix(c(28*2,25*2), nrow=1, ncol=2))
puntito <- makeClass(r, val=99, npatch = 1, rast = TRUE,
size = 1,
pts = matrix(c(28*2,25*2), nrow=1, ncol=2))
puntito_v <- rasterToPoints(puntito, fun=function(x){x==99})
puntitoBuffer <- gBuffer(puntito_v,
width = c(round(sqrt(200/pi)), round(sqrt(200/pi))))
pB <- dismo:::.generateCircles(puntito_v, d=round(sqrt(200/pi)), lonlat = FALSE, crs=NA)
puntito <- as.data.frame(cbind(x=50, y=56))
coordinates(puntito)<-~x+y
pun <- circles(puntito, d=round(sqrt((size_pp-100)/pi)), lonlat = FALSE)
pp_value <- 1
pine_limit <- rasterToPolygons(myl, fun=function(x){x==pp_value}, dissolve = TRUE)
# get centroid
centroids <- as.data.frame(getSpPPolygonsLabptSlots(pine_limit))
names(centroids) <- c('x', 'y')
coordinates(centroids) <- ~x+y
pun <- circles(centroids, d=round(sqrt((size_pp)/pi)), lonlat = FALSE)
plot(myl)
plot(pun, add=TRUE)
pp <- makeClass(r, val=1, npatch = 1, rast=TRUE,
size = size_pp,
pts = matrix(c(28*2,25*2), nrow=1, ncol=2))
pol_pine <- rasterToPolygons(pp, fun=function(x){x==pp_value}, dissolve = TRUE)
centroids <- as.data.frame(getSpPPolygonsLabptSlots(pol_pine))
names(centroids) <- c('x', 'y')
coordinates(centroids) <- ~x+y
pun <- circles(centroids, d=round(sqrt(((size_pp-50))/pi)), lonlat = FALSE)
pp <- crop(pp, pun@polygons)
plot(pp)
plot(pun, add=TRUE)
### 3 ### PLOT LANDSCAPE CREADO
colour_tree_density <- '#238b45'
colores <- c('lightgoldenrod1', # Crops
'green', # Natural forests
'gray99', # Other
colour_tree_density) # Pine plantation
## Legend
myKey <- list(text = list(lab = c("crop", "natural forest","other", "pine")),
rectangles=list(col = colores),
space='bottom', columns=4)
levelplot(myl, att='landuse', scales=list(draw=FALSE),
col.regions = colores, colorkey=FALSE, key = myKey)
### 4 ### VALORES RIQUEZA INICIALES
# nf_value <- 2
#
pp_value <- 1
pine_limit <- rasterToPolygons(myl, fun=function(x){x==pp_value}, dissolve = TRUE)
levelplot(myl, att='landuse', scales=list(draw=FALSE),
col.regions = colores, colorkey=FALSE, key = myKey) +
layer(sp.polygons(pine_limit, fill=NA))
myr_range <- as.data.frame(
cbind(value = c(0,1,2,3),
lowRich = c(0, 12.82, mean(13.72, 15.62), 1),
upRich = c(0, 13.34, mean(16.11, 19.66), 2)))
### Raster de distancias. Distancias desde los bordes de las manchas de Natural forest hacia el retso de celdas del espacio
dist_raster <- dist2nf(myl, nf_value = 2)
## Visualiza el resultado de las distancias
# Esto para eliminar la tabla de atributos
# dist_raster <- setValues(raster(dist_raster), dist_raster[])
#
# levelplot(dist_raster, scales=list(draw=FALSE), margin=FALSE,
# par.settings = RdBuTheme) +
# layer(sp.polygons(nf_edges, fill=NA))
# Computa riquiza inicial
mapa_riqueza <- initRichness(r = myl,
draster = dist_raster,
r_range = myr_range,
treedensity = den_pp,
pastUse = pastUse,
rescale = FALSE)
# nf_ri <- calc(stack(myl, mapa_riqueza), function())
# rich_pp <- reactive({
# calc(stack(landscapeInit(), rasterRich()), fun=function(x) ifelse(x[1] == pp_value, x[1]*x[2], NA))
# })
nf_value <- 2
rich_nf <- calc(stack(myl, r), fun=function(x) ifelse(x[1] == nf_value, (x[1]/nf_value)*x[2], NA))
plot(rich_nf)
cellStats(rich_nf, min)
cellStats(rich_nf, mean)
cellStats(rich_nf, max)
crop_value <- 3
rich_crop <- calc(stack(myl, r), fun=function(x) ifelse(x[1] == crop_value, (x[1]/crop_value)*x[2], NA))
plot(rich_crop)
cellStats(rich_crop, min)
cellStats(rich_crop, mean)
cellStats(rich_crop, max)
pp_value <- 1
rich_pp <- calc(stack(myl, r), fun=function(x) ifelse(x[1] == pp_value, (x[1]/crop_value)*x[2], NA))
plot(rich_pp)
cellStats(rich_pp, min)
cellStats(rich_pp, mean)
cellStats(rich_pp, max)
rich_nf <- reactive({
rich_nf <- calc(stack(myl, r), fun=function(x) ifelse(x[1] == nf_value, (x[1]/nf_value)*x[2], NA))
list(
min = cellStats(rich_nf, min),
mean = cellStats(rich_nf, mean),
max = cellStats(rich_nf, max)
)
})
# Natural forest (2)
aux <- r_range[which(r_range$value == 2), ]
r[r == 2] <- sample(aux$lowRich:aux$upRich, ncell(r[r==2]), replace = TRUE)
plot(r)
r[r == 0] <- NA
r[r == 1] <- -100
r[r == 2] <- -200
r[r == 3] <- -300
aux <- r_range[which(r_range$value == 2), ]
r[r == -200] <- sample(aux$lowRich:aux$upRich, ncell(r[r==-200]), replace = TRUE)
aux <- r_range[which(r_range$value == 3), ]
r[r == -300] <- sample(aux$lowRich:aux$upRich, ncell(r[r==-300]), replace = TRUE)
ftreeden <- exp(-0.5*((treedensity - 0.22)/1504.1)^2)
sh <- calc(draster, fun=function(x){1.7605 - 0.0932*(sqrt(sqrt(x)))})
# Create a stack with the shanon diversity raster and landuse raster, and then compute values for pine plantations
s <- calc(stack(r, sh), fun=function(x) ifelse(x[1] == -100 , (x[1]/-100)*x[2], NA))
# Scale the distance effect from 0 to 1
sh_scaled <- (s - cellStats(s, "min"))/(cellStats(s, "max") - cellStats(s, "min"))
## ~ PastUSE
### Past Land Use
fplu <- ifelse(pastUse == 'Oak', .9999,
ifelse(pastUse == 'Shrubland', .4982,
ifelse(pastUse == 'Crop', .0279, .0001)))
###### ---------------------------------------
## Combine factor to correct pine plantations (OJO!!!!!!!!!!!!!! PESOS)
f_pine <- (sh_scaled*0.35) + (.45*ftreeden + .2*fplu)
aux <- r_range[which(r_range$value == 1), ]
r[r == -100] <- sample(aux$lowRich:aux$upRich, ncell(r[r==-100]), replace = TRUE)
r <- calc(stack(r, f_pine),
fun = function(x) ifelse(x[1] == 1, x[1]*x[2], x[1]))
plot(r)
plot(sh)
plot(s)
r
# nf_ri <- calc(stack(myl, mapa_riqueza), function())
# rich_pp <- reactive({
# calc(stack(landscapeInit(), rasterRich()), fun=function(x) ifelse(x[1] == pp_value, x[1]*x[2], NA))
# })
nf_value <- 2
rich_nf <- calc(stack(myl, mapa_riqueza), fun=function(x) ifelse(x[1] == nf_value, (x[1]/nf_value)*x[2], NA))
plot(rich_nf)
cellStats(rich_nf, mean)
crop_value <- 3
rich_crop <- calc(stack(myl, mapa_riqueza), fun=function(x) ifelse(x[1] == crop_value, (x[1]/crop_value)*x[2], NA))
plot(rich_crop)
cellStats(rich_crop, min)
###### Disperse module
v <- disper(x=myl, xr=mapa_riqueza, nf_value = 2, pp_value = 1)
# levelplot(dist_raster, scales=list(draw=FALSE), margin=FALSE,
# par.settings = RdBuTheme)
levelplot(stack(v[['msb']],
v[['mmb']],
v[['mma']]),
margin=FALSE, par.settings = RdBuTheme)
# percentage disperser
per_sb <- 0.1
per_mb <- 0.3
per_ma <- 0.6
# propaguleInputs
piB <- (3.7)/10
piM <- (0.2)/10
vv <- disper_time(msb = v[['msb']],
mmb = v[['mmb']],
mma = v[['mma']],
x = myl, xr=mapa_riqueza, pp_value = 1,
per_sb = per_sb, per_mb = per_mb, per_ma = per_ma,
propaguleInputBird = piB,
propaguleInputMammal = piM,
time_span = 10)
# crs(vv) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#writeRaster(vv, '/Users/ajpelu/Desktop/stackTEST/out.tif', bylayer=TRUE)
msb = v[['msb']]
mmb = v[['mmb']]
mma = v[['mma']]
x = myl
xr=mapa_riqueza
pp_value = 1
propaguleInputBird = piB
propaguleInputMammal = piM
time_span = 10
levelplot(stack(vv),
margin=FALSE, par.settings = RdBuTheme)
x=myl
xr=mapa_riqueza
nf_value = 2
pp_value = 1
# Cosas que organizar para el adjacency module
# https://stackoverflow.com/questions/45338384/calculate-the-length-of-shared-boundaries-between-multiple-polygons
ad <- c(0.606, 10.2, 24.2, 33.4, 21.6, 7.53, 39.0, 44.1, 99.5)
seedlim <- c(0.797, 0.800, 0.789, 0.666, 0.666, 0.526, 0.468, 0.395, 0.392)
df_ad <- as.data.frame(cbind(ad, seedlim))
model_ad <- lm(seedlim~ad, data=df_ad)
adj <- c(0,25,50,75,100)
adj <- as.data.frame(adj) %>%
mutate(sl = 0.736658946 -0.004037077 * adj,
sli = 1/sl,
zadj = (sl - min(sl)) / (max(sl) - min(sl)),
zadji = (sli - min(sli)) / (max(sli) - min(sli)),
zadji1 = zadji + 1)
## ~ TreeDensity
### Fraction of Potential Richness (tree Density Eq. 3 Gomez Aparicio et al. 2009)
treedensity <- den_pp
ftreeden <- exp(-0.5*((treedensity - 0.22)/1504.1)^2)
## ~ PastUSE
### Past Land Use
fplu <- ifelse(pastUse == 'Oak', .9999,
ifelse(pastUse == 'Shrubland', .4982,
ifelse(pastUse == 'Crop', .0279, .0001)))
# Get contribution factor of tree Density and Past Land Use
f2 <- (.6*ftreeden + .1*fplu)
sh <- calc(dist_raster, fun=function(x){1.7605 - 0.0932*(sqrt(sqrt(x)))})
# Create a stack with the shanon diversity raster and landuse raster, and then compute values for pine plantations
s <- calc(stack(myl, sh), fun=function(x) ifelse(x[1] == 1 , x[1]*x[2], NA))
sh_scaled <- (s - cellStats(s, "min"))/(cellStats(s, "max") - cellStats(s, "min"))
r_range <- myr_range
ola <- initRichness(r, r_range, treedensity, pastUse, rescale=TRUE)
sh <- st_as_sf(aux_shape)
Touching_List <- st_touches(sh)
line <- st_intersection(sh[1,], sh[Touching_List[[1]][1],])
l_line <- st_length(line)
sh
Shapefile.sf <- st_as_sf(aux_shape)
# ---- Touching list ----cl
Touching_List <- st_touches(Shapefile.sf)
# ---- Polygons perimeters ----
perimeters <- st_length(Shapefile.sf)
from <- 1
lines <- st_intersection(Shapefile.sf[from,], Shapefile.sf[Touching_List[[from]],])
lines <- st_cast(lines) # In case of multiple geometries (ex. from=71)
l_lines <- st_length(lines)
library(ks)
library(KernSmooth)
library(raster)
set.seed(0)
GPS <- data.frame(lon=runif(100), lat=runif(100)*2)
est <- bkde2D(GPS, bandwidth=0.1, gridsize = c(400L, 400L))
names(est) <- c('x', 'y', 'z')
est.raster <- raster(est)
plot(est.raster)
??
https://stackoverflow.com/questions/35752736/creating-a-raster-file-from-a-kernel-smooth-with-differing-x-y-resolutions
https://cran.r-project.org/web/packages/splancs/splancs.pdf
# ver libro virgilio Applied Spatial Data Analysis with R kernel2d
library(dplyr)
library(ggplot2)
library(fitdistrplus)
# --------------------
# Arrendajo from Gomez 2003
dis <- c(10, 30, 50, 70, 90, 150, 250, 350, 450, 550, 700, 900)
per <- c(1.07, 1.40, 2.34, 3.78, 9.82, 33.01, 16.11, 11.40, 10.06, 10.33, 4.89, 0.77)
df <- data.frame(cbind(dis, per))
df <- df %>%
mutate(k = -0.1 + 0.02*log(dis) + 0.04*(log(dis^2)),
u = 1/k)
df %>% ggplot(aes(x=dis, y=u)) + geom_point() + geom_smooth()
# Fit a lognormal and weibull
arre_lnorm <- fitdist(df$per, "lnorm")
arre_w <- fitdist(df$per, "weibull")
# --------------------
grid <- seq(0,1500,1)
# Arrendajo (Pons & Pausas 2007)
plot(dlnorm(grid,meanlog = 3.844, sdlog = 0.851), col='blue',
type="l")
# kernel para mamĂferos González-Varo et al. 2012 http://onlinelibrary.wiley.com/doi/10.1111/1365-2656.12024/full
lines(grid,
ifelse(grid <= 400,
dweibull(grid, shape = 1.385, scale = 137),
dlnorm(grid, meanlog = 6.621, sdlog = 0.297)),
type="l",xlab="x",ylab="f(x)", col='red')
lines(grid,
dlnorm(grid, meanlog = log(200),
sdlog= arre_lnorm$estimate["sdlog"]), col='blue',
lty=2)
lines(grid,
dweibull(grid,
shape = arre_w$estimate["shape"],
scale = 88),
col='blue', lty=6)
lines(grid,
dlnorm(grid,meanlog = arre_lnorm$estimate["meanlog"],
sdlog= arre_lnorm$estimate["sdlog"]), col='blue',
lty=2)
dis <- data.frame(d = seq(1,1000, b=1))
dis <- dis %>%
mutate(k = -0.1 + 0.02*d + 0.04*(d^2),
u = 1/k,
u2 = 1/(-0.1 - 0.02*(log(d)) + 0.04*((log(d))^2)),
gu = 1/(0.14 - 2.17*(log(d)) - 2.50*((log(d))^2) + 0.7*((log(d))^3)))
dis %>% filter(u > 0) %>%
ggplot(aes(x=d, y=gu)) + geom_line()
#### Proposal
data <- data.frame(
dis = seq(0,1500,1))
data <- data %>%
mutate(
sb = dlnorm(dis, meanlog = log(51), sdlog = .9),
mb = dlnorm(dis, meanlog = log(151), sdlog = 1),
ma = ifelse(dis <= 400,
dweibull(dis, shape = 1.385, scale = 137),
dlnorm(dis, meanlog = 6.621, sdlog = 0.297)),
gg = dlnorm(dis, meanlog = log(450), sdlog = 1))
data %>% ggplot(aes(x=dis)) +
geom_line(aes(x=dis, y=sb), color='blue') +
geom_line(aes(x=dis, y=mb), color='black') +
geom_line(aes(x=dis, y=ma), color='red') +
geom_line(aes(x=dis, y=gg), color='green') +
scale_x_continuous(breaks = seq(0,1500, 100)) +
theme_bw() +
theme(panel.grid = element_blank()) +
ylab('f(x)') +
annotate("text",
x= 1000, y = 0.010,
label=("blue = sb, red = ma, black=mb, green = gg"))
plot(columbus, border="grey")
plot(col.gal.nb, coords, add=TRUE)
plot(dxx, coords, add=TRUE, col="red")
title(main=paste("Differences (red) in Columbus GAL weights (black)",
"and polygon generated queen weights", sep="\n"))
xxx <- poly2nb(columbus, queen=FALSE)
dxxx <- diffnb(xxx, col.gal.nb)
plot(columbus, border = "grey")
plot(col.gal.nb, coords, add = TRUE)
plot(dxxx, coords, add = TRUE, col = "red")
title(main=paste("Differences (red) in Columbus GAL weights (black)",
"and polygon generated rook weights", sep="\n"))
cards <- card(xx)
maxconts <- which(cards == max(cards))
if(length(maxconts) > 1) maxconts <- maxconts[1]
fg <- rep("grey", length(cards))
fg[maxconts] <- "red"
fg[xx[[maxconts]]] <- "green"
plot(columbus, col=fg)
title(main="Region with largest number of contiguities")
example(nc.sids)
system.time(xxnb <- poly2nb(nc.sids))
plot(nc.sids)
plot(xxnb, coordinates(nc.sids), add=TRUE, col="blue")
require("rgdal")
require("rgeos")
library(sp)
library(raster)
download.file("https://www.dropbox.com/s/vbxx9dic34qwz63/Polygons.zip?dl=1", "Polygons.zip")
unzip("./Polygons.zip")
Shapefile <- readOGR(".","Polygons")
plot(edges, col=edges@data)
Touching_List <- gTouches(edges, byid = TRUE, returnDense=FALSE)
# Touching_DF <- setNames(utils::stack(lapply(Touching_List, as.character)), c("TOUCHING", "ORIGIN"))
# ---- Calculate perimeters of all polygons ----
perimeters <- sp::SpatialLinesLengths(as(edges, "SpatialLines"))
# ---- Example with the first object of the list and first neighbor ----
from <- 1
to <- 1
line <- rgeos::gIntersection(edges[from,], edges[Touching_List[[from]][to],])
l_line <- sp::SpatialLinesLengths(line)
plot(edges[c(from, Touching_List[[from]][to]),])
plot(line, add = TRUE, col = "red", lwd = 2)
# ---- Example with the first object of the list and all neighbors ----
from <- 1
lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE)
l_lines <- sp::SpatialLinesLengths(lines)
plot(Shapefile[c(from, Touching_List[[from]]),])
plot(lines, add = TRUE, col = 1 + 1:length(Touching_List[[from]]), lwd = 2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.