Testing the display of AgTech cow GPS data to demonstrate spatial data usage

Using R, sf, ggplot, ggmap, mapbox, and leaflet for static and dynamic map plots, demonstrating workflow data from Snowflake, repo in Github, and published on RStudio connect at: https://connect.dairynz.co.nz/connect/#/apps/6f1136bf-9f62-4873-b91f-1caadfd904f0

This code currently lives in gps-platemeter repo (see Paul Edwards) Borrowing some code from Cow GPS repo (see Wayne Hoffman)

knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
library(tidyverse)
library(sf)
library(lubridate)

library(ggmap)

library(leaflet)
# Using Mapbox for leaflet satellite background
# https://walker-data.com/mapboxapi/
#install.packages("mapboxapi") #once only
library(mapboxapi)
# once only, install token to environment
#mb_access_token("pk.eyJ1IjoibWFya2JuZWFsIiwiYSI6ImNrYW05ZGRqNDFiNnEyeWw2eGQyeDlxaWEifQ.jPGwcOvKFcC08HwbkwbztQ", install = TRUE)

# Mapbox token - Mark's, for colleagues use, do not abuse or set up for high volume usage!
#pk.eyJ1IjoibWFya2JuZWFsIiwiYSI6ImNrYW05ZGRqNDFiNnEyeWw2eGQyeDlxaWEifQ.jPGwcOvKFcC08HwbkwbztQ

library(DBI)
library(odbc)
library(dbplyr)

Get gps data from Agtech devices in Snowflake

And then wrangle and filter to a neat section of the data.

driver    <- "dnz-economics-sf" #The name in OBDC 64 bit connection, a user DSN that points to SnowflakeDSIIDriver
database  <- "SANDPIT_01"
#database  <- "SANDPIT" #for testing when i had the wrong role
schema    <- "AGTECH360"
warehouse <- "COMPUTE_WH"
con <- DBI::dbConnect(odbc::odbc(), 
                      driver,
                      database=database,
                      schema=schema,
                      warehouse=warehouse
                      )
# 
# dbListFields(con, "ANIMAL_LOCATIONS")
# current_locations <- dbReadTable(con, "ANIMAL_LOCATIONS")

# dbListTables(con)
dbListFields(con, "ANIMAL_LOCATIONS_HISTORIC")
## [1] "ANIMALID"  "EID"       "VISUALTAG" "DATE"      "TAGNAME"   "VALUE"    
## [7] "ALARM"     "LONGITUDE" "LATITUDE"
# dbExistsTable(con, "ANIMAL_LOCATIONS_HISTORIC")
locations <- dbReadTable(con, "ANIMAL_LOCATIONS_HISTORIC")

locations <- locations %>%
  filter(!is.na(LONGITUDE)) %>% 
  filter(LONGITUDE != 0) %>% 
  filter(!is.na(LATITUDE)) %>% 
  filter(LATITUDE != 0) %>% 
  mutate(LONGITUDE = as.numeric(LONGITUDE),
         LATITUDE = as.numeric(LATITUDE),
         ANIMALID = as.numeric(str_sub(ANIMALID,start =  2,end = -2)))
# summary(locations)

# Working out why simple as.numeric didn't work - It was because string included " at start and end.
# str_length(locations$ANIMALID[1])
# str_length(noquote(locations$ANIMALID))
# str_sub = str_sub(locations$ANIMALID[1],start =  2,end = -2)
# str_length(str_sub)


agtech_location_data <- locations # %>%
    # dplyr::filter(TAGNAME == "Location") %>%
    # collect()

# agtech_location_data

#lets fix the date data
# str_length(agtech_location_data$DATE)[1]
# str_sub(agtech_location_data$DATE[1],start =  2,end = -2)
agtech_location_data <- agtech_location_data %>%
    mutate(date_time = as_datetime(DATE)) %>% # str_sub(DATE,start =  2,end = -2))
    mutate(NZTime = with_tz(date_time,"Pacific/Auckland")) %>%
    mutate(TAGNAME = "agtech")
# class(agtech_location_data$date_time)

## Include tag serial for consistency previous work

