## ------------------------------------------------------------------------
library(sf) # new package for geographic vector datasets (replaces sp)
library(mapview) #mapView
library(RANN) # nn2
library(tidyverse) # contains dplyr, magrittr
## ---- echo=FALSE---------------------------------------------------------
#knitr::knit_exit()
## ---- echo=FALSE---------------------------------------------------------
setwd("\\\\pinea\\mlc\\Aulas\\CURSOS_R\\isa-setembro-2018-grandes-cdg\\aulas")
library(knitr)
## ------------------------------------------------------------------------
flam<-sf::st_read(dsn=file.path(getwd(),"datasets","Flam2015v4.shp"))
urb<-sf::st_read(dsn=file.path(getwd(),"datasets","urb12015v4.shp"))
## ---- echo=TRUE----------------------------------------------------------
table(st_geometry_type(flam))
table(st_geometry_type(urb))
## ------------------------------------------------------------------------
caop<-sf::st_read(dsn=file.path(getwd(),"datasets","Cont_AAD_CAOP2015.shp"))
caop.proj<-sf::st_transform(caop, crs=sf::st_crs(flam))
regiao<-sf::st_union(caop.proj[caop.proj$Concelho=="GUIMARÃES",])
## ------------------------------------------------------------------------
its<-st_intersects(flam,regiao,sparse=FALSE)
flams<-flam[its[,1],]
its<-st_intersects(urb,regiao,sparse=FALSE)
urbs<-urb[its[,1],]
## ------------------------------------------------------------------------
mapview(flams,col.regions="green")+mapview(urbs,col.regions="gray")
## ------------------------------------------------------------------------
system.time({int<-st_intersection(flams,urbs)}) # 90.74 for regiao==Guimarães
mapview(int,color="red")+mapview(flams,col.regions="green")+mapview(urbs,col.regions="gray")
table(st_geometry_type(int))
int[st_geometry_type(int)=="MULTILINESTRING",]
## ------------------------------------------------------------------------
table(st_geometry_type(int))
int[st_geometry_type(int)=="MULTILINESTRING",]
## ------------------------------------------------------------------------
xyflam<-st_coordinates(flams)
xyurb<-st_coordinates(urbs)
head(xyflam,3)
## ------------------------------------------------------------------------
nrow(flams) # number of features
range(xyflam[,"L2"]) # L2 indicates the index of the feature
## ------------------------------------------------------------------------
xyurb<-cbind(idx=1:nrow(xyurb),xyurb) # idx keeps the index of each original urb vertex
head(xyurb,3)
## ------------------------------------------------------------------------
tol<-3 # tolerance
knn<-nn2(data=xyflam[,c("X","Y")],query=xyurb[,c("X","Y")],k=1,searchtype = "radius",radius=tol)
tail(knn$nn.idx,3)
xyint<-xyurb[knn$nn.idx>0,]
## ------------------------------------------------------------------------
sfint<-st_as_sf(as.data.frame(xyint), coords = c("X","Y") , crs = st_crs(urb)$proj4string)
mapview(sfint,color="red")+mapview(flams,col.regions="green")+mapview(urbs,col.regions="gray")
## ------------------------------------------------------------------------
jump <- diff(sfint$idx)>1 # indicates if the vertex index "jumps" more than 1 position
sfint<-cbind(sfint,segment=cumsum(c(0,jump))) # adds columns lineidx
## ------------------------------------------------------------------------
fastint<-sfint %>% group_by(L2,L1,segment) %>% summarize(do_union=FALSE) %>% st_cast("LINESTRING")
mapview(sfint,color="blue")+mapview(fastint,color="red")+mapview(flams,col.regions="green")+mapview(urbs,col.regions="gray")
## ---- eval=TRUE----------------------------------------------------------
system.time({
xyflam<-st_coordinates(flam)
xyurb<-st_coordinates(urb)
xyurb<-cbind(idx=1:nrow(xyurb),xyurb)
tol<-3 # tolerance to consider vertices coincident
knn<-nn2(data=xyflam[,c("X","Y")],query=xyurb[,c("X","Y")],k=1,searchtype = "radius",radius=tol)
xyint<-xyurb[knn$nn.idx>0,]
sfint<-st_as_sf(as.data.frame(xyint), coords = c("X","Y") , crs = st_crs(urb)$proj4string)
jump <- diff(sfint$idx)>1
sfint<-cbind(sfint,segment=cumsum(c(0,jump)))
fastint<-sfint %>% group_by(L2,L1,lineidx) %>% summarize(do_union=FALSE) %>% st_cast("LINESTRING")
})