In JBrenn/SpatialInterpol: Spatial Interpolation of Soil Variables | Province of Bolzano

---

This R Markdown document is made interactive using Shiny. To learn more, see Interactive Documents.

---

if(!require("SpatialInterpol"))
{
  if(!require("devtools"))
  {
    install.packages(devtools)
    require("devtools")
  }
  install_github("SpatialInterpol", "JBrenn")
  require("SpatialInterpol")
}

if(!require("leaflet"))
{
  install.packages(leaflet)
  require("leaflet")
}

if(!require("raster"))
{
  install.packages(raster)
  require("raster")
}

if(!require("rgdal"))
{
  install.packages(rgdal)
  require("rgdal")
}

wpath <- "/home/jbre/Schreibtisch/krig/"
coordsys = "+proj=utm +zone=32 ellps=WGS84"

masterfile <- dir(file.path(wpath,"master"), full.names = T)
data <- read.csv(masterfile, header = T)

Interactive 1

inputPanel(

  numericInput(inputId = "long", label = "area of interest | longitude", value = 11.14),
  numericInput(inputId = "lat", label = "area of interest | latitude", value = 46.47),

  numericInput(inputId = "intA", label = "area of interest | A", value = 11000, min = 0, max = 100000, step = 1), 

  selectInput(inputId = "variable", label = "discover variable", choices = dimnames(data)[[2]][c(7,8,10:17)], selected = "Humus____"),

  radioButtons(inputId = "mask", label = "use raster mask", choices = c("YES","NO"), selected = "YES"),

  actionButton("go1", "Apply changes")

)

---

Interactive Maps

p1 <- eventReactive(input$go1, {

  # lat / long
  xy <- project(xy = cbind(data$x_Coord, data$y_Coord), proj = coordsys, inv = T)
  long <- xy[,1]; lat <- xy[,2]

    xy_sp <- SpatialPointsDataFrame(coords =  cbind(data$x_Coord, data$y_Coord), proj4string = CRS(coordsys), data = data[,c("FID",input$variable)])

  xy <- project(xy = cbind(input$long, input$lat), proj = coordsys, inv = F)
  mat <- cbind(long = c(xy[1,1],xy[1,1]+input$intA,xy[1,1]+input$intA,xy[1,1]), 
                lat = c(xy[1,2],xy[1,2],xy[1,2]+input$intA, xy[1,2]+input$intA))

  P1  <- Polygon(mat)
  Ps1 <- Polygons(list(P1), ID = "a")
  SPs <- SpatialPolygons(list(Ps1), proj4string = CRS(coordsys))

  pointsIN <- over(xy_sp,SPs)
  pointsIN <- which(!is.na(pointsIN))
  pt.poly <- xy_sp[pointsIN,]  
  pointsIN_latlong <- project(xy = pt.poly@coords, proj = coordsys, inv = T)

  latlong <- as.data.frame(project(mat,  proj = coordsys, inv = T))

  m <- leaflet() %>%  
  addTiles(group = "OSM (default)") %>%
  addCircleMarkers(lng=long, lat = lat, radius = 3, group = "observations", opacity = .5, color = "orange",
                   popup = paste(data$ComuneCata, data$TipoColtur, "| ", input$variable, data[,input$variable]),
             clusterOptions = markerClusterOptions()) %>%
  addPolygons(data = latlong, lng = ~long, lat = ~lat, fill = F, weight = 2, color = "#FFFFCC", group = "area of interest")  %>% 
  addCircleMarkers(lng=input$long, lat=input$lat, radius=2.5, fill = T, color = "red", popup = "input point") %>%
  addCircleMarkers(lng=pointsIN_latlong[,1], lat=pointsIN_latlong[,2], radius=1.5, fill = F, color = "grey", group = "points of interest")

  m

})

renderLeaflet({
  p1()
})

---

