Modalities of conversion from NCDF to raster with R

Question

I am using environmental data layers from Copernicus and as they are available in NCDF format I'm converting them in raster to be able to manipuplate them with QGIS. I use the following script :

## Output: a GeoTIFF file
file.tiff <- 'currents_2021.tiff'

## Import netCDF
currents_2021 <- raster(file.nc)

## Save to disk as GeoTIFF
writeRaster(currents_2021, filename = file.tiff, format = 'GTiff', overwrite = T)

## For multiple files, could use a for loop
## Input directory
dir.nc <- '~/ENSAT/Mobilité/Stage 2A/Archipelagos/Seabirds project/Environmental variables/Currents/2021'
files.nc <- list.files(dir.nc, full.names = T, recursive = T)

## Output directory
dir.output <- '~/ENSAT/Mobilité/Stage 2A/Archipelagos/Seabirds project/Environmental variables/SST/requested_files'

It works but the problem is that I have datas that are over 1 year (and collected every day), so for each latitude and longitude, several values of the variable considered are associated - theorically 365, one per day. I don't know how the raster function deals with that and which value it takes for each coordinate. I wouldlike to have an average of the year, so in each case of the raster an average of all the value associated, but I don't know if this script does it automatically. Does anyone knows that ? Thanks ! :)

dimfalk · Answer 1 · 2022-07-28T08:11:45.900

?raster illustrates that raster::raster() has a band argument: Integer. The layer to use in a multi-layer file.

Since you are reading netCDF files (with 365 layers in your case), it seems like raster::raster() is only able to read one layer at once. If you don't provide any band argument, band = 1 is chosen per default.

library(raster)
library(ncdf4)

raster::raster("chirps-v2.0.2006.days_p25.nc")
#> class      : RasterLayer 
#> band       : 1  (of  365  bands)
#> dimensions : 400, 1440, 576000  (nrow, ncol, ncell)
#> resolution : 0.25, 0.25  (x, y)
#> extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : chirps-v2.0.2006.days_p25.nc 
#> names      : Climate.Hazards.group.InfraRed.Precipitation.with.Stations 
#> z-value    : 2006-01-01 
#> zvar       : precip

raster::raster("chirps-v2.0.2006.days_p25.nc", band = 40)
#> class      : RasterLayer 
#> band       : 40  (of  365  bands)
#> dimensions : 400, 1440, 576000  (nrow, ncol, ncell)
#> resolution : 0.25, 0.25  (x, y)
#> extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : chirps-v2.0.2006.days_p25.nc 
#> names      : Climate.Hazards.group.InfraRed.Precipitation.with.Stations 
#> z-value    : 2006-02-09 
#> zvar       : precip

The object returned is a single-band RasterLayer. You can import a multi-layered file (such as data coming from a netCDF container) at once making use of raster::stack() instead:

nc_data <- raster::stack("chirps-v2.0.2006.days_p25.nc")
nc_data
#> class      : RasterStack 
#> dimensions : 400, 1440, 576000, 365  (nrow, ncol, ncell, nlayers)
#> resolution : 0.25, 0.25  (x, y)
#> extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> names      : X2006.01.01, X2006.01.02, X2006.01.03, X2006.01.04, X2006.01.05, X2006.01.06, X2006.01.07, X2006.01.08, X2006.01.09, X2006.01.10, X2006.01.11, X2006.01.12, X2006.01.13, X2006.01.14, X2006.01.15, ...

As you can see, the result consists of 365 layers this time. You can now subset the RasterStack object to query individual layers if you like, ...

nc_data[[40]]
#> class      : RasterLayer 
#> band       : 40  (of  365  bands)
#> dimensions : 400, 1440, 576000  (nrow, ncol, ncell)
#> resolution : 0.25, 0.25  (x, y)
#> extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : chirps-v2.0.2006.days_p25.nc 
#> names      : X2006.02.09 
#> zvar       : precip

... and you can now simply make use of raster::mean() on your stack to calculate average values over all 365 layers, corresponding to the yearly average you requested, and write the resulting single-band raster to disk using raster::writeRaster():

nc_mean <- raster::mean(nc_data)
nc_mean
#> class      : RasterLayer 
#> dimensions : 400, 1440, 576000  (nrow, ncol, ncell)
#> resolution : 0.25, 0.25  (x, y)
#> extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : memory
#> names      : layer 
#> values     : 0, 71.93983  (min, max)

If you want to find out how to read and write netCDF files using terra, see my other post. Just make sure to change the file extension in terra::writeRaster() from asc to tif.

Modalities of conversion from NCDF to raster with R

1 Answers1