5

EDIT: An example setup can be found here: The main file is extract_cordex.R

I need to analyze climate data as part of my thesis, however I am not particularly trained in working with climate data and especially the remapping part, so I'd be glad if someone could look at what I am doing and in case that is correct, someone might find the functions useful.

I am interested in climate projections for a catchment located in north western India. My first step was to retrieve netcdf files for precipitation and the relevant variables for calculating the reference evapotranspiration (mean/min/max temperature, relative humidity, cloud cover/sunshine hours, windspeed) for the RCP45 and RCP85 experiments from the Earth System Grid Federation.

The domain is EAS44, the model is MPI-ESM, I chose downscaled data as part of CORDEX. An example filename is:

pr_EAS-44_MPI-M-MPI-ESM-LR_rcp45_r1i1p1_CLMcom-CCLM5-0-2_v1_mon_200601-201012.nc

The ncdump header is:

ncdump pr_EAS-44_MPI-M-MPI-ESM-LR_rcp45_r1i1p1_CLMcom-CCLM5-0-2_v1_mon_200601-201012.nc
netcdf pr_EAS-44_MPI-M-MPI-ESM-LR_rcp45_r1i1p1_CLMcom-CCLM5-0-2_v1_mon_200601-201012 {
dimensions:
    bnds = 2 ;
    rlon = 203 ;
    rlat = 167 ;
    time = UNLIMITED ; // (60 currently)
variables:
    char rotated_pole ;
        rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
        rotated_pole:grid_north_pole_latitude = 77.61f ;
        rotated_pole:grid_north_pole_longitude = -64.78f ;
    double rlon(rlon) ;
        rlon:axis = "X" ;
        rlon:standard_name = "grid_longitude" ;
        rlon:long_name = "longitude in rotated pole grid" ;
        rlon:units = "degrees" ;
    double lon(rlat, rlon) ;
        lon:standard_name = "longitude" ;
        lon:long_name = "longitude" ;
        lon:units = "degrees_east" ;
    double rlat(rlat) ;
        rlat:axis = "Y" ;
        rlat:standard_name = "grid_latitude" ;
        rlat:long_name = "latitude in rotated grid" ;
        rlat:units = "degrees" ;
    double lat(rlat, rlon) ;
        lat:standard_name = "latitude" ;
        lat:long_name = "latitude" ;
        lat:units = "degrees_north" ;
    double time(time) ;
        time:units = "days since 1949-12-01 00:00:00" ;
        time:standard_name = "time" ;
        time:long_name = "time" ;
        time:calendar = "proleptic_gregorian" ;
        time:bounds = "time_bnds" ;
    double time_bnds(time, bnds) ;
        time_bnds:long_name = "time bounds" ;
    float pr(time, rlat, rlon) ;
        pr:standard_name = "precipitation_flux" ;
        pr:long_name = "Precipitation" ;
        pr:units = "kg m-2 s-1" ;
        pr:cell_methods = "time: mean" ;
        pr:coordinates = "lon lat" ;
        pr:grid_mapping = "rotated_pole" ;
        pr:missing_value = 1.e+20f ;
        pr:_FillValue = 1.e+20f ;

// global attributes:
        :Conventions = "CF-1.4" ;
        :conventionsURL = "http://www.cfconventions.org" ;
        :creation_date = "2017-09-07T11:06:59" ;
        :contact = "http://coastmod.hzg.de" ;
        :experiment_id = "rcp45" ;
        :driving_model_id = "MPI-M-MPI-ESM-LR" ;
        :driving_model_ensemble_member = "r1i1p1" ;
        :driving_experiment_name = "rcp45" ;
        :frequency = "mon" ;
        :institute_id = "CLMcom" ;
        :rcm_version_id = "v1" ;
        :model_id = "CLMcom-CCLM5-0-2" ;
        :project_id = "CORDEX" ;
        :CORDEX_domain = "EAS-44" ;
        :product = "output" ;
        :experiment = "rcp45" ;
        :driving_experiment = "MPI-M-MPI-ESM-LR, rcp45, r1i1p1" ;
        :Institution = "Helmholtz-Zentrum Geesthacht" ;
        :references = "http://cordex.clm-community.eu" ;
        :institution_run_id = "r15m4r4" ;
        :institution_data_path = "/hpss/arch/gg0302/g260068/CORDEX-EA/r15m4r4" ;
        :source = "Climate Limited-area Modelling Community (CLM-Community)" ;
        :title = "CLMcom-CCLM5-0-2 model output prepared for CORDEX rcp45" ;
        :comment = "Please use the following reference for this climate 

I processed the data using the following steps and functions:

1. Remap the data to lon lat and the catchment area

The first thing I did was writing an R function that takes the lon lat shapefile of my catchment and produces a grid that is later used for interpolation using cdo:

cdo_create_lonlat_grid <- function(shape, grid_steps, grid_name = 'grid'){

  xfirst <- bbox(shape)[1,1]
  xlast <- bbox(shape)[1,2]
  yfirst <- bbox(shape)[2,1]
  ylast <- bbox(shape)[2,2]

  x_grid_stepsize <- (xlast-xfirst)/grid_steps
  y_grid_stepsize <- (ylast-yfirst)/grid_steps
  grid <- c(
    'gridtype  = lonlat',
    paste0('gridsize  = ', grid_steps^2),
    'xname     = lon',
    'xlongname = longitude',
    'xunits    = degrees_east',
    'yname     = lat',
    'ylongname = latitude',
    'yunits    = degrees_north',
    paste0('xsize     = ', grid_steps),  
    paste0('ysize     = ', grid_steps), 
    paste0('xfirst    = ', xfirst),
    paste0('xinc      = ', x_grid_stepsize),
    paste0('yfirst    = ', yfirst),
    paste0('yinc      = ', y_grid_stepsize)
  )
  write_lines(grid, grid_name)
}

That should produce an equally spaced grid that has the extents of the provided shapefile.

Then I create a list of all nc files and use cdo remapbil and the grid to bilinearly interpolate every nc file from the rotated grid to the normal lon-lat grid. I use the foreach package for parallel computation.

cdo_remap_bilinear_batch <- function(nc_filenames, outputfolder, grid_name, overwrite = FALSE){
  # register one core less than the sum of present cores for computation
  library(doMC)
  library(foreach)
  library(purrr)
  library(dplyr)
  library(stringr)

  registerDoMC(cores = detectCores() - 1)

  foreach(i = seq_along(nc_filenames)) %dopar% {
    outputfile <- str_c(outputfolder, '/', basename(nc_filenames[i]))
    command <- paste0(
      'cdo remapbil,',
      grid_name,
      ' ',
      nc_filenames[i],
      ' ',
      outputfile
    )
    if (!file.exists(outputfile)) {
      system(command)
    } else if (file.exists(outputfile) & overwrite) {
      system(command)
    }
    outputfile
  } %>% 
    map_chr(~ .)

}

2. Summarise and extract the data

The former function returns a character vector containing paths to the remapped nc files. The following function then reads one file into a brick (kind of a list of rasters), masks each rasterlayer with the shapefile and then calculates the mean of the masked layer. Additional information is extracted from the information embedded in the ǹc file and its name. The return value is a tibble in long format. For those unfamiliar with the pipe operator %>%: It takes the value left of the operator and puts it as first argument into the function at the right of the operator. This enables function chaining akin to method chaining on objects with the .. In practice I load all of the packages outside of the functions.

extract_cordex_data <- function(nc_file, shape){

  library(raster)
  library(tibble)
  library(dplyr)
  library(stringr)

  raster_brick <- nc_file %>% brick

  name_components <- nc_file %>% basename %>% str_split('_')

  climate_scenario <- name_components[[1]][4]

  parameter_scenario <- name_components[[1]][5]

  variable_name <- name_components[[1]][1]

  timestep <- name_components[[1]][8]

  # crop should be redundant

  spatial_subset <- raster_brick %>% crop(shape) %>% mask(shape) 

  # calculate summary for each timestep, ignoring NAs 

  values <- cellStats(spatial_subset, stat = 'mean', na.rm = TRUE)

  # extract dates of timesteps and convert to date format

  dates <- values %>% names %>% str_sub(2) %>% ymd

  # create ouput table for the current nc file

  tibble(
    date = dates,
    value = values,
    variable = variable_name,
    variable_long_name = raster_brick@title,
    unit = raster_brick@data@unit,
    climate_scenario = climate_scenario,
    parameter_scenario = parameter_scenario,
    timestep = timestep
  ) 

}
  • Interesting. Can you upload example data? – aldo_tapia Feb 9 '18 at 15:26
  • I would say that your approach is fine, it's similar to what I do (I did not test each script directly). However, I usually find that processing heavy data with R or Python is quite a bit slower than using CDO or NCO. This is true especially if you work on long time periods and multiple models at the same time. Moreover raster is very slow with files that are chunked, because it completely ignores chunking. Lastly, a tibble output is very inefficient in memory usage, compared to reading only what you need each time you need it. But if this works for you, well done! – AF7 Feb 10 '18 at 9:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.