R/find_risk.R

In mur4mat/floodshp: Handling ESRI Shapefile Flood Risk

#' Get risk value from flood shapefile
#'
#' @param long Longitude
#' @param lat Latitude
#' @param region_col_name Name of column containing region name from shapefileBounds
#' @param risk_col_name Name of column containing risk value from shapefile
#' @param shapefile_divided Flood shapefile divided via \code{\link{divideNChunks}} function
#' @param shapefile_bounds Shapefile containing teritorial division (i.e. administrative)
#'
#' @return character
#' @export
#'
#' @importFrom sp SpatialPointsDataFrame
#' @importFrom sp CRS
#' @importFrom sp over
#' @examples
#' \dontrun{
#' risk_level <- getRiskLevel(long = -0.08189974, lat = 51.31244, region_col_name = "NAME",
#' risk_col_name = "PROB_4BAND, shapefile_divided = flood_shp_divided, shapefile_bounds = shapefile_bounds)
#'}
getRiskLevel <- function(long, lat, region_col_name, risk_col_name, shapefile_divided, shapefile_bounds){

  # Defining SpatialPointsDataFrame containing queried coordinates
  geocoordinates <- SpatialPointsDataFrame(coords = data.frame(long = as.numeric(long),
                                                               lat = as.numeric(lat)),
                                           data = data.frame('point'),
                                           proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"),
                                           match.ID=TRUE)

  # Check region name of given coordinates
  region_name <- as.character(over(geocoordinates, shapefile_bounds)[[region_col_name]])
  if(is.na(region_name)) stop(simpleError("Given long value out of shapefile boundaries."))

  # Subset region polygon
  region <- shapefile_bounds[which(shapefile_bounds[[region_col_name]]==region_name),]
  region <- spTransform(region, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

  # Subset flood polygons only from selected region
  flood_polygons <- shapefile_divided[[region_name]][region,]

  # Get risk value
  # If given coordinates do not lay over flood polygons, then NA
  tryCatch(
    risk <- as.character(over(geocoordinates, flood_polygons)[[risk_col_name]]),
    error=function(e){
      risk <- NA
    }
  )

  return(risk)
}

#' Generate interactive flood risk map
#'
#' @param long Longitude
#' @param lat Latitude
#' @param region_col_name Name of column containing region name from shapefileBounds
#' @param risk_col_name Name of column containing risk value from shapefile
#' @param shapefile Flood shapefile
#' @param shapefile_divided Flood shapefile divided via \code{\link{divideNChunks}} function
#' @param shapefile_bounds Shapefile containing teritorial division (i.e. administrative)
#'
#' @return A \link{leaflet} map widget
#' @export
#'
#' @import leaflet
#'
#' @examples
#' \dontrun{
#' generateFloodRiskMap(long = -0.08189974, lat = 51.31244, region_col_name = "NAME", risk_col_name = "PROJ_4BAND",
#'                      shapefile = shapefile, shapefile_divided = shapefile_divided, shapefile_bounds = shapefile_bounds)
#'}
generateFloodRiskMap <- function(long, lat, region_col_name, risk_col_name, shapefile,
                                 shapefile_divided, shapefile_bounds){
  content <- paste("Long: ",
                   long,
                   "<br/>Lat: ",
                   lat,
                   "<br/>Flood risk: ",
                   getRiskLevel(long = -0.08189974,
                                lat = 51.31244,
                                region_col_name = "NAME",
                                risk_col_name = "PROB_4BAND",
                                shapefile_divided = shapefile_divided,
                                shapefile_bounds = shapefile_bounds))

  map <- leaflet()  %>%
    addTiles() %>%
    setView(lng = long, lat = lat, zoom = 14) %>%
    addPolygons(data = shapefile, weight = 5, color = '#9C7AFF') %>%
    addPopups(long, lat, content,
              options = popupOptions(closeButton = FALSE)) %>%
    setMaxBounds(shapefile_bounds@bbox[1], shapefile_bounds@bbox[2],
                 shapefile_bounds@bbox[3], shapefile_bounds@bbox[4])

  return(map)
}