## ------------------------------------------------------------------------
library(raster)
library(rgdal) # read/write with GDAL
library(sf) # vector data (replaces sp)
library(mapview) #mapView
library(RANN) # nn2
library(data.table) # fread
library(velox) # alternative to raster
library(fasterize) # fast resampling
library(viridisLite) # colors (used in mapview)
## ---- echo=FALSE---------------------------------------------------------
#knitr::knit_exit()
## ---- echo=FALSE---------------------------------------------------------
setwd("\\\\pinea\\mlc\\Aulas\\CURSOS_R\\isa-setembro-2018-grandes-cdg\\aulas")
library(knitr)
## ------------------------------------------------------------------------
system.time({flam.sp<-rgdal::readOGR(dsn=file.path(getwd(),"datasets"),layer="Flam2015v4",verbose = FALSE)}) # 21.27
flam.sp
## ------------------------------------------------------------------------
system.time({flam.sf<-sf::st_read(dsn=file.path(getwd(),"datasets","Flam2015v4.shp"),quiet=TRUE)}) # 4.32
print(flam.sf,n=5)
## ----eval=FALSE----------------------------------------------------------
## as(x, "Spatial") # x is a sf object
## st_as_sf(x) # x is a sp object
## ------------------------------------------------------------------------
areas<-st_area(flam.sf)
bigflam<-flam.sf[which.max(areas),]
## ---- eval=TRUE----------------------------------------------------------
mapview(bigflam)
## ---- eval=TRUE----------------------------------------------------------
caop<-sf::st_read(dsn=file.path(getwd(),"datasets","Cont_AAD_CAOP2015.shp"),quiet=TRUE)
print(caop,n=3)
# reproject caop to the CRS of flam.sf
caop.proj<-sf::st_transform(caop, crs=sf::st_crs(flam.sf))
concelho<-caop.proj[caop.proj$Concelho=="MONCHIQUE",]
flam.concelho<-sf::st_intersection(flam.sf,concelho)
mapview(concelho,col.regions="blue")+mapview(flam.concelho,col.regions="green")
## ------------------------------------------------------------------------
system.time({s<-raster::stack(as.list(file.path(getwd(),"datasets",
c("LC82030332014151LGN00_sr_band3.tif",
"LC82030332014151LGN00_sr_band4.tif",
"LC82030332014151LGN00_sr_band5.tif"))))})
inMemory(s) #FALSE
system.time({b<-raster::brick(s)}) # takes some time
inMemory(b) #TRUE
writeRaster(b,"b.tif",overwrite=TRUE) # creates multilayer file
b2<-brick("b.tif")
inMemory(b2) #FALSE
b
## ------------------------------------------------------------------------
names(b)<-c("b3","b4","b5")
b[b<=0 | b>10000]<-NA # non valid values
ndvi<-(b$b5-b$b4)/(b$b5+b$b4)
writeRaster(ndvi,file="ndvi.tif",overwrite=TRUE)
## ------------------------------------------------------------------------
system.time({s<-raster::stack(as.list(file.path(getwd(),"datasets",c("LC82030332014151LGN00_sr_band3.tif","LC82030332014151LGN00_sr_band4.tif", "LC82030332014151LGN00_sr_band5.tif"))))})
vxs<-velox(s)
names(vxs$rasterbands)<-c("b3","b4","b5") # each band is a matrix
## ------------------------------------------------------------------------
v<-function(x) {x[x<=0 | x>10000]<-NA; return(x)}
b5<-v(vxs$rasterbands$b5) # assign NA to non valid values
b4<-v(vxs$rasterbands$b4) # idem
vals<-(b5-b4)/(b5+b4)
vxndvi<-velox(vals,extent=vxs$extent,res=vxs$res,crs=vxs$crs)
## ------------------------------------------------------------------------
r<-vxndvi$as.RasterLayer(band=1)
mapview(r)
## ------------------------------------------------------------------------
# extract values
ndvi<-raster("ndvi.tif")
ext<-as.vector(extent(ndvi))
N<-10000; locs<-cbind(runif(N,min=ext[1],max=ext[2]),runif(N,min=ext[3],max=ext[4]))
system.time({x<-raster::extract(ndvi,y=locs)}) # 4.94 for N=10000
sf.locs<-st_as_sf(data.frame(x=locs[,1],y=locs[,2]), coords = c("x","y") , crs = as.character(crs(ndvi)))
system.time({y<-vxndvi$extract_points(sf.locs)})
n<-4;setequal(round(x,n),round(y,n))
## ------------------------------------------------------------------------
# read coordinates
system.time({xy<-coordinates(ndvi)}) # 1.34
system.time({xy<-vxndvi$getCoordinates()}) # 0.44
system.time({x<-raster::crop(ndvi,0.5*extent(ndvi))}) # 2.66
system.time({vxndvi$crop(0.5*extent(ndvi))}) # 0.16
## ------------------------------------------------------------------------
vxndvi$extent # after crop
# original version
vxndvi<-velox(vals,extent=vxs$extent,res=vxs$res,crs=vxs$crs)
vxndvi$extent #
# make a copy that is not referenced to vxndvi
vxcopy<-vxndvi$copy()
## ------------------------------------------------------------------------
cores<-c("red","yellow","green","blue")
cores<-colorRampPalette(c(rgb(0,0,1,alpha=1), rgb(0,0,1,alpha=0)), alpha = TRUE)(8)
cores<-viridisLite::inferno(n=10,alpha=1,begin=0.2,end=1,direction= -1)
pie(rep(1,length(cores)),col=cores)
## ------------------------------------------------------------------------
cores<-colorRampPalette(c(rgb(1,1,0,0.5), rgb(0,1,0,0.5)), alpha = TRUE)(10) # any number >=4
mapview(ndvi,col.regions=cores,at=c(0.1,.3,.5,.7,1))
## ------------------------------------------------------------------------
#system.time({df<-read.csv(file=file.path(getwd(),"datasets","viirs-fires-2017-world.csv"))}) # 148.41
system.time({af<-data.table::fread(file.path(getwd(),"datasets","viirs-fires-2017-world.csv"))}) # 14.34
af
## ------------------------------------------------------------------------
afpt<-af[latitude>36 & latitude<43 & longitude < -6 & longitude > -10 ]
data.table::fwrite(afpt,file=file.path(getwd(),"datasets","viirs-fires-2017-portugal.csv"))
## ------------------------------------------------------------------------
af[,footprint:=scan*track] # create new variable "footprint"
af[order(frp,decreasing=TRUE)] # order table by frp value
af[,sum(frp),by=acq_date] # returns table with 365 rows
af[,s:=sum(frp),by=acq_date] # adds new column to af
af[,.N,by=confidence] # a bit faster than >table(af$confidence)
af[,.(confidence,frp)] # select columns
## ------------------------------------------------------------------------
DT<-data.table(x=1:3,y=c(4,8,12))
newDT<-DT
newDT[,x:=NULL]
DT # is modified by reference
## ------------------------------------------------------------------------
DT<-data.table(x=1:3,y=c(4,8,12))
newDT<-copy(DT)
newDT[,x:=NULL]
DT # is not modified
## ------------------------------------------------------------------------
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)
## ------------------------------------------------------------------------
afaf<-af[latitude > -35 & latitude < 0 & longitude > 0 & longitude < 40] # S Africa
latlon<-CRS("+proj=longlat +datum=WGS84")
y<-raster(crs=latlon,resolution=0.01,xmn=0,xmx=40,ymn=-35,ymx=0)
system.time({r<-rasterize(x=data.frame(lon=afaf$longitude,lat=afaf$latitude),y=y,field=afaf$frp,fun=max)}) # 124.16
mapview(r)
## ------------------------------------------------------------------------
system.time({flam.sf<-sf::st_read(dsn=file.path(getwd(),"datasets","Flam2015v4.shp"))})
ext<-st_bbox(flam.sf)
RES<-50
r<-raster(xmn=ext$xmin,xmx=ext$xmax,ymn=ext$ymin,ymx=ext$ymax,crs=st_crs(flam.sf)$proj4string,resolution=RES)
rflam<-fasterize(flam.sf,r,field="FID_COS")
#
caop<-sf::st_read(dsn=file.path(getwd(),"datasets","Cont_AAD_CAOP2015.shp"),quiet=TRUE)
caop.proj<-sf::st_transform(caop, crs=sf::st_crs(flam.sf))
concelho<-caop.proj[caop.proj$Concelho=="MONCHIQUE",]
mapview(rflam)+mapview(sf::st_crop(flam.sf,sf::st_bbox(concelho)))