Commit 9d660fad authored by Conor Anderson's avatar Conor Anderson

Initial implementation of data cropping

Squashed commit of the following:

commit fe81678aa7a6cccc82381b58c1dac9f3b4d8bdea
Author: Conor Anderson <[email protected]>
Date:   Fri May 15 15:03:30 2020 -0400

    Initial complete implementation of crop

commit 062210c8
Author: Conor Anderson <[email protected]>
Date:   Fri May 15 11:28:31 2020 -0400

    Plumber changes to enable cropping

commit 9db63d19
Merge: df4f5e45 01dff937
Author: Conor Anderson <[email protected]>
Date:   Fri May 15 11:26:24 2020 -0400

    Merge branch 'crop' of gitlab.com:ConorIA/conjuntool into crop

commit df4f5e45
Merge: d20e129b 24bd33a3
Author: Conor Anderson <[email protected]>
Date:   Fri May 15 11:24:44 2020 -0400

    Merge branch 'devel' into crop

commit 01dff937
Merge: 55c65ee3 caae4217
Author: Conor Anderson <[email protected]>
Date:   Tue Jan 22 10:21:13 2019 -0500

    Merge remote-tracking branch 'origin/devel' into crop

commit d20e129b
Merge: 87f42f34 59c8e3e8
Author: Conor Anderson <[email protected]>
Date:   Fri Oct 12 19:43:19 2018 -0400

    Merge branch 'devel' into crop

commit 87f42f34
Merge: dde1388e e9b8051a
Author: Conor Anderson <[email protected]>
Date:   Fri Oct 12 19:38:48 2018 -0400

    Merge branch devel

    Merge branch 'devel' into crop

    # Conflicts:
    #	server.R
    #	ui.R

commit dde1388e
Author: Conor Anderson <[email protected]>
Date:   Thu Oct 11 18:36:41 2018 -0400

    Crop check-in

commit 55c65ee3
Author: Conor Anderson <[email protected]>
Date:   Thu Sep 27 23:22:36 2018 -0400

    Start crop improvements
