pasture9.rmd

---
title: "Pasture Potential"
runtime: shiny
output: html_document
---

<!-- Simon Woodward, DairyNZ, 2018 -->
<!-- now use joined.rds data, including elevation nitrogen etc -->
<!-- # Note: Unlike Shiny Apps, Interactive R Markdown Documents are do not require ui and server,  -->
<!-- #       the whole document is treated as a server. -->
<!-- # You can't guarantee what order the reactive elements will execute I think -->


```{r setup, include=FALSE} 
knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(leaflet) # interactive map
library(ggmap) 
library(cowplot)
library(geosphere)
library(proj4) # warning masked by rgdal
library(KernSmooth) # for contouring
library(sp) # spatial polygons
library(quantreg)
#library(broom) # used to convert spatial data to tibble
library(rgdal) 

#### these objects are available to all sessions ####
# this is good place to put global data or utility functions
# you can change them (for every session) but you have to use the <<- operator

# read multiyear data and tidy up 
data_all <- readRDS("joined.rds") %>%
  mutate(pasture_eaten=as.numeric(pasture_and_crop_eaten_t_dm_ha),
         region=as.factor(region),
         season=as.factor(season),
         slope=slope.x,
         aspect=aspect.x,
         nitrogen_applied=as.numeric(nitrogen_applied_for_year_l2),
         topo=topo_position_index.y,
         soil=as.character(nzsc_order.x))  
data_all <- data_all %>%
  select(farm_number, pasture_eaten, region, season, long, lat, soil, elev, 
         slope, aspect, topo, nitrogen_applied) 
temp <- nrow(data_all)
data_all <- data_all %>% drop_na() # drops a lot of rows!
cat(file=stderr(), paste("Rows dropped with missing data =", temp-nrow(data_all), "\n"))
cat(file=stderr(), paste("Rows remaining =", nrow(data_all), "\n"))

summary(data_all)
data_all <- data_all %>%
  mutate(elev_fact=factor(x=case_when(elev<=50 ~ "Low (0-50m)", 
                                      elev<=200 ~ "Middle (50-200m)", 
                                     TRUE ~ "High (200m+)"),
                          levels=c("Low (0-50m)", "Middle (50-200m)", "High (200m+)")))

# create season and soil list
season_all <- as.list(sort(unique(as.character(data_all$season))))
names(season_all) <- season_all
soil_all <- as.list(sort(unique(data_all$soil)))
names(soil_all) <- soil_all
names(soil_all)[soil_all=="L"] <- "Allophanic"
names(soil_all)[soil_all=="A"] <- "Anthropic"
names(soil_all)[soil_all=="B"] <- "Brown"
names(soil_all)[soil_all=="Z"] <- "Podzol"
names(soil_all)[soil_all=="M"] <- "Pumice"
names(soil_all)[soil_all=="W"] <- "Raw"
names(soil_all)[soil_all=="G"] <- "Gley"
names(soil_all)[soil_all=="N"] <- "Granular"
names(soil_all)[soil_all=="E"] <- "Melanic"
names(soil_all)[soil_all=="R"] <- "Recent"
names(soil_all)[soil_all=="S"] <- "Semiarid"
names(soil_all)[soil_all=="U"] <- "Ultic"
names(soil_all)[soil_all=="O"] <- "Organic"
names(soil_all)[soil_all=="X"] <- "Oxidic"
names(soil_all)[soil_all=="P"] <- "Pallic"
elev_all <- levels(data_all$elev_fact)
names(elev_all) <- elev_all

# create data point map/contours for plotting on leaflet
# https://gis.stackexchange.com/questions/168886/r-how-to-build-heatmap-with-the-leaflet-package
data_pts <- unique(data_all[c("long", "lat")])
kde <- bkde2D(data.matrix(data_pts), bandwidth=c(0.1, 0.1), gridsize=c(100,100))
CL <- contourLines(kde$x1 , kde$x2 , kde$fhat) # contour lines (list)
LEVS <- as.factor(sapply(CL, `[[`, "level")) # contour levels (vector)
NLEV <- length(levels(LEVS)) # number of levels 
pgons <- lapply(1:length(CL), function(i) # convert to polygons (ID=i is actually the level)
    Polygons(list(Polygon(cbind(CL[[i]]$x, CL[[i]]$y))), ID=i))
spgons = SpatialPolygons(pgons)
spgon_cols <- topo.colors(NLEV, NULL)[LEVS]

# convert spatialPolygons to data frame for use in ggplot
#spgonsdf <- broom::tidy(spgons, region=ID)

# calculate NZTM2000 coordinates for farm locations
proj4string <- "+proj=tmerc +lat_0=0.0 +lon_0=173.0 +k=0.9996 +x_0=1600000.0 +y_0=10000000.0 +datum=WGS84 +units=m"
nzgd <- data.matrix(data_all[,c("long", "lat")])
nztm <- proj4::project(xy=nzgd, proj=proj4string)
temp <- proj4::project(xy=nzgd, proj=proj4string, inverse=TRUE)
data_all$east <- nztm[,1]
data_all$north <- nztm[,2]

# define some constants 
trim <- 0.0 # rqss fails near tails if insufficient data
probs <- seq(trim, 1-trim, 0.02) # probabilities for sampcdf
nprobs <- length(probs)
windows <- c(60,40,20)
nmin <- 4L # minimum number of farms in a window for analysis

# gets a list of default ggplot colours
gg_colour_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
window_cols <- gg_colour_hue(length(windows))

# normal list with same elements, useful for testing
# my <- list(name="You Are Here", 
#                      # long=172.833333, lat=-41.5, # Nelson
#                      long=175.619105, lat=-40.386396, # Massey
#                      # long=174.865530, lat=-41.259256, # Wellington
#                      # long=175.352116, lat=-37.781841, # DairyNZ
#                      distortion=1, 
#                      east=NA, north=NA, 
#                      season_here=list(NA), 
#                      season=list(NA),
#                      soil_here=list(NA), 
#                      soil=list(NA),
#                      elev_here=list(NA), 
#                      elev=list(NA),
#                      data=list(NA),
#                      breaks=list(NA),
#                      recalc=0L
#                      )

```

Pasture is a fundamental component of profitable dairy systems. In general, the more pasture you can grow and feed to your animals, the more profitable your system will be. This tool allows you to compare your pasture produced and consumed with pasture produced and consumed on other farms in your region. This can indicate whether there is potential for you to increase pasture production and intake, and hence profitability.

# Where are You?

Please select which season of data to analyse, and click your location on the map. The contours show the availability of data.

```{r input_pane, echo=FALSE, eval=TRUE}

#### initialise session ####
# useful trigger object to allow manual control (not used yet)
# https://www.r-bloggers.com/dynamically-generated-shiny-ui/
makeReactiveTrigger <- function() {
   rv <- reactiveValues(a = 0)
   list(
     depend = function() {
       rv$a
       invisible()
     },
     trigger = function() {
       rv$a <- isolate(rv$a + 1)
     }
   )
}
# usage
# my_trigger <- makeReactiveTrigger() # create trigger object
# my_trigger$depends() # put this where you want a dependency
# my_trigger$trigger() # put this where you want to trigger the dependency

# collect info about the current location and selections
# these variables are available in the reactive conext
# "my" is not itself a reactive object, it's a list of reactive objects
my <- reactiveValues(name="You Are Here", 
                     # long=172.833333, lat=-41.5, # Nelson
                     long=175.619105, lat=-40.386396, # Massey
                     # long=174.865530, lat=-41.259256, # Wellington
                     # long=175.352116, lat=-37.781841, # DairyNZ
                     distortion=1, 
                     east=NA, north=NA, 
                     season_here=list(NA), 
                     season_default=list(NA), 
                     season=list(NA),
                     soil_here=list(NA), 
                     soil_default=list(NA), 
                     soil=list(NA),
                     elev_here=list(NA), 
                     elev_default=list(NA), 
                     elev=list(NA),
                     data=list(NA),
                     breaks=list(NA),
                     recalc=0L
                     )

# set leaflet info
v <- reactiveValues(zoom=5, minzoom=5, maxzoom=15, long=NA, lat=NA) 

# initialise map centre
isolate({
  v$long <- my$long
  v$lat <- my$lat
})

#### ui elements ### 

# titlePanel("Where are You?")

output$season_selector <- renderUI({
  cat(file=stderr(), paste("render season selector"), "\n")
  selectInput("season", h4("Season?"), my$season_here, selected=my$season_default)
})

output$soil_selector <- renderUI({
  cat(file=stderr(), paste("render soil selector"), "\n")
  selectInput("soil", h4("Soils?"), my$soil_here, selected=my$soil_default, 
              selectize=FALSE, multiple=TRUE)
})

output$elev_selector <- renderUI({
  cat(file=stderr(), paste("render elev selector"), "\n")
  selectInput("elev", h4("Elevation?"), my$elev_here, selected=my$elev_default,
              selectize=FALSE, multiple=TRUE)
})

# Make your initial map
# https://stackoverflow.com/questions/34348737/r-leaflet-how-to-click-on-map-and-add-a-circle
output$map <- renderLeaflet({
  cat(file=stderr(), paste("render leaflet"), "\n")
  
    isolate({ # prevent redraw if arguments change

      leaflet(spgons, options=leafletOptions(minZoom=v$minzoom, maxZoom=v$maxzoom)) %>%
        setView(v$long, v$lat, zoom=v$zoom)  %>%
        addTiles() %>% # default map
        addPolygons(data=spgons, color=spgon_cols, weight=0, options=pathOptions(clickable=FALSE)) %>%
        addMarkers(my$long, my$lat, "layer1", options=pathOptions(clickable=FALSE)) %>%
        addCircles(my$long, my$lat, layerId=as.character(windows), radius=windows*1000, 
                   color=window_cols, weight=4, fill=NA, options=pathOptions(clickable=FALSE)) 
    })

  }) # end renderLeaflet

#### ui layout ####

sidebarLayout(
  
  sidebarPanel(
    cat(file=stderr(), paste("render sidebar"), "\n"),
    uiOutput("season_selector"), # this control is created in the server
    uiOutput("soil_selector"), # this control is created in the server
    uiOutput("elev_selector") # this control is created in the server
    # plotOutput("plot6", height=150) # soil histogram
    # actionButton("go", "Go!") # manual recalc button
  ), # end sidebarPanel
  
  mainPanel(
    cat(file=stderr(), paste("render main panel"), "\n"),
    leafletOutput("map"),
    tags$head(tags$style(
			"#map {
			cursor: pointer;
			}")) 
  ), # end mainPanel
  
  position="right"
  
) # end sidebarLayout


#### react to mouse clicks ####
# see also https://rstudio.github.io/leaflet/shiny.html
observeEvent(input$map_click, {
  
  cat(file=stderr(), "\n")
  cat(file=stderr(), paste("observed map_click"), "\n")    

  click <- input$map_click
  my$long <- click$lng
  my$lat <- click$lat

  # mark map
  leafletProxy("map", deferUntilFlush=FALSE) %>%
    addMarkers(my$long, my$lat, "layer1", options=pathOptions(clickable=FALSE)) %>%
    addCircles(my$long, my$lat, layerId=as.character(windows), radius=windows*1000, 
               color=window_cols, weight=4, fill=NA, options=pathOptions(clickable=FALSE)) 
  
}) # end observe mouse click

#### react to change of location ####
observeEvent(c(my$long, my$lat), {

  cat(file=stderr(), "\n")
  cat(file=stderr(), paste("change of location =", my$long, my$lat), "\n")

  # calculate aspect ratio near my farm
  nzgd <- data.matrix(tibble(long=c(my$long, my$long, my$long-0.5, my$long+0.5),
                             lat=c(my$lat-0.5, my$lat+0.5, my$lat, my$lat)))
  nztm <- proj4::project(xy=nzgd, proj=proj4string)
  my$distortion <- (max(nztm[,2])-min(nztm[,2]))/(max(nztm[,1])-min(nztm[,1]))
  
  # location for map centre
  nzgd <- data.matrix(c(my$long, my$lat))
  nztm <- proj4::project(xy=nzgd, proj=proj4string)
  my$east <- nztm[,1]
  my$north <- nztm[,2]

  # filter by distance #
  # we need to use rowwise()  because distm is not vectorised, I think, although rowwise() is deprecated
  # http://www.expressivecode.org/2014/12/17/mutating-using-functions-in-dplyr/
  data <- data_all %>% 
    rowwise() %>% 
    mutate(dist = distm(c(my$long, my$lat), c(long, lat), fun=distHaversine), # this needs rowwise()
           dist = dist/1000) %>% # km
    filter(dist < max(windows))

  my$data <- data
  
  # calculate width of histogram for region
  my$breaks <- seq(floor(min(data$pasture_eaten)), ceiling(max(data$pasture_eaten)), 1)

  # what season are available here   
  i <- sort(unique(as.character(data$season)))
  my$season_here <- season_all[match(i, season_all)]
  n <- unlist(map(my$season_here, function(u) sum(u==data$season)))
  names(my$season_here) <- paste(names(my$season_here), " (", n, " Farms)", sep="")
  cat(file=stderr(), paste("my$season_here =", length(my$season_here)), "\n")
  cat(file=stderr(), paste(names(my$season_here)), "\n")

  # reset selections
  my$season_default <- tail(my$season_here, 1)
  my$season <- my$season_default
  my$soil <- list(NA)
  my$elev <- list(NA)
  
}) # end reaction to location changing

observeEvent(input$season, {
  req(input$season, input$season!="NA")
  if (my$season != input$season) {
    my$season <- input$season
    my$soil <- list(NA)
    my$elev <- list(NA)
  }
})

#### react to change of season ####
observeEvent(my$season, {

  cat(file=stderr(), "\n")
  cat(file=stderr(), paste("change of my$season =", my$season), "\n")
  req(my$season, my$season!="NA")

  data <- my$data %>%
    filter(season == my$season)
  
  # make soil list
  i <- sort(unique(data$soil))
  my$soil_here <- soil_all[match(i, soil_all)]
  n <- unlist(map(my$soil_here, function(u) sum(u==data$soil)))
  names(my$soil_here) <- paste(names(my$soil_here), " (", n, " Farms)", sep="")
  cat(file=stderr(), paste("my$soil_here =", length(my$soil_here)), "\n")
  cat(file=stderr(), paste(names(my$soil_here)), "\n")

  # reset selections
  my$soil_default <- my$soil_here
  my$soil <- my$soil_default 
  my$elev <- list(NA)

}) # end reaction to season changing

observeEvent(input$soil, {
  req(input$soil, input$soil!="NA")
  if (any(my$soil != input$soil)) {
    my$soil <- input$soil
    my$elev <- list(NA)
  }
})

#### react to change of soil ####
observeEvent(my$soil, {
  
  cat(file=stderr(), "\n")
  cat(file=stderr(), paste("change of my$soil = "), paste(my$soil), "\n")
  req(my$soil, my$soil!="NA")
  
  data <- my$data %>%
    filter(season == my$season) %>%
    filter(soil %in% my$soil)

  # make elev list
  # data <- data_all[1:10,]
  i <- unique(data$elev_fact)
  my$elev_here <- elev_all[which(elev_all %in% i)] # retains sorting
  n <- unlist(map(my$elev_here, function(u) sum(u==data$elev_fact)))
  names(my$elev_here) <- paste(my$elev_here, " (", n, " Farms)", sep="")
  cat(file=stderr(), paste("my$elev_here =", length(my$elev_here)), "\n")
  cat(file=stderr(), paste(names(my$elev_here)), "\n")

  # reset selections
  # my$elev_default <- my$elev_here[which(n==max(n))]
  my$elev_default <- my$elev_here
  my$elev <- my$elev_default 

}) # end reaction to soil changing

observeEvent(input$elev,  {
  req(input$elev, input$elev!="NA")
  if (any(my$elev != input$elev)) {
    my$elev <-  input$elev
  }
})

#### react to change of elevation ####
observeEvent(my$elev,  {

  cat(file=stderr(), "\n")
  cat(file=stderr(), paste("change of my$elev ="), paste(my$elev), "\n")
  req(my$elev, my$elev!="NA")

  my$recalc <- my$recalc + 1L # trigger recalc
  
})

```

<!-- # Your Neighbourhood -->

```{r output_calculations, eval=TRUE, echo=FALSE, warning=TRUE}

#### react to change of inputs ####
calc <- eventReactive(my$recalc,  {

  req(my$season!="NA", my$soil!="NA", my$elev!="NA")
  req(my$long, my$lat, my$season, my$soil, my$elev)
  
  cat(file=stderr(), "\n")
  cat(file=stderr(), paste("analyse for location =", my$long, my$lat), "\n")
  cat(file=stderr(), paste("analyse for season =", my$season), "\n")
  cat(file=stderr(), paste("analyse for soil ="), paste(my$soil), "\n")
  cat(file=stderr(), paste("analyse for elev ="), paste(my$elev), "\n")

  my$name <- paste("You (", season_all[my$season[[1]]], ")", sep="")

  data <-  my$data %>% 
    filter(season == my$season) %>% 
    filter(soil %in% my$soil) %>%
    filter(elev_fact %in% my$elev)
  
  cat(file=stderr(), paste("nrow(data) = ", nrow(data)), "\n")

  # circle function
  circle_fun <- function(centre=c(0,0), r=1, npoints=100){
    tt <- seq(0, 2*pi, length.out=npoints)
    xx <- centre[1] + r * cos(tt)
    yy <- centre[2] + r * sin(tt)
    return(tibble(x=xx, y=yy))
  }    
    
  # prepare empty data frames for loop
  farms <- tibble(x=numeric(), y=numeric(), east=numeric(), north=numeric(), long=numeric(), lat=numeric(),
                  pasture=numeric(), dist=numeric(), window=numeric(), radius=factor())
  sampcdf <- tibble(probs=numeric(), quants=numeric(), radius=factor())
  samppdf <- tibble(pasture=numeric(), window=numeric(), radius=factor(),
                    q=numeric(), qr=numeric(), qrlower=numeric(), qrupper=numeric())
  circles <- tibble(east=numeric(), north=numeric(), radius=factor())
  
  # loop through decreasing window sizes 
  for (window in windows) {

    # select data within window
    data_window <- data %>% filter(dist < window)
    n <- nrow(data_window)
    code <- paste(window," km\n(", format(n, width=3), " Farms)", sep="")

    cat(file=stderr(), paste("window = ", window, " km"), "\n")

    # calculate circle
    circle <- circle_fun(centre=c(my$east, my$north), r=window*1000, npoints=100)
    nztm <- data.matrix(circle[,c("x", "y")])
    nzgd <- proj4::project(xy=nztm, proj=proj4string, inverse=TRUE)
    circle$long <- nzgd[,1]
    circle$lat <- nzgd[,2]

    # save selected farms for plot
    if (n >= 1) {
      farms <- rbind(farms, tibble(east=data_window$east, north=data_window$north,
                                 long=data_window$long, lat=data_window$lat,
                                 pasture=data_window$pasture_eaten,
                                 dist=data_window$dist, window=window, radius=as.factor(code)))
    }
    
    circles <- rbind(circles, tibble(long=circle$long, lat=circle$lat, radius=as.factor(code)))

    # save sample quantiles if enough data to be sensible
    if (n >= nmin) {
      
      # calculate quantile
      qr1 <- rq(formula=pasture_eaten ~ 1, tau=0.9, data=data_window) # linear quantile regression
      se_method <- "boot" # how condience intervals are calculated, some methods more robust
      yqr1<- predict(qr1, tibble(east=my$east, north=my$north), interval="confidence", level=0.95, se=se_method)
      q90 <- quantile(data_window$pasture_eaten, 0.9, type=1) # also calc simple q90
  
      cat(file=stderr(), paste("yqr1 =", yqr1), "\n")
      # cat(file=stderr(), paste("q90 =", q90), "\n") # should be the same

      quants <- quantile(data_window$pasture_eaten, probs=probs, type=8) # see documentation for type=?
      sampcdf <- rbind(sampcdf, tibble(probs=probs, quants=quants, radius=as.factor(code)))
      samppdf <- rbind(samppdf, tibble(pasture=data_window$pasture_eaten, window=window, radius=as.factor(code),
                       q=q90, qr=yqr1[1], qrlower=yqr1[2], qrupper=yqr1[3]))
    } # if n >= nmin

  } # next window size

  # biggest circle
  circle <- circle_fun(centre=c(my$east, my$north), r=max(windows)*1000, npoints=100)
  nztm <- data.matrix(circle[,c("x", "y")])
  nzgd <- proj4::project(xy=nztm, proj=proj4string, inverse=TRUE)
  circle$long <- nzgd[,1]
  circle$lat <- nzgd[,2]

  # return results as function for testing
  #calc <- function() list(data=data, circles=circles, circle=circle, farms=farms, sampcdf=sampcdf, samppdf=samppdf)
  
  # return results in a list
  return(list(data=data, circles=circles, circle=circle, farms=farms, 
              sampcdf=sampcdf, samppdf=samppdf))
  
}) # end reaction to elev chaning, calculation of calc <- list(results)

```

This is the distribution of pasture produced and consumed within various distances of your location. The 90th percentile is also shown.

```{r output_plot, eval=TRUE, echo=FALSE}

# create histograms
output$stacked_histogram <- renderPlot({

  samppdf <- calc()$samppdf # get data for histogram when calc() changes

  isolate({
    
    cat(file=stderr(), paste("render stacked histograms"), "\n")
    
    title_string <- paste("Pasture Eaten near", my$name)
    
    # create empty plot
    stacked_histogram <- ggplot() +
      labs(title=title_string, y="Number of Farms", 
           x="Pasture Eaten (tDM "*ha^-1~y^-1*")", colour="Radius (km)") +
      theme_cowplot() +
      #scale_y_continuous(breaks=c()) + # remove y-scale when too many facets
      panel_border(colour="black") +  
      theme(legend.position="none")
    
    # add histograms to empty plot
    if (nrow(samppdf)>0) { 
      
      breaks <- my$breaks  
      
      # cat(file=stderr(), paste("xlim =", min(breaks), max(breaks)), "\n")
      
      stacked_histogram <- stacked_histogram +
        geom_rect(data=samppdf, mapping=aes(xmin=qrlower, xmax=qrupper, ymin=0, ymax=Inf), fill="lightcyan") +
        geom_histogram(data=samppdf, mapping=aes(x=pasture, colour=radius), fill=NA, size=1.1, binwidth=1) +
        geom_vline(data=samppdf, mapping=aes(xintercept=qr), size=1.5, colour="lightcyan4", alpha=0.2) +
        geom_vline(data=samppdf, mapping=aes(xintercept=q), size=1.5, colour="black") +
        geom_text(data=samppdf, mapping=aes(x=q, y=4, label="90th"), hjust=0, nudge_x=0.1) +
        facet_grid(radius ~ ., as.table=FALSE) + # as.table=FALSE reverses the order
        theme(strip.background=element_blank(), strip.text.y=element_text(angle=0)) +
        scale_x_continuous(breaks=breaks) +
        coord_cartesian(xlim=c(min(breaks),max(breaks)))
  
    } # end add histograms to empty plot
    
  }) # end isolate
  
  stacked_histogram

}) # end renderPlot
  
# show histogram
fluidPage({
    fluidRow(
      plotOutput("stacked_histogram")
    ) # end fluidRow
}) # end fluidPage

```

<!-- this comment was important!!!! -->