p2 <- eventReactive(input$go2, {
  # lat / long
  xy <- project(xy = cbind(data$x_Coord, data$y_Coord), proj = coordsys, inv = T)
  long <- xy[,1]; lat <- xy[,2]

   xy_sp <- SpatialPointsDataFrame(coords =  cbind(data$x_Coord, data$y_Coord), proj4string = CRS(coordsys), data = data[,c("FID",input$variable)])

  xy <- project(xy = cbind(input$long, input$lat), proj = coordsys, inv = F)
  mat <- cbind(long = c(xy[1,1],xy[1,1]+input$intA,xy[1,1]+input$intA,xy[1,1]), 
                lat = c(xy[1,2],xy[1,2],xy[1,2]+input$intA, xy[1,2]+input$intA))
  P1  <- Polygon(mat)
  Ps1 <- Polygons(list(P1), ID = "a")
  SPs <- SpatialPolygons(list(Ps1), proj4string = CRS(coordsys))

  pointsIN <- over(xy_sp,SPs)
  pointsIN <- which(!is.na(pointsIN))
  pt.poly <- xy_sp[pointsIN,]  
  pointsIN_latlong <- project(xy = pt.poly@coords, proj = coordsys, inv = T)

  # write out in master files
  df_out <- data.frame(pt.poly@data, pt.poly@coords)
  names(df_out) <- c(names(pt.poly@data), "x_Coord", "y_Coord")

  sim_path <- paste("sim", input$long, input$lat, input$intA, sep="_")
  area     <- "sim"
  dir.create(file.path(wpath, sim_path, "master"), recursive = T)
  write.csv(x = df_out, file = file.path(wpath, sim_path, "master", paste("masterfile_",area,".txt",sep="")), row.names = F, quote = F)

  # use mask?
  if (input$mask == "YES" & !file.exists(file.path(wpath, sim_path, "mask", "mask.tif")))
  {
    mask_org_file <- grep(x = dir(file.path(wpath, "mask")), pattern = ".tif")[1]
    mask_org_file <- dir(file.path(wpath, "mask"), full.names = T)[mask_org_file]
    mask <- raster(mask_org_file)
    crs(mask) <- coordsys
    mask_masked <- mask(mask, SPs)
    dir.create(file.path(wpath, sim_path, "mask"), recursive = T)
    writeRaster(x = mask_masked, filename = file.path(wpath, sim_path, "mask", "mask.tif"))
  }

  if (input$mask == "YES") {
    rastermask <- "mask/mask.tif"            
  } else {
    rastermask = NA
  }

if (input$doglobal == "YES" & !any(grep("global", dir(file.path(wpath, sim_path, input$variable, "maps")))))
{
  global <- OrdKrig(wpath = file.path(wpath, sim_path), datafolder = "master", rastermask = rastermask, local = FALSE, variable = input$variable, inverseDistWeigths = FALSE, npix = input$npix, cutoff = input$cutoff, anis_deg = input$anis_deg, anis_ax = input$anis_ax, psill = input$psill, nugget = input$nugget, nmax = input$nmax, nmin = input$nmin, omax = input$omax, model = input$model, validation = FALSE)

r_global <-  global[[area]]$map_pred
}

if (any(grep("global", dir(file.path(wpath, sim_path, input$variable, "maps"))))) {
  which2load <- grep("global_predict", dir(file.path(wpath, sim_path, input$variable, "maps")))
  files <- dir(file.path(wpath, sim_path, input$variable, "maps"), full.names = T)
  r_global <- raster(files[which2load])
}  

local <- OrdKrig(wpath = file.path(wpath, sim_path), datafolder = "master", rastermask = rastermask, local = TRUE, variable = input$variable, inverseDistWeigths = FALSE, npix = input$npix, cutoff = input$cutoff, anis_deg = input$anis_deg, anis_ax = input$anis_ax, psill = input$psill, nugget = input$nugget, nmax = input$nmax, nmin = input$nmin, omax = input$omax, model = input$model, validation = FALSE)

r_local <-  local[[area]]$map_pred

which2load <- grep("local_predict", dir(file.path(wpath, sim_path, input$variable, "maps")))
files2load <-  dir(file.path(wpath, sim_path, input$variable, "maps"), full.names = T)[which2load]
r_local_list <- lapply(files2load, raster)
names(r_local_list) <- dir(file.path(wpath, sim_path, input$variable, "maps"))[which2load]
r_local_val  <- lapply(r_local_list, values)

latlong <- as.data.frame(project(mat,  proj = coordsys, inv = T))

if (input$doglobal == "YES") {
  all_val <-  c(values(r_global),unlist(r_local_val))
} else {
  all_val <-  unlist(r_local_val)
}

pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), all_val, na.color = "transparent")

