otherR_functions/visualizations/map_rd_sa.R

library(leaflet)
library(mapview)
map_rd <- function(aws_name, sensor_name, addition = "", buffers, vegetation_height_raster, vegetation_radius, AHN3 = FALSE, aws.df = AWS.df, class = "", zoom = 8){
#   tag.map.title <- tags$style(HTML("
#   .leaflet-control.map-title { 
#     transform: translate(-50%,20%);
#     position: fixed !important;
#     left: 50%;
#     text-align: center;
#     padding-left: 10px; 
#     padding-right: 10px; 
#     background: rgba(255,255,255,0.75);
#     font-weight: bold;
#     font-size: 28px;
#   }
# "))
#     
  
  
  single_aws <- select_single_aws(aws_name = aws_name, sensor_name = sensor_name, aws.df = aws.df)[["aws.df"]]
  aws_name_trim <- getAWS_name_trim(aws.df = aws.df, aws_name = aws_name, addition = addition)
  
  #print(class(vegetation_height_raster))
  if(AHN3 == TRUE){
    AHN <- "AHN3"
  } else {
    AHN <- "AHN2"
  }
  
  if(missing(buffers)){
    circles <- FALSE
  } else {
    circles <- TRUE
    if(100 %in% buffers == FALSE){
      buffers <- append(100, buffers)
    }
  }
  
  if(missing(vegetation_height_raster)){
    vegetation <- FALSE
  } else {
    vegetation <- TRUE
    #vegetation_height_raster <- raster(paste0("output/", aws_name_trim, "/vegetation_height/", aws_name_trim, "_", AHN, "_", vegetation_radius, "m_height_difference.tif"))
    path <- paste0("output/", aws_name_trim, "/solar_shadow_angles/rasters/", AHN, "/Shadows/DeBilt_AHN3_sa_99.7438895423655.tif")
    print(path)
    shadow_angles_raster <- raster(path)
  }
  #vegetation <- FALSE
  #resolutions <- c(3440.640, 1720.320, 860.160, 430.080, 215.040, 107.520, 53.760, 26.880, 13.440, 6.720, 3.360, 1.680, 0.840, 0.420)
  # epsg28992 <- leafletCRS(crsClass= 'L.Proj.CRS', code = 'EPSG:28992',
  #                         proj4def = '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs',
  #                         resolutions = resolutions,
  #                         bounds = c(-285401.92, 22598.08, 595401.9199999999, 903401.9199999999)
  #                         #, origin = c(-305401.92, 22598.08))  #c(7.1389803, 48.1833305)
  #                         )                 
  overlays <- c("AWS", "AHN2", "AHN3")
  luchtfoto <- "https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wmts/2018_ortho25/EPSG:3857/{z}/{x}/{y}.jpeg"
  map <- leaflet(options = leafletOptions(
    #crs = epsg28992, 
    minZoom = 4, maxZoom = 20
    )) %>%
    setView(single_aws[1,"LON"], single_aws[1,"LAT"], zoom = zoom) %>%
    addTiles(luchtfoto, attribution = "KNMI, PDOK, AHN, 2018") %>%
    addScaleBar(position = "bottomright", options = scaleBarOptions(maxWidth = 100, metric = TRUE, imperial = FALSE, updateWhenIdle = TRUE)) %>%
    # addMeasure(
    #   position = "bottomleft",
    #   primaryLengthUnit = "meters",
    #   primaryAreaUnit = "sqmeters",
    #   activeColor = "#3D535D",
    #   completedColor = "#7D4479") %>%
    addMouseCoordinates() %>%
    addWMSTiles(
      "https://geodata.nationaalgeoregister.nl/ahn2/wms?"
      , layers = "ahn2_05m_ruw"
      , options = WMSTileOptions(format = "image/png", transparent = FALSE)
      , group = "AHN2"
    )%>%
    addWMSTiles(
      "https://geodata.nationaalgeoregister.nl/ahn3/wms?"
      , layers = "ahn3_05m_dsm"
      , options = WMSTileOptions(format = "image/png", transparent = FALSE)
      , group = "AHN3"
    ) #%>%
    #addControl(title, position = "bottomleft", className="map-title")
  
    # addWMSTiles(
    #   "https://geodata.nationaalgeoregister.nl/beta/bgt/wms?"
    #   , layers = "BGT WMS"
    #   , options = WMSTileOptions(format = "image/png", transparent = FALSE)
    #   , group = "BGT"
    # ) %>%
    #addMarkers(data = )
  
    # %>% addLegend("bottomright", pal = pal, values = ~gdp_md_est,
    #            title = "Est. GDP (2010)",
    #            labFormat = labelFormat(prefix = "$"),
    #            opacity = 1
    # )
  
  ## Add Land use (BGT)
  if(circles == TRUE){  
    overlays <- append(overlays, c("Land use (BGT)", "Circle areas"))
    my_bgt <- st_read(dsn = paste0("data/BGT/", aws_name_trim, "/BGT_", aws_name_trim, ".shp"))
    
    bgt <- st_transform(my_bgt, epsg_rd)
    my_bgt_wgs <- st_transform(bgt, "+init=epsg:4326")
    
    fillColours <- as.character(my_bgt_wgs$hexColour)
    bgt_objects <- as.character(my_bgt_wgs$object)
    bgt_objects <- unique(bgt_objects)
    bgt_colours <- unique(fillColours)
    id <- 1:nrow(my_bgt_wgs)
    bgt_counter <- 1
    for(i in id){
      map <- addPolygons(map, data = my_bgt_wgs[id[bgt_counter],], smoothFactor = 0.5,
                         opacity = 1.0, color = 'black', weight = 0.4,
                         fillColor = fillColours[id[bgt_counter]], fillOpacity = .7, 
                         group = "Land use (BGT)"
                         )
      bgt_counter <- bgt_counter + 1
    }
  }
  #add height difference raster
  if(vegetation == TRUE){
    overlays <- append(overlays, "Shadow angles")
    pal <- colorNumeric(c("#fafa00", "#000000"), values(shadow_angles_raster), na.color = "transparent")
    map <- addRasterImage(map, shadow_angles_raster, colors = pal, opacity = 1.0, group = "Shadow angles")
  }
  if(circles == TRUE){
    #add circles
    for (b in 1:length(buffers)){
      map <- addCircles(map, data = single_aws,
                        lng = ~LON, lat = ~LAT,
                        radius = buffers[b], 
                        opacity = 1.0, color = 'black', weight = 2.0,
                        fill = TRUE, fillOpacity = 0,
                        label = paste(buffers[b], "m"),
                        group = "Circle areas")
    }
  }
  
  if(circles == TRUE) {
    aws_label <- paste0(aws_name, " (land use class ", class,")")
  } else {
    aws_label <- paste0(aws_name)
  }
  map <-     addCircleMarkers(map, data = single_aws[1,], lng = single_aws[1,"LON"], lat = single_aws[1,"LAT"],
                              radius = 3,  stroke = TRUE, color = "red", weight = 1.0,
                              fillColor = "red", fillOpacity = 1.0,
                              group = "Temperature sensor"
  ) %>%
    addCircleMarkers(lng = 5.17969000, lat = 52.09883000,
                     radius = 3,  stroke = TRUE, color = "#45f215", weight = 1.0,
                     fillColor = "#45f215", fillOpacity = 1.0,
                     group = "Anemometer"
    )
  #add legend and control
  map <- addLegend(map, position = "topright", values = aws_name, color = "red", labels = aws_label,
              title = " Air temperature sensor", group = "AWS")
  map <- addLegend(map, position = "topright", values = aws_name, color = "#45f215", labels = aws_label,
                   title = "Wind mast", group = "Anemometer")
  if(circles == TRUE){
    map <- addLegend(map, position = "topright",
                    values = bgt_objects, colors = bgt_colours, labels = bgt_objects,
                    title = "Land use", group = "Land use (BGT)")
  }
  if(vegetation == TRUE){
    map <- addLegend(map, position = "topright", values = values(shadow_angles_raster), pal = pal,
                title = "Shadow angle (degrees)", group = "Shadow angles")
  }

  map <- addLayersControl(map,
    position = "topleft",
    baseGroups = c("Satellite"),
    overlayGroups = overlays,
    options = layersControlOptions(collapsed = FALSE)) %>%
  setView(single_aws[1,"LON"], single_aws[1,"LAT"], zoom = zoom)
  map <- hideGroup(map, "AHN2")
  map <- hideGroup(map, "AHN3")
  return(map)
}

shadow_angle_map <- map_rd(aws_name = "De Bilt",
       sensor_name = temperature_sensor_name,
       zoom = 18,
       #buffers = c(50, 30, 10),
       vegetation_height_raster = TRUE,
       vegetation_radius = 10,
       AHN3 = TRUE)

shadow_angle_map

mapshot(shadow_angle_map, file = "shadow_angles_map.png", url = NULL, remove_url = TRUE, remove_controls = c("zoomControl", "layersControl", "homeButton"))
veg_map <- function(){
  ras <- raster(paste0("output/DeBilt/vegetation_height/DeBilt_AHN3_10m_height_difference.tif"))
  breaks <- c(0,0.1,0.4,0.6)
  pal <- colorNumeric(c("#99d8c9", "#2ca25f"), values(ras), na.color = "transparent")

 raster::plot(ras)
 
 #ggsave(filename = "C:/Users/3691233/Dropbox/Apps/ShareLaTeX/Automated Temperature SC  MSc Thesis/Figures/Results/vegetation_height.png", plot = last_plot())
}

veg_map()
?colorNumer
Jellest/wmoSC documentation built on July 31, 2019, 12:34 p.m.