ANIMALID <-c(6258, 6259, 6262, 6265, 6264, 6263, 6256, 6260, 6261, 6266, 6257, 6267, 0002)

##AnimalID  0002 dummy ID believed to relate to Tag Serial 2555 (broken)

TAG_SERIAL <-c(2578, 2579, 2586, 2590, 2583, 2589, 2575, 2581, 2584, 2576, 2577, 2503, 2555)

Agtech_Tags <- data.frame(TAG_SERIAL, ANIMALID)

# Agtech_Tags

agtech_location_data <- full_join(agtech_location_data, Agtech_Tags)
#by = c("AnimalID")

#Filtered for only Bradshaw data

AgTech_Bradshaw_Data <- agtech_location_data %>%
    #mutate(TAG_SERIAL = as.factor(TAG_SERIAL)) %>%
    filter(NZTime >="2021-04-05 06:00:00" & NZTime <="2021-04-09 06:00:00") %>%
    mutate(time = format(NZTime, "%H:%M:%S")) %>%
    mutate(Day = day(NZTime), month (NZTime)) %>%
    mutate(farm_day = ifelse(time > "00:00:00" & time <"06:00:00", Day -1, Day)) %>%
    mutate(grazing_break = ifelse(time > "06:00:00" & time < "15:00:00", 'AM', 'PM')) %>%
    mutate(Cal_day = paste0(farm_day, " ", 'April 2021')) %>%
    dplyr::select(ANIMALID, TAGNAME, LATITUDE, LONGITUDE, DATE,
                   NZTime, time, Day, farm_day, Cal_day, grazing_break)  #Data available for 4 weeks



locations_sf <- st_as_sf(AgTech_Bradshaw_Data, coords = c("LONGITUDE", "LATITUDE"), crs = 4326)

p <- ggplot(locations_sf) +
  geom_sf(aes(colour = as.factor(ANIMALID))) +
  facet_wrap(~ farm_day)
p

ggsave("plot of agtech locations.png", p)

Get paddock boundaries from shapefile

# Get boundaries and include in plots

paddock_gis <- st_read("jnkn_geo.shp", quiet=TRUE)
#st_crs(paddock_gis)
paddock_gis$paddock <- paddock_gis$field_data %>%
    str_replace_all("[^[:alnum:]]", " ") %>%
    str_replace("^.+name ", "") %>%
    str_trim()

ggplot()+
  geom_sf(data = paddock_gis, aes(colour = paddock))

Including ggplots

You can plot in ggplot with a Google satellite background:

Bradshaw_map <- get_map(location=c(lon = 172.552, lat = -43.346), zoom=14, maptype = 'hybrid')

ggmap(Bradshaw_map)

ggmap(Bradshaw_map)+
  geom_sf(data = paddock_gis, aes(colour = paddock), inherit.aes = FALSE) +  
  # The inherit.aes = FALSE is important when mixing ggmap and sf objects
  geom_sf(data = locations_sf, aes(colour = as.factor(ANIMALID)), inherit.aes = FALSE)+
  facet_wrap(~ farm_day)

ggsave("AgtechTestPeriodBradshaw.jpeg")

Including Leaflet Plots

You can also embed leaflet plot for zoom and scroll of plots. However, for true interactivity (eg map changes on date or time selection), you (probably?) need to go to Shiny.

#Bradshaws
#-43.34634731910686, 172.55232788540212

m <- leaflet(locations_sf) %>% 
  setView(lng = 172.552, lat = -43.346, zoom = 14)

# m %>% addTiles() #using open street maps

# marker options (points from sf object are markers in leaflet terminology)
# lots of marker options!
# https://rstudio.github.io/leaflet/markers.html

m %>%
  addMarkers(options = markerOptions()) %>% 
  addPolygons(data = paddock_gis, color = "blue") %>% 
  addMapboxTiles(style_id = "satellite-streets-v11",
                 username = "mapbox") 
# get style id from here:
# https://docs.mapbox.com/api/maps/styles/

Note that the echo = FALSE parameter could be added to the code chunk to prevent printing of the R code that generated the plot.