parent 24bd33a3
......@@ -6,7 +6,8 @@ COPY api/plumber.R /api/plumber.R
RUN echo 'deb http://deb.debian.org/debian bullseye main' > /etc/apt/sources.list
RUN apt-get update &&\
apt-get install -y --no-install-recommends curl libnetcdf-dev &&\
apt-get -y dist-upgrade &&\
apt-get install -y --no-install-recommends curl libgdal-dev libnetcdf-dev &&\
apt-get clean && rm -rf /tmp/* /var/lib/apt/lists/*
RUN cd /api &&\
......
library(base64enc)
library(dplyr)
library(jsonlite)
library(lubridate)
......@@ -6,6 +7,7 @@ library(ncdf4.helpers)
library(plumber)
library(purrr)
library(raster)
library(rgdal)
library(storr)
library(stringr)
library(tibble)
......@@ -44,7 +46,7 @@ get_choices <- function(cmip, var, lim) {
} else {
names(choices) <- c("Variable", "Period", "Model", "Scenario", "Ensemble", "Grid", "Range", "Start", "End")
}
choices <- choices %>%
mutate(Model = gsub("([0-9])-([0-9])", paste0("\\1", ".", "\\2"), gsub("([0-9])-([0-9])", paste0("\\1", ".", "\\2"), Model)),
Scenario = sub("rcp([0-9])([0-9])", paste0("RCP", "\\1", ".", "\\2"), Scenario),
......@@ -55,11 +57,11 @@ get_choices <- function(cmip, var, lim) {
#* @filter checkAuth
function(req, res){
auth = sub("Basic ", "", req$HTTP_AUTHORIZATION)
if (!(rawToChar(jsonlite::base64_dec(auth)) %in% paste(users, passwords, sep = ":"))) {
if(req$PATH_INFO == "/download" || (rawToChar(jsonlite::base64_dec(auth)) %in% paste(users, passwords, sep = ":"))) {
plumber::forward()
} else {
res$status <- 401 # Unauthorized
return(list(error="Authentication required"))
} else {
plumber::forward()
}
}
......@@ -124,13 +126,13 @@ function(cmip = "CMIP5", var, lim) {
#* @post /timeseries
function(key) {
time_tmp <- NULL
if (st_point$exists(key)) {
message("Hit the cache for ", key)
time_tmp <- st_point$get(key)
## Invalidate old cache entries.
if (is.null(attr(time_tmp, "cache_ver")) || attr(time_tmp, "cache_ver") != cache_ver) {
message("Cached version is invalid. Deleting.")
......@@ -138,11 +140,11 @@ function(key) {
time_tmp <- NULL
}
}
if (is.null(time_tmp)) {
message("Missed the cache for ", key)
filename <- str_extract(key, ".*\\.nc4?")
# Get the time
components <- get.split.filename.cmip(filename)
cmip <- ifelse(has_name(components, "grid"), "CMIP6", "CMIP5")
......@@ -161,17 +163,17 @@ function(key) {
} else {
nc_time <- as.yearmon(format(nc_time, format = "%Y-%m-%d hh:mm:ss"))
}
cells <- unlist(strsplit(sub(paste0(filename, "_"), "", key), "_"))
lon_cell <- cells[2]
lat_cell <- cells[1]
# Now load only that grid data
nc_var <- nc.get.var.subset.by.axes(nc_nc, components[['var']], axis.indices = list(X = lon_cell, Y = lat_cell))
# Close the nc connection
nc_close(nc_nc); rm(nc_nc)
# Note, does this still apply?
if (dim(nc_var)[1] > 1 || dim(nc_var)[2] > 1) {
warning("We got two latitud or longitude cells. We'll average across them.")
......@@ -179,7 +181,7 @@ function(key) {
} else {
nc_var <- nc_var[1, 1, ]
}
time_tmp <- tibble(Var = components[['var']],
Model = gsub("([0-9])-([0-9])", paste0("\\1", ".", "\\2"), gsub("([0-9])-([0-9])", paste0("\\1", ".", "\\2"), components[['model']])),
Scenario = sub("rcp([0-9])([0-9])", paste0("RCP", "\\1", ".", "\\2"), components[['emissions']]),
......@@ -188,7 +190,7 @@ function(key) {
Year = format(as.yearmon(Time), format = "%Y"),
Month = format(as.yearmon(Time), format = "%m"),
Value = nc_var)
rm(nc_var, nc_time); gc()
attr(time_tmp, "cache_ver") <- cache_ver
st_point$set(key, time_tmp)
......@@ -204,57 +206,57 @@ function(key) {
#* @post /mapseries
function(key) {
if (st_avg$exists(key)) {
if (debug_flag) message("Hit the cache.")
map_data <- st_avg$get(key)
} else {
if (debug_flag) message("Missed the cache.")
filename <- str_extract(key, ".*\\.nc4?")
unlisted <- unlist(str_split(key, "_"))
components <- get.split.filename.cmip(filename)
cmip <- ifelse(has_name(components, "grid"), "CMIP6", "CMIP5")
nc_nc <- nc_open(file.path(file_dir, cmip, components['var'], "verified", filename), readunlim = FALSE)
nc_time <- try(nc.get.time.series(nc_nc, v = components['var'], time.dim.name = "time"))
if (inherits(nc_time, "try-error")) {
warning(paste(fs[f], "failed"))
break
}
nc_time <- as.yearmon(format(nc_time, format = "%Y-%m-%d hh:mm:ss"))
# Get the time that we are interested in
index_start <- min(which(format(nc_time, format = "%Y") == unlisted[length(unlisted)-1]))
index_end <- max(which(format(nc_time, format = "%Y") == unlisted[length(unlisted)]))
nc_lat <- ncvar_get(nc_nc, "lat")
nc_lon <- ncvar_get(nc_nc, "lon")
lon_bnds <- try(ncvar_get(nc_nc, "lon_bnds"), silent = TRUE)
if (inherits(lon_bnds, "try-error")) {
lon_bnds <- ncvar_get(nc_nc, "lon_bounds")
}
lat_bnds <- try(ncvar_get(nc_nc, "lat_bnds"), silent = TRUE)
if (inherits(lat_bnds, "try-error")) {
lat_bnds <- ncvar_get(nc_nc, "lat_bounds")
}
nc_var <- nc.get.var.subset.by.axes(nc_nc, components['var'], axis.indices = list(T = index_start:index_end))
# Close the nc connection
nc_close(nc_nc)
rm(nc_nc)
map_data_mat <- apply(nc_var, c(1,2), mean)
map_data <- raster(t(map_data_mat),
xmn = mean(lon_bnds[,1]), xmx = mean(lon_bnds[,ncol(lon_bnds)]),
ymn = mean(lat_bnds[,1]), ymx = mean(lat_bnds[,ncol(lat_bnds)]),
crs="+proj=longlat +datum=WGS84")
# This step fails if the model outputs the corner instead of the centre point.
map_data <- try(rotate(map_data))
if(inherits(map_data, "try-error")) {
......@@ -263,12 +265,12 @@ function(key) {
map_data <- raster(t(map_data_mat), xmn = (nc_lon[2] - nc_lon[1])/2, xmx = max(nc_lon) + (nc_lon[2] - nc_lon[1])/2, ymn = min(nc_lat), ymx = max(nc_lat), crs="+proj=longlat +datum=WGS84")
map_data <- rotate(map_data)
}
map_data <- flip(map_data, 'y')
st_avg$set(key, map_data)
if (debug_flag) message("Cached data.")
}
serialize(map_data, NULL)
}
......@@ -279,7 +281,7 @@ function(key) {
#* @post /cacheget
function(key){
serialize(
(if (st_meta$exists(key))st_meta$get(key) else "NULL"),
(if (st_meta$exists(key)) st_meta$get(key) else "NULL"),
NULL)
}
......@@ -292,3 +294,80 @@ function(req, key){
st_meta$set(key = key, value = as_tibble(fromJSON(req$postBody)))
}
#* Crop netcdf file and return the filename for download
#* @param key The model key to process
#* @post /crop
function(key, req, res){
message("The model is ", key)
shp <- unserialize(base64decode(req$postBody))
filename <- str_extract(key, ".*\\.nc4?")
unlisted <- unlist(str_split(key, "_"))
components <- get.split.filename.cmip(filename)
cmip <- ifelse(has_name(components, "grid"), "CMIP6", "CMIP5")
nc_nc <- nc_open(file.path(file_dir, cmip, components[['var']], "verified", filename), readunlim = FALSE)
var_details <- nc_nc$var[[components['var']]]
nc_close(nc_nc)
chbr = brick(file.path(file_dir, cmip, components['var'], "verified", filename))
if (attr(extent(chbr), "xmax") > 180) chbr = rotate(chbr)
if (is(shp, "SpatialPolygonsDataFrame")) {
message("Working with a shapefile")
if (!identical(crs(shp), crs(chbr))) shp = spTransform(shp, crs(chbr))
chbr.crop <- crop(chbr, shp)
rm(chbr); invisible(gc())
chbr.out = mask(chbr.crop, shp)
rm(chbr.crop); invisible(gc())
} else if (is(shp, "Extent")) {
message("Working with an Extent")
chbr.out <- crop(chbr, shp)
rm(chbr); invisible(gc())
} else {
stop("An unknown error occurred")
}
ddir <- file.path(file_dir, "dist")
if(!dir.exists(ddir)) dir.create(ddir, recursive = TRUE)
fileout <- file.path(ddir, paste("cropped",
format(Sys.time(), format = "%Y%m%d_%H%M%S"),
filename, sep = "_"))
writeRaster(chbr.out, fileout, overwrite=TRUE, format="CDF",
varname=var_details$name,
varunit=var_details$units,
longname=var_details$longname,
xname=var_details$dim[[1]]$name,
yname=var_details$dim[[2]]$name,
zname=var_details$dim[[3]]$name)
message("File written to ", fileout)
res$status <- 200
res$body <- jsonlite::toJSON(auto_unbox = TRUE, list(
status = 200,
message = "Data was processed successfully. The filename is included in this body.",
filename = fileout
))
return(res)
}
# If the URL gets called the browser will automatically download the file.
#* @serializer contentType list(type="application/octet-stream")
#* @param filepath The file to return for download
#* @get /download
function(filepath, res) {
filepath <- URLdecode(filepath)
if (!file.exists(filepath)) {
res$status <- 404
res$body <- jsonlite::toJSON(auto_unbox = TRUE, list(
status = 404,
message = "File does not exist")
)
return(res)
}
readBin(filepath, "raw", n = file.info(filepath)$size)
}
crop_req <- function(key, shp) {
if (debug_flag) message("Asking plumber to crop ", key)
r <- POST(URLencode(paste0(plumber_address, "/crop?key=", key)),
config = list(add_headers(accept = "application/octet-stream")),
authenticate(plumber_user, plumber_password),
body = base64encode(serialize(shp, NULL)), encode = "json")
if (r$status_code != 200) {
stop("There was an error")
}
res <- parse_json(rawToChar(content(r)))
URLencode(res$filename, reserved = TRUE)
}
This diff is collapsed.
On this page you can crop the GCM data from a single model. You can choose to upload an ESRI shapefile (with the necessary auxilary files) in a single ZIP archive, or you can manually choose a single point, or two corners of a rectangle object.
Click on the map to choose coordinates. You can verify or modify the coordinates in the "Manual Coordinates" panel.
Some notes:
You should _not_ provide both a Shapefile _and_ manual input. However, if in doubt, keep the following guidelines in mind:
- Uploading a Shapefile will reset the manual input to the default (invalid) coordinates. The Shapefile will be used.
- If you set coordiantes manually _after_ uploading a Shapefile, we will remove the Shapefile and prefer the manual input.
- If you change your mind again, click the map once more. When the box disappears, the Shapefile will be used.
- If you want to use a different Shapefile, uploading it will replace any previously uploaded Shapefiles.
**This functionality remains a work in progress and is only lightly tested. In particular, the only models that are available are those that are packaged as a single _.nc_ file from 2006 to 2100. Models that are packaged by decade, for example, are not available, nor are the historical runs. Sometime the extent that is passed up to raster doesn't work. If that happens, try again, or try using a shapefile.**
......@@ -76,6 +76,13 @@ singleModelSelect <- function(input, output, session, choices, coords, projectio
choices <- choices()
coords <- coords()
projection <- projection()
if (debug_flag) {
message("get_metadata\n\n",
"coords: ", coords,
"\nmodel: ", input$model_in,
"\nscenario: ", input$scen_in,
"\nensemble: ", input$run_in)
}
get_metadata(choices, coords, baseline = NULL, projection = projection,
model = input$model_in, scenario = input$scen_in, ensemble = input$run_in)
}))
......
## -- Load required packages -- ##
library("shiny")
library("base64enc")
library("broom")
library("crul")
library("doParallel")
library("dplyr")
library("DT")
library("gdalUtils")
library("ggplot2")
library("googleway")
library("htmlwidgets")
library("httr")
library("jsonlite")
library("leaflet")
library("lubridate")
library("mapview")
......@@ -16,12 +20,12 @@ library("purrr")
library("raster")
library("RColorBrewer")
library("readr")
library("rgdal")
library("stringr")
library("tibble")
library("tidyr")
library("rgdal")
library("zoo")
library("crul")
validate <- shiny::validate
## -- Import modules and borrowed functions -- ##
lapply(c(dir("modules", full.names = TRUE), dir("R", full.names = TRUE),
......@@ -455,4 +459,144 @@ shinyServer(function(input, output, session) {
}
)
### Resample Data (WIP)
crop_model_in <- callModule(singleModelSelect, "crop_model_in",
choices = reactive(get_choices(input$var_filter_crop, lim = TRUE)),
coords = reactiveVal(NULL),
projection = reactiveVal(NULL))
output$resample_map <- renderLeaflet({
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
setView(lng = 0, lat = 0, zoom = 3)
})
unzip_shp <- reactive({
print(sprintf("The Shapefile is %s", input$shapefile$datapath))
zip_dir <- tempdir()
print(sprintf("Unzipping to %s", zip_dir))
unzip(input$shapefile$datapath, exdir = zip_dir)
shp <- readOGR(zip_dir)
ext <- extent(shp)
list(shp = shp, ext = ext)
})
map_shapefile <- function() {
# Clear the manual inputs
updateNumericInput(session, "Lat1", value = -999)
updateNumericInput(session, "Lat2", value = -999)
updateNumericInput(session, "Lon1", value = -999)
updateNumericInput(session, "Lon2", value = -999)
shp_data <- unzip_shp()
shp <- shp_data$shp
ext <- shp_data$ext
proxy <- leafletProxy("resample_map")
proxy %>% clearMarkers() %>% clearShapes()
proxy %>% addPolygons(data = shp) %>%
flyToBounds(attr(ext, "xmin"), attr(ext, "ymin"),
attr(ext, "xmax"), attr(ext, "ymax"))
}
# Show Shapefile on map
observeEvent(input$shapefile, {
map_shapefile()
})
# Show popup on click
observeEvent(input$resample_map_click, {
click <- input$resample_map_click
proxy <- leafletProxy("resample_map")
if (input$Lat1 == -999 && input$Lon1 == -999) {
updateNumericInput(session, "Lat1", "Lat1", value = click$lat); updateNumericInput(session, "Lon1", "Lon1", value = click$lng)
text <- paste("Lattitude ", click$lat, "Longtitude ", click$lng)
proxy %>% clearMarkers() %>% clearShapes()
proxy %>% addMarkers(click$lng, click$lat, text)
}
else if(input$Lat1 != -999 && input$Lon1 != -999 && (input$Lat2 == -999 || input$Lon2 == -999)) {
updateNumericInput(session, "Lat2", value = click$lat)
updateNumericInput(session, "Lon2", value = click$lng)
ext <- list(xmin = min(click$lng, input$Lon1),
xmax = max(click$lng, input$Lon1),
ymin = min(click$lat, input$Lat1),
ymax = max(click$lat, input$Lat1))
text <- paste("Lattitude ", click$lat, "Longtitude ", click$lng)
proxy %>% addMarkers(click$lng, click$lat, text) %>%
addRectangles(click$lng, click$lat, input$Lon1, input$Lat1) %>%
flyToBounds(ext$xmin, ext$ymin, ext$xmax, ext$ymax)
}
else {
updateNumericInput(session, "Lat1", value = -999)
updateNumericInput(session, "Lat2", value = -999)
updateNumericInput(session, "Lon1", value = -999)
updateNumericInput(session, "Lon2", value = -999)
proxy %>% clearMarkers() %>% clearShapes()
if (!is.null(input$shapefile)) {
map_shapefile()
}
}
})
output$latlon_checkbox <- renderUI({
model <- try(req(crop_model_in()))
shp <- try(req(input$shapefile))
content <- '<b>To do:</b>'
if (all(c(input$Lat1, input$Lat2, input$Lon1, input$Lon2) %in% -999) &&
is.null(input$shapefile)) {
content <- c(content, '<i class="far fa-square"></i> Select at least one point or upload a .shp')
} else {
content <- c(content, '<i class="far fa-check-square"></i> Select at least one point or upload a .shp')
}
if (inherits(model, "try-error")) {
content <- c(content, '<i class="far fa-square"></i> Choose a mod/var/scen to process.')
} else {
content <- c(content, '<i class="far fa-check-square"></i> Choose a mod/var/scen to process.</i>')
}
HTML(paste(content, collapse = "<br>"))
})
output$crop_action <- renderUI({
if ((!all(c(input$Lat1, input$Lat2, input$Lon1, input$Lon2) %in% -999) ||
!is.null(input$shapefile)) &
!is.null(crop_model_in())) {
actionButton("crop_go", "Process!")
} else {
NULL
}
})
get_crop_link <- eventReactive(input$crop_go, {
message("Button pressed")
meta <- crop_model_in()
key = unlist(meta$Files)
if (!is.null(input$shapefile)) {
shp_data <- unzip_shp()
shp <- shp_data$shp
} else {
xmin <- isolate(min(input$Lon1, input$Lon2))
xmax <- isolate(max(input$Lon1, input$Lon2))
ymin <- isolate(min(input$Lat1, input$Lat2))
ymax <- isolate(max(input$Lat1, input$Lat2))
shp <- extent(c(xmin, xmax, ymin, ymax))
}
saveRDS(list(key, shp), "debug.rds")
crop_req(key, shp)
})
output$cropped_data_link = renderUI({
link <- get_crop_link()
link <- paste0(plumber_address, "download?filepath=", link)
tags$a(href = link, target = "blank", "Click here for your data")
})
output$source <- renderUI({
desc <- system("git log devel -n 20 --format='%cd %s' --date=short", intern = TRUE)
commits <- system("git log devel -n 20 --format='%h'", intern = TRUE)
changes <- paste0(desc, " (",
sprintf('<a href="https://gitlab.com/ConorIA/conjuntool/commit/%s" target="_blank">%s</a>', commits, commits),
") <br>")
HTML(changes)
})
})
......@@ -21,11 +21,17 @@ dashboardPage(
menuItem("GCM Plots", tabName = "plot", icon = icon("line-chart")),
menuItem("GCM Anomalies", tabName = "anomalies", icon = icon("thermometer-three-quarters")),
menuItem("Overlay Map", tabName = "map", icon = icon("globe")),
menuItem("Crop Data", tabName = "crop", icon = icon("crop")),
menuItem("Source Code", tabName = "source", icon = icon("git"))
)
),
dashboardBody(
dashboardBody( tags$head(
tags$style(HTML("
.shiny-output-error-validation {
color: gray;
}
"))
),
tabItems(
## README
tabItem(tabName = "README",
......@@ -58,7 +64,6 @@ dashboardPage(
textInput("plot_city_in", "Location", value = "UTSC, Scarborough, ON", width = NULL, placeholder = NULL),
leafletOutput("plot_map", height = 200)
),
box(
title = "Output",
width = 9,
......@@ -153,6 +158,64 @@ dashboardPage(
), # End tab 3
# Tab 4
tabItem(tabName = "crop",
fluidRow(
box(width = 10,
includeMarkdown("include/crop.md")
),
box(width = 2,
uiOutput("latlon_checkbox"),
uiOutput("crop_action"),
conditionalPanel("input.crop_go", htmlOutput("cropped_data_link"))
)
),
fluidRow(
box(
width = 3,
selectInput("var_filter_crop", "Select Variables", c("tas", "tasmax", "tasmin", "pr")),
singleModelSelectInput("crop_model_in", "Select Model")
#sliderInput("crop_year_in", label = "Years to Include", min = 2006,
# max = 2100, value = c(2011, 2100), step = 1,
# round = FALSE, sep = "", ticks = FALSE)
),
box(
width = 9,
leafletOutput("resample_map")
)
),
fluidRow(
box(
width = 12,
title = "Upload a Shapefile",
collapsible = TRUE,
collapsed = FALSE,
fileInput("shapefile", "Choose a zip file containing .shp and other files",
multiple = FALSE,
accept = ".zip")
)
),
fluidRow(
box(
width = 12,
title = "Manual Coordinates",
collapsible = TRUE,
collapsed = FALSE,
column(
width = 6,
numericInput("Lat1", "Latitide Bound 1", value = -999, width = -999),
numericInput("Lon1", "Longitude Bound 1", value = -999, width = -999)
),
column(
width = 6,
numericInput("Lat2", "Latitude Bount 2", value = -999, width = -999),
numericInput("Lon2", "Longitude Bound 2", value = -999, width = -999)
)
)
)
), # End tab 4
# Tab 5
tabItem(tabName = "source",
fluidRow(
box(width = 12,
......@@ -168,7 +231,7 @@ dashboardPage(
)
)
)
) # End tab 4
) # End tab 5
) # End tabItems
) # End dashboardBody
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment