## ------------------------------------------------------------------------
library(raster)
library(rgdal)
library(sf) # new package for geographic vector datasets (replaces sp)
library(mapview) #mapView
library(RANN) # nn2
library(data.table) # fread
library(fasterize)
library(tidyverse) # contains dplyr
library(knitr)
library(date)
library(igraph)
library(randomcoloR)
library(concaveman)
## ---- echo=FALSE---------------------------------------------------------
#knitr::knit_exit()
## ---- echo=FALSE---------------------------------------------------------
setwd("\\\\pinea\\mlc\\Aulas\\CURSOS_R\\isa-setembro-2018-grandes-cdg\\aulas")
## ------------------------------------------------------------------------
af<-data.table::fread(file.path(getwd(),"datasets","viirs-fires-2017-world.csv"))
afpt<-af[latitude>36 & latitude<43 & longitude < -6 & longitude > -10 ] # Continental Portugal
latlon<-CRS("+proj=longlat +datum=WGS84")
y<-raster(crs=latlon,resolution=0.01,xmn=-10,xmx=-6,ymn=36,ymx=42)
r<-rasterize(x=data.frame(lon=afpt$longitude,lat=afpt$latitude),y=y,field=afpt$frp, fun=max)
mapview(r)
## ------------------------------------------------------------------------
dt<-date::as.date(x=afpt$acq_date, order = "ymd") %>% as.integer()
hist(dt, xlabel="Julian dates (starting January 1, 1960", ylabel="number of active fires", main="VIIRS active fires: Portugal 2017")
## ------------------------------------------------------------------------
sfpt<- as.data.frame(afpt) %>% st_as_sf(coords = c("longitude","latitude") , crs = 4326) %>% st_transform(crs=3763)
## ------------------------------------------------------------------------
d<-2000 # meters
t<-2 # days
K<-20 # neighbors to consider
## ------------------------------------------------------------------------
knn<-nn2(data=st_coordinates(sfpt),k=K,searchtype="radius",radius=d)
head(knn$nn.idx[,1:10],5)
head(knn$nn.dists[,1:5],3)
## ------------------------------------------------------------------------
neigh<-cbind(1:nrow(knn$nn.idx),knn$nn.idx)
head(neigh[,1:10],4)
edges<-data.table(v1=rep(neigh[,1],ncol(neigh)-1),v2=as.vector(neigh[,-1]))
dim(edges)
## ------------------------------------------------------------------------
edges[v2==0] # non edges
edges<-edges[v2!=0] # remove non edges
## ------------------------------------------------------------------------
edges[,d1:=dt[v1]] # acquisition dates
edges[,d2:=dt[v2]] # acquisition dates
edges<-edges[abs(d1-d2)<=t]
## ------------------------------------------------------------------------
xy<-st_coordinates(sfpt) # for clarity, but optional
edges[,x1:=xy[v1,"X"]] # coordinates, optional
edges[,y1:=xy[v1,"Y"]] # coordinates, optional
edges[,x2:=xy[v2,"X"]] # coordinates, optional
edges[,y2:=xy[v2,"Y"]] # coordinates, optional
edges[1:3,]
## ------------------------------------------------------------------------
max(edges[,sqrt((x1-x2)^2+(y1-y2)^2)])
max(edges[,abs(d1-d2)])
sum(duplicated(edges)) # check that there are no duplicated edges
## ------------------------------------------------------------------------
Gall<-igraph::graph_from_data_frame(d=edges,directed=FALSE)
G<-igraph::simplify(Gall) # removes loops and multiple edges
## ------------------------------------------------------------------------
v<-as.integer(vertex_attr(G)$name) # vertex indices which match row indices of edges
range(v)== range(edges[,.(v1,v2)]) # confirm index numering
## ------------------------------------------------------------------------
CL<-igraph::components(G)
head(CL$membership) # returns vertice name and vertice cluster
head(names(CL$membership)) # names (character)
## ------------------------------------------------------------------------
head(CL$csize)
idxcomp<-2500
# vertices of the component
vids<-as.integer(names(CL$membership)[CL$membership==idxcomp])
# subgraph
S<-igraph::induced_subgraph(G,vids=which(CL$membership==idxcomp))
## ------------------------------------------------------------------------
# colors for subgraph: depending on FRP
cores<-colorRampPalette(c(rgb(1,1,0), rgb(1,0,0)))(length(V(S)))[order(sfpt$frp[vids])]
S <- set.vertex.attribute(S, "color", value=cores)
# color is a new attribute
vertex.attributes(S)
# plot using xy coordinates, create layout for ploting
lxy<-st_coordinates(sfpt)[vids,]
## ----eval=FALSE----------------------------------------------------------
## plot(S,layout=lxy,vertex.label=sfpt$acq_date[vids],vertex.size=30,vertex.color=cores)
## plot(S,layout=lxy,vertex.label=sfpt$acq_time[vids],vertex.size=30,vertex.color=cores)
## ------------------------------------------------------------------------
CL<-igraph::components(G)
sfpt$clust<-CL$membership
## ------------------------------------------------------------------------
n <- 20
palette <- distinctColorPalette(n)
cores<-rep(palette,times=max(sfpt$clust)/n)
mapview(sfpt,color=cores[sfpt$clust],cex=4)
## ------------------------------------------------------------------------
# Louvain
CLlouvain<-igraph::cluster_louvain(G) # returns community object; can use positive weights
sfpt2<-sfpt # make a copy of sfpt for mapview
sfpt2$clust<-CLlouvain$membership
mapview(sfpt,color=cores[sfpt$clust],cex=4)+mapview(sfpt2,color=cores[sfpt2$clust],cex=4)
## ------------------------------------------------------------------------
sfmpt<-sfpt2 %>% group_by(clust) %>% summarize(date=modal(acq_date), N=length(acq_date),do_union=TRUE) %>% st_cast(to="MULTIPOINT")
# number points in each feature
dimfeat<-sapply(st_geometry(sfmpt),function(x) length(x))
sfmpt3<-sfmpt[dimfeat>=6,] # two coordinates per point
## ------------------------------------------------------------------------
fconv<-function(feat)
{
xy<-st_coordinates(feat)[,c("X","Y")];
pol<-xy[chull(xy),];
st_polygon(list(rbind(pol,pol[1,])))
}
## ------------------------------------------------------------------------
sfconv<-sfmpt3 # make a copy
# replace geometry by a polygon geometry applying convex hull to the groups of points
st_geometry(sfconv)<-st_sfc(lapply(st_geometry(sfmpt3),FUN=fconv), crs=st_crs(sfmpt)$proj4string)
mapview(sfconv)+mapview(sfpt2,color=cores[sfpt2$clust],cex=4)
## ------------------------------------------------------------------------
fconc<-function(feat)
{
xy<-st_coordinates(feat)[,c("X","Y")];
pol<-concaveman(xy,concavity=3);
st_polygon(list(rbind(pol,pol[1,])))
}
sfconc<-sfmpt3 # make a copy
st_geometry(sfconc)<-st_sfc(lapply(st_geometry(sfmpt3),FUN=fconc), crs=st_crs(sfmpt)$proj4string)
## ----echo=FALSE, eval=FALSE----------------------------------------------
## aa<-st_read(file.path(getwd(),"datasets","Areas_Ardidas_2017_Portugal_Continental_v1.shp"),quiet=TRUE)
## save(aa,sfpt2,sfconc,sfconv,file="patches.RData")
## setwd("\\\\pinea\\mlc\\Aulas\\CURSOS_R\\isa-setembro-2018-grandes-cdg\\aulas")
## load("patches.RData")
## library("htmlwidgets")
## library(mapview)
## library(sf)
## library("randomcoloR")
## n <- 20
## palette <- distinctColorPalette(n)
## cores<-rep(palette,times=max(sfpt2$clust)/n)
## aa<-st_read(file.path(getwd(),"datasets","Areas_Ardidas_2017_Portugal_Continental_v1.shp"),quiet=TRUE)
## m<-mapview(sfpt2,color=cores[sfpt2$clust],cex=4)+mapview(aa,col.regions=distinctColorPalette(20))+mapview(sfconc,col.regions="gray")
## # create html file from map m
## saveWidget(m@map, file="areas_ardidas_pt_2017.html", selfcontained = FALSE)
## ------------------------------------------------------------------------
aa<-st_read(file.path(getwd(),"datasets","Areas_Ardidas_2017_Portugal_Continental_v1.shp"),quiet=TRUE)
n <- 20
palette <- distinctColorPalette(n)
cores<-rep(palette,times=max(sfpt2$clust)/n)
mapview(sfpt2,color=cores[sfpt2$clust],cex=4,layer.name="active fires")+mapview(aa,col.regions=distinctColorPalette(20),layer.name="burned areas")+mapview(sfconc,col.regions="gray", layer.name="s-p patches")
## ------------------------------------------------------------------------
latlon<-CRS("+proj=longlat +datum=WGS84")
r<-raster(crs=latlon,resolution=0.0005,xmn=-10,xmx=-6,ymn=36,ymx=42)
r<-projectRaster(from=r, crs=CRS(st_crs(aa)$proj4string))
sfconc<-st_transform(sfconc,crs=st_crs(aa)$proj4string)
## ------------------------------------------------------------------------
system.time({
sfconc$new<-1
rconc<-fasterize(sf=sfconc,raster=r,field="new",fun="max",background=0)
aa$new<-1
raa<-fasterize(sf=aa,raster=r,field="new",fun="max",background=0)
})
## ------------------------------------------------------------------------
DT<-data.table(aa=values(raa),conc=values(rconc))
DT[,.N,by=.(aa,conc)]