m <- leaflet() %>%  
  addTiles(group = "OSM (default)") %>%
   addLegend(pal = pal, values = all_val, title = input$variable, position = "topleft")  %>% 
  #addCircleMarkers(lng=long, lat = lat, radius = 3, group = "observations", opacity = .5, color = "orange",
  #                 popup = paste(data$ComuneCata, data$TipoColtur, "| ", input$variable, data[,input$variable]),
  #           clusterOptions = markerClusterOptions()) %>%
  addPolygons(data = latlong, lng = ~long, lat = ~lat, fill = F, weight = 2, color = "#FFFFCC", group = "area of interest")  %>% 
  addCircleMarkers(lng=input$long, lat=input$lat, radius=2.5, fill = T, color = "red", popup = "input point") %>%
  addCircleMarkers(lng=pointsIN_latlong[,1], lat=pointsIN_latlong[,2], radius=1.5, fill = F, color = "grey", group = "points of interest", opacity = .2)

overlayGroups <- c("area of interest", "points of interest")

for(i in names(r_local_list))
{
  groupname <- strsplit(split = "_krige_", x = i)[[1]][2]
  groupname <- substr(groupname, 1, nchar(groupname)-4)

  m <- m  %>% 
    addRasterImage(x = r_local_list[[i]],  colors = pal, opacity = 0.7, group = groupname)
    overlayGroups <- c(overlayGroups, groupname)
}

if (input$doglobal == "YES") {
  m <- m  %>% 
  addRasterImage(x = r_global, colors = pal, opacity = 0.8, group = "global kriging")
  overlayGroups <- c(overlayGroups, "global kriging")
}

m <- m %>%
      addLayersControl(
    overlayGroups = overlayGroups, options = layersControlOptions(collapsed = FALSE))
m  
})

renderLeaflet({
  p2()
})

Interactive 2

inputPanel(

  numericInput(inputId = "cutoff", label = "cutoff", value = 400, min = 0, max = 10000, step = 1),
  numericInput(inputId = "anis_deg", label = "anis_deg", value = 0, min = 0, max = 360, step = 1), 
  numericInput(inputId = "anis_ax", label = "anis_ax", value = 0.5, min = 0, max = 1, step = 0.01), 
  numericInput(inputId = "psill", label = "psill", value = 1),
  numericInput(inputId = "nugget", label = "nugget", value = 1),
  numericInput(inputId = "nmax", label = "nmax", value = 100),
  numericInput(inputId = "nmin", label = "nmin", value = 25),
  numericInput(inputId = "omax", label = "omax", value = 5),

  selectInput(inputId = "model", label = "model", choices = c("Exp",  "Sph", "Gau", "Mat"), selected = "Sph"),
  numericInput(inputId = "npix", label = "resolution", value = 100, min = 0, max = 1000, step = 1),

  radioButtons(inputId = "doglobal", label = "use global kriging?", choices = c("YES","NO"), selected = "NO"),

  actionButton("go2", "Apply changes")

)

---

JBrenn/SpatialInterpol documentation built on May 7, 2019, 7:39 a.m.