1

I'm processing yearly multilayer netCDF files with daily precipitation data from CHIRPS. I have the files for the whole world, each file about 1.2gb large. I need to calculate indices from the precipitation data for each cell in the raster for a specific region. In order to do that I'm trying to crop the files to get a rectangular shape above my area of interest using the raster R package.

This is the code I'm using, exemplary for the first file.

library(ncdf4)
library(raster)
library(rgdal)

# Crop extent
crop_extent <- as(raster::extent(79, 89, 25, 31), "SpatialPolygons")
proj4string(crop_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# Set directory with original files
setwd("~/data")

# Read file
chirps81 <- stack("chirps-v2.0.1981.days_p05.nc")
chirps81crop <-crop(chirps1981, crop_extent)

# Write cropped file back to different folder
setwd("~/croppeddata")
writeRaster(chirps81crop, "chirps81crop.nc", overwrite=TRUE)

For some reason however while writing the file the layers lose their name. In the original files and after cropping the names have layer names of the format "X1981.01.01". But after writing and reading the netCDF file with new file <- stack("chirps81crop.nc") the layer names are changed to the format 'X1' up to 'X365'. I think it should be fine working with it, assuming that the order of the layers didn't get mixed up but I don't understand what is happening to the layer names and if this happens because there is something wrong with the code.

avocado1
  • 253
  • 2
  • 8

2 Answers2

2

It's the writeRaster() function that is losing the layer names, not the crop operation. It is possible to use lower level ncdf functions to assign a numeric value (not a string unfortunately) to each layer which will then show up in the name of the layers after reading. Taking inspiration from the example here, I created some code that shows this.

library(ncdf4)
library(raster)
library(rgdal)

# Crop extent
crop_extent <- as(raster::extent(5.74, 5.75, 50.96, 50.97), "SpatialPolygons")
proj4string(crop_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# make a sample file
r <- raster(system.file("external/test.grd", package="raster"))
r.latlon <- projectRaster(r, crs = proj4string(crop_extent))
writeRaster(x=r.latlon, filename = 'test.nc', format = 'CDF', overwrite=TRUE)

# read the sample as a 2 layer stack and crop it
test <- stack('test.nc', 'test.nc')
writeRaster(test, 'teststack.nc', overwrite=TRUE, format='CDF')
testcrop <- crop(test, crop_extent)
names(testcrop)
# [1] "test.1" "test.2"

# write the cropped file and make the zname equal to Layer
writeRaster(testcrop, 'testcrop.nc', overwrite=TRUE, format='CDF', zname='Layer')
# open the cdf file directly
nc <- nc_open('testcrop.nc', write = T)
# give the layers numbers starting from 10 so 
# we can see them easily
layers = 1:nlayers(testcrop) + 10
layers
# [1] 11 12
ncvar_put(nc, 'Layer', layers)
nc_close(nc)

newtestcrop <- stack('testcrop.nc')
names(newtestcrop)
# [1] "X11" "X12"
nc <- nc_open('testcrop.nc', write = F)
layers = ncvar_get(nc, 'Layer')
layers
# [1] 11 12
nc_close(nc)

So it is possible to get names with numbers under your control when writing the raster, but I don't know enough about your environment to determine if this will help since it might be tricky to map the names you need to a single unambiguous number.

Andrew Chisholm
  • 6,362
  • 2
  • 22
  • 41
1

I hope you don't mind me offering a non-R solution, but this task is much easier from the command line using CDO: 

cdo sellonlatbox,79,89,25,31 chirps-v2.0.1981.days_p05.nc cropped_file.nc

Which indices did you want to calculate? I suspect it is possible to calculate those quickly and easily with CDO functions too...

ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86
  • Hi Adrian, thanks for the answer. Mmh cdo does seem interesting, I will look into it. I have files though for the entire timeframe of the Chirps data set, so it's around 40 multilayer raster files that have to be cropped, stacked and used to calculate SPIs for different time periods (at least 3 month and 6 month SPI) which finally I have to merge with a shapefile. – avocado1 Feb 27 '20 at 14:19
  • This would be the R equivalent: `chirps81 <- brick("chirps-v2.0.1981.days_p05.nc")` ; `chirps81crop <-crop(chirps1981, extent(79, 89, 25, 31), filename="chirps81crop.nc")` (not that bad either) – Robert Hijmans Feb 27 '20 at 22:31
  • I tried this before with the extra line of code: `writeRaster(chirps81, "chirps81crop.nc", format="CDF")` and it only wrote the first band of the layer. I think it has to be a raster stack in order for all bands to be written back to file. Anyways I did find a solution using `chirps_stack <- stack(chirps_list)`; `chirps_cropped <- crop(chirps_stack, crop_extent)` and then write it back to netCDF. Now I have a file with continuous layer names from 1 to 12783 which is the exact number of days from 1981 to end 2015, which is pretty much exactly what I need. – avocado1 Feb 27 '20 at 22:56