ErikKusch / KrigR

An R Package for downloading, preprocessing, and statistical downscaling of the European Centre for Medium-range Weather Forecasts ReAnalysis 5 (ERA5) family provided by the European Centre for Medium‐Range Weather Forecasts (ECMWF).
MIT License
104 stars 24 forks source link

Unable to compile and write kriged product #57

Closed ndichistan closed 1 year ago

ndichistan commented 1 year ago

Hi Erik,

I am having some issues writing the kriged product. I always have this error when I run the script below script using Rscript

##Historical Baseline-1
train_HIST_25km <- stack("tas_day_CNRM-ESM2-1_ssp245_r1i1p1f2_gr_all_20152024.nc")
train_HIST_25km <- crop(train_HIST_25km,Extent_crop)
train_HIST_25km <- mask(train_HIST_25km, Extent_crop)

## Covariate Data
GMTED_DE_9km <- download_DEM(
  Train_ras = train_HIST_25km,
  Target_res = 0.008334,
  Shape = Extent_crop,
  Keep_Temporary = TRUE,
  Dir = Dir.Covariates
)
## Kriging
Output_HIST_9km <- krigR(
  Data = train_HIST_25km,
  Covariates_coarse = GMTED_DE_9km[[1]], 
  Covariates_fine = GMTED_DE_9km[[2]],  
  Keep_Temporary = FALSE,
  Cores = 10,
  Dir = Dir.Exports,  
  FileName = "CNRM-ESM2-1_ssp_245_tas_1km_1.nc", 
  nmax = 40
)
Commencing Kriging
  |======================================================================| 100%
Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for function 'rast': cannot open the connection
Calls: krigR ... hdr -> .writeHdrRaster -> file -> .handleSimpleError -> h
In addition: There were 50 or more warnings (use warnings() to see the first 50)
Execution halted

I can provide more information if necessary. Thank you.

ErikKusch commented 1 year ago

Hiya,

I'm afraid I've not seen this error before. Could you:

  1. Make sure you run the latest version of KrigR by reinstalling it
  2. Make the script self-contained so I can test it on my machine or
  3. Provide the data files necessary for this script

It looks like kriging works fine, but saving doesn't. Are you working with a network drive?

ndichistan commented 1 year ago

Hi Erik,

  1. I already have the latest version of the KrigR installed.
  2. I have made a self-contained script as suggested.
  3. I have added the data file necessary for the script.
  4. I am running the script on a server using Rscript.
#rm(list=ls())
#gc() #free up memory and report the memory usage.

library(KrigR)
library(rosm)
library(gimms)
library(cowplot)
library(viridis)
library(ggplot2)
library(tidyr)
library(ncdf4)
library(raster)
library(rgdal)
library(tidync)
library(parallel)
library(rnaturalearth)
library(rgeos)
library(sf)
library(terra)

setwd('/local1/Krig/test10yr')

####### create functions
source("FUN_Plotting.R")  #Same as from the tutorial

##create directories
Dir.Base <- getwd() # identifying the current directory
Dir.Data <- file.path(Dir.Base, "Data") # folder path for data
Dir.Shapes <- file.path(Dir.Data, "Shapes") # folder path for shapefiles
Dir.Covariates <- file.path(Dir.Base, "Covariates")
dir.create(Dir.Covariates)
Dirs <- sapply(c(Dir.Data, Dir.Shapes), function(x) if (!dir.exists(x)) dir.create(x))
Dir.StateExt <- file.path(Dir.Data, "State_Extent")
dir.create(Dir.StateExt)
Dir.StatePipe <- file.path(Dir.Data, "State_Pipe")
Dir.GMTED2010 <- file.path(Dir.Covariates, "GMTED2010")
Dir.Exports <- file.path(Dir.Base, "Exports")
dir.create(Dir.StatePipe)
dir.create(Dir.GMTED2010)
dir.create(Dir.Exports)

####Define Shape file for the project
download.file("https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces.zip",
              destfile = "highres.zip")
