Converting large number of NetCDFs to single raster using R

Question

There are ~87,500 NetCDF files of a global dataset that must be stacked and the values averaged to create 1 raster. My machine does not like this... I've written nested for loops that, despite strategically placing rm()'s and gc()'s, still cripple my machine and cause the conversion process to halt. Here is a link to 10 on which the code should work fine: https://www.dropbox.com/scl/fo/oayqkfxv01f24h4b29s9p/h?dl=0&rlkey=g2m031253uf7ngi0o1jsb1epz

The code chunk that follows is what converts the netcdfs to rasters and is at the core of the messier code with for loops at the bottom of this post.

require(ncdf4)
require(maps)
require(rgdal)
require(raster)
require(terra)


dir<-"/your working directory where the netcdfs from the dropbox link are stored/"
files<-list.files(path=dir, pattern='.nc', full.names = TRUE)
setwd(dir)

for(k in 1:length(files)){
print(k)
  if(k==1){
    nc = nc_open(file)
    total_par = ncvar_get(nc, "total_par")
    z = rast(t(total_par))
    ext(z)=c(-180, 180, -90, 90)
    crs(z) = "epsg:4326"
    temp <- brick(z)
  }
  if(k!=1){
    nc = nc_open(file)
    total_par = ncvar_get(nc, "total_par")
    z = rast(t(total_par))
    ext(z)=c(-180, 180, -90, 90)
    crs(z) = "epsg:4326"
    temp2 <-brick(z)
    temp <- stack(temp,temp2)
  }
}

