Load libraries and read in the tracking data

library(sf)
library(rgdal)
library(mapview)
library(leaflet.extras)
library(lubridate)
library(dplyr)
library(raster)
library(ggplot2)
library(concaveman)

mapviewOptions(basemaps =c("Esri.WorldImagery",  "OpenTopoMap", "OpenStreetMap.BlackAndWhite", "OpenStreetMap.HOT"))

pnts<-read.csv("/home/chris/Behav_clean_reproj.csv")
pnts<-pnts[!is.na(pnts$x_UTM),]
pnts$Date<-dmy(pnts$Date)
pnts$Date_time<-dmy_hm(pnts$Date_time)
pnts$hr<-hour(pnts$Date_time)
coordinates(pnts) = ~x_UTM+y_UTM  # Set coordinates
pnts<-st_as_sf(pnts)  # Make sf object
st_crs(pnts) <- 32647## Set CRS
pnts_3857<-st_transform(pnts,3857) ## The rasters are in Google mercator, which is useful for visualisation and reasonably accurate for measurements near the equator, although not perfect. This could be corrected by reprojection to UTM. Anyway, we'll transform to this for now.

Processing digital surface raster

Aggregating the dsm to 25m squares and calculating mean, min, max, sd and interquartile ranges.

dsm<-raster("/home/rstudio/shiny/data_loggers/data_loggers2/insol_tiffs/dsm.tif")

dsm_sd<-raster::aggregate(dsm,fact=25,fun=sd,na.rm=TRUE)
dsm_mean<-raster::aggregate(dsm,fact=25,fun=mean,na.rm=TRUE)
dsm_max<-raster::aggregate(dsm,fact=25,fun=max,na.rm=TRUE)
dsm_min<-raster::aggregate(dsm,fact=25,fun=min,na.rm=TRUE)

### Interquartile range
iqr<-function(x, ...){
  x1<-quantile(x,0.75,na.rm=TRUE)
  x2 <- quantile(x,0.25,na.rm=TRUE)
  x1-x2
}
dsm_iqr<-raster::aggregate(dsm,fact=25,fun=iqr,na.rm=TRUE)

dsm_brick<-brick(dsm_mean,dsm_sd,dsm_iqr,dsm_max,dsm_min)
names(dsm_brick)<-c("dsm_mean","dsm_sd","dsm_iqr","dsm_max","dsm_min")

Plot the rasters as a mapview

m<-mapview(dsm_brick)
m@map %>% addFullscreenControl()

Extract values

pnts_sp<-as(pnts_3857, "Spatial") ## Change to sp

d<-cbind(pnts_3857,raster::extract(dsm_brick,pnts_sp))

Look at interquartile range (canopy heterogeneity) and feeding over time.

d1<-filter(d,Activity%in%c("eyl","fruit"))
g1<-ggplot(d1,aes(x=hr,y=dsm_iqr))
g1 +geom_point() + geom_smooth() +facet_wrap("Activity")

Nothing much to make of this.

Tree model

If there is any possibility of discriminating any acativity using these variables something should come out on a tree model, which could then be refined.

library(rpart)

mod<-rpart(Activity~dsm_mean+dsm_sd +dsm_min + dsm_max + hr,data=d)

library(partykit)
plot(as.party(mod), digits=3)

Nothing here either, apart from the known differences at the start and end of the day when moving to the sleep trees.

Sorry, but there is just not much to see.