unzip("highres.zip")
Shape_shp <- readOGR("ne_10m_admin_1_states_provinces.shp")
#Shape_shp <- Shape_shp[Shape_shp$name_en %in% c("Texas"), ]
Extent <- extent(c(-97.5, -93.5, 28.5, 33.8)) # roughly the extent of Saxony
ggplot() +
  geom_polygon(data = Shape_shp, aes(x = long, y = lat, group = group), colour = "brown", fill = "black") +
  theme_bw() +
  labs(x = "Longitude", y = "Latitude")

####define extent
e <- extent(c(-97.5, -93.5, 28.5, 33.8)) # roughly the extent of Saxony
Extent <- as(e, 'SpatialPolygons')
crs(Extent) <- "+proj=longlat +datum=WGS84 +no_defs"
shapefile(Extent, 'Extent.shp', overwrite = T)
#bmaps.plot(bbox = Extent, type = "AerialWithLabels", quiet = TRUE, progress = "none")

Extent_crop <- crop(x = Shape_shp, 
                    y = Extent)
ggplot() +
  geom_polygon(data = Extent_crop, aes(x = long, y = lat, group = group), colour = "darkorange", fill = "lightgrey") +
  theme_bw() +
  labs(x = "Longitude", y = "Latitude")

##### Statrting Kriging ####################
### TAS

KrigStart1 <- Sys.time()
KrigStart1
##Historical Baseline-1
train_HIST_25km <- stack("/local1/Krig/test10yr/tas_day_CNRM-ESM2-1_historical_19701979.nc")
train_HIST_25km <- crop(train_HIST_25km,Extent_crop)
train_HIST_25km <- mask(train_HIST_25km, Extent_crop)

## Covariate Data
GMTED_DE_9km <- download_DEM(
  Train_ras = train_HIST_25km,
  Target_res = 0.008334,
  Shape = Extent_crop,
  Keep_Temporary = TRUE,
  Dir = Dir.Covariates
)

## Kriging
Output_HIST_9km <- krigR(
  Data = train_HIST_25km,
  Covariates_coarse = GMTED_DE_9km[[1]], 
  Covariates_fine = GMTED_DE_9km[[2]],  
  Keep_Temporary = FALSE,
  Cores = 10,
  Dir = Dir.Exports,  
  FileName = "CNRM-ESM2-1_hist_tas_1km_3.nc", 
  nmax = 40
)

##plotting Historical
Plot_Krigs(Output_HIST_9km[[1]],
           Shp = Extent_crop,
           Dates = "CNRM-ESM2-1 Historical 1km", columns = 2)

KrigStop1 <- Sys.time()
KrigTime1 <- KrigStop1 - KrigStart1
KrigTime1

tas_day_CNRM-ESM2-1_historical_19701979.nc

ErikKusch commented 1 year ago

Perfect! I have requested access to the data in you GDrive.

ErikKusch commented 1 year ago

Thanks for the access to the data and the well-put-together script. I am running the Kriging now on my end. Since this is a rather large file (many layers) you are trying to Krig, this will take some time. I have already successfully managed to krig and save the first 10 layers. I will keep you updated.

ErikKusch commented 1 year ago

Ok - so I am afraid I cannot reproduce this issue. I have now kriged the complete data set and it was saved successfully to my hard drive.

Maybe, to troubleshoot this, you might want to run the krigR() function with Cores = 1. I imagine that the error you receive may trace back to an incomplete writing of one of the temporary files saved to the Kriging_CNRM-ESM2-1_hist_tas_1km_3.nc directory in your Dir.Exports.

However, since I have the data ready, it would be a waste to keep it from you. Please download the kriged product, standard error, as well as temporary files from here. This link will be active for the next seven days.

Note that this is not a precedent for me as a kriging service (my server admins would take me out back if I offered that). :-D

ndichistan commented 1 year ago

Thanks for taking some time to look into this. May I ask if this test was done on a server or a local machine?

PS. Thanks for the data.

ErikKusch commented 1 year ago

I ran it on both - on my local machine (a mac), I kriged half the data (to give the laptop at least part of it's well-deserved weekend. The full data I am sharing with you was kriged on a windows server.

ndichistan commented 1 year ago

Thanks for all your help. I'll check the machine setup to see if it is responsible. Thanks again for the help.