What I am hoping someone can provide is a solution for converting a large number of NetCDFs into one raster that does not consume so much memory and cause R to terminate. I think loops (at least how I've written them) are probably not the best way forward. Is anyone aware of a more eloquent technique or a way to fix my crashing loops?

Below is is what code I've written for this task. Note, you'll need to recreate the working directory architecture (I can add that in for you if requested). Also, here is a link to about 5000 NetCDFs (if you're feeling bold): https://www.dropbox.com/scl/fo/t5q2xjhixkn3najmnxcv7/h?dl=0&rlkey=5i5fuuuzie49f4zx5ng45ahxo - note the code below works for about 3000 NetCDFs and then R crashes:

# this code calls in ncdf files containing 3 hr daily PAR measurements for the global surface at 10km horizontal resolution. Once called, the netcdfs are converted to rasters, stacked, and then averaged to create a global PAR surface over the time period 1982-2018. Because of memory constraints on my machine, I had to add some loop structures that saved the netcdf to tifs then purged the memory. The loops index where in the process the next iteration needed to pick up from.
require(ncdf4)
require(maps)
require(rgdal)
require(raster)
require(terra)
years <- seq(1984,2018) # this is the year range over which the data was collected and used to call data organized in files by year
quot = 0 # this is a step trigger that causes the loop to progress between year steps - it only increases once all NetCDFs for a year have been stacked, averaged, converted to a raster, and saved locally. Quot is determined by dividing the count of csvs created iteratively by the for loop and stored in a special working directory. When quot = 1, the year step progresses forward one (at least that was the goal) because the while loop is exited. After this progression, the working directory is deleted and the build up of quot begins again.
file_start_k = 0 # this is used to select the proper starting file following the data and garbage dump
yr_ct_dir<-"/Volumes/WD_BLACK/PAR/counts/" # csvs are saved in this file to signify the progression of years through the loop iterations
for(year in years){
#find all ncs in your directory
dir<-paste("/Volumes/WD_BLACK/PAR/",year,sep="")

#get a list of all files with total_par in the name in your directory
files<-list.files(path=dir, pattern='.nc', full.names = TRUE)
#for calls
setwd(dir)
k = 1 # organization index
# this cycles through all files of indicated extension and pulls necessary data
while(quot < 1){
for (file in files){ #this loop translates netcdf stored info into rasters
 
 if(file_start_k == 0){
    print(k) #sometimes, you just gotta know...
    files <- files[(file_start_k+1):500] 
  }
  if(file_start_k > 0 && file_start_k < 2000){ # this determines which file within files to start with following the garbage dump
    files <- files[(file_start_k+1):(file_start_k+500)] 
  }
  if(file_start_k == 2000){
    files <- files[(file_start_k+1):length(files)] 
  } 
  if(k==1){
    nc = nc_open(file)
    total_par = ncvar_get(nc, "total_par")
    z = rast(t(total_par))
    ext(z)=c(-180, 180, -90, 90)
    crs(z) = "epsg:4326"
    temp <- brick(z)
  }
  if(k!=1){
    nc = nc_open(file)
    total_par = ncvar_get(nc, "total_par")
    z = rast(t(total_par))
    ext(z)=c(-180, 180, -90, 90)
    crs(z) = "epsg:4326"
    temp2 <-brick(z)
    temp <- stack(temp,temp2)
  }
  k = k + 1
  if (k > length(files)){ # this loop is necessary because of memory constraints on my machine
    #removes NAs and averages rasters by year
    write.csv(1,paste(yr_ct_dir,Sys.time(),".csv",sep="")) # this acts as within year index as reference for next iteration year used
    write.csv(1,paste(dir,"/",Sys.time(),".csv",sep="")) # this acts as within year index to reference starting file
    rs1 <- calc(temp,mean,na.rm=TRUE) # this creates a single tif which is the average of all files in the 500 stack
    writeRaster(rs1,paste(dir, "/",year,"_",file_start_k,".tif",sep="")) # this saves the tif locally
    # have to clear memory. so, built indexing system to recall where the process was
    # clean up
    rm(list=c("temp","temp2","total_par","z","nc", "rs1"))
    gc()
    gc()
    gc()

    # this chunk of black magic points the code to the proper directories and files to pick up from where the garbage dump above occurred.
    year_csvs<-list.files(path="/Volumes/WD_BLACK/PAR/counts/", pattern='.csv', full.names = TRUE) # this calls the above csvs for counting
    quot <- trunc(length(year_csvs)/5,0) #number of counts divided by 5. after 5, the year advances by one. removing the remainder provides the additive factor that modifies the minimum year of the year range such that it progresses every 5th 500th iteration
    # year = year + quot  # takes above additive factor and corrects age when necessary
    dir <- paste("/Volumes/WD_BLACK/PAR/",year,sep="")
    files_csvs <- list.files(path=dir, pattern='.csv', full.names = TRUE)
    file_start_k <- length(files_csvs) * 500 # used to strip the files already ready from the files variable used to call each file 
    k = 1
    
    if(quot == 1){
      unlink("/Volumes/WD_BLACK/PAR/counts/*.csv")
    }
    }
  
}
  }
}

Answer 1

Try this. It uses the online mean computation algorithm and doesn't bother to convert the data to raster since R can work with the matrix data straight out of NetCDF - all the NetCDFs are assumed to be on an identical grid. It also closes the NetCDF handle once its been used.

mean_nc <- function(files){
    n = 0
    first=TRUE
    for(file in files){
        nc = nc_open(file)
        total_par = ncvar_get(nc, "total_par")
        nc_close(nc) # release
        if(first){
            m = matrix(0, nrow=nrow(total_par), ncol=ncol(total_par))
            first = FALSE
        }
        n = n + 1
        m = (1/n) * total_par + ((n-1)/n) * m
    }
    return(m)
}

Note: I've realised the "first" mechanism here which I wrote to setup a zeros matrix with the right shape inside the loop on the first instance isn't required like that. If you have a simple m=0 outside the loop then on the first calculation it will get added to the values matrix and the updated m will suddenly become a matrix of the right dimension. I've not tested this and there's not much of an efficiency gain in doing it compared to the existing code, but I thought I'd mention it in case someone else noticed.

I've tested it on four from your download list:

> ncs = list.files(".",pattern=".nc$")
> mn = mean_nc(ncs)
> image(mn)

(white lines are plotting artefacts)

and I would suggest doing some hand-calculation on a few pixels to make sure the mean is correct. This also assumes there's no missing data values. There's a lot of zeroes...

You might also want to add a progress bar with utils::txtProgressBar or otherwise see how its going...

On my 2014 PC (which was speedy then) it takes about 36s to process 400 files (actually 100 copies of each of the 4 I've downloaded), which works out at about 2 hours for 87,500...