Generating data files netCDF

This example generates files of daily station records from the Global Historical Climate Network

## Save the GHCN station data as a netCDF4 file
## Go through different elements to generate one netCDF file for each. 
## Read 50 stations at the time and append to avoid memory issues.

library(esd)
nmin <- 75                       ## Pick only stations with minimum 75 years of data
ele <- c(111,121,601)            ## Elemennts: Tmin, Tmax, Precip
SS <- select.station(src='ghcnd',nmin=nmin,ele=ele)
stids <- SS$station_id
ns <- length(stids)

print(paste('Put',dim(SS)[1],'stations world-wide in netCDF file')) 
cntrs <- rownames(table(SS$country))
cntrs <- gsub(" ",".",cntrs)
cntrs <- gsub("[","",cntrs,fixed=TRUE)
cntrs <- gsub("]","",cntrs,fixed=TRUE)
cntrs <- gsub(",",".",cntrs,fixed=TRUE)

is <- seq(1,ns,by=50)
if (max(is) < ns) is <- c(is,ns)

param <- tolower(as.character(ele2param(ele,src='ghcnd')[5]))
print(param)
print(paste(ns,"stations and",length(cntrs),'countries with more than',
            nmin,'years of data'))
fname <- paste(param,'ghcnd','nc',sep='.')
print(fname)
append <- file.exists(fname)
  
for (id in is) {
    iii <- seq(id,id+49,length=50)
    iii <- iii[iii <= ns]
    print('Read data');print(id)
   
    x <- try(station(SS[iii,]))
    
    if (!inherits(x,'try-error')) {
        print(loc(x))
          units <- switch(toupper(param),'SD'='cm',
                          'CC'='octas','RR'='mm/day','FX'='m/s',
                          'DD'='degree','FG'='m/s','PP'='hPa',
                          'SS'='hours','HU'='percent')
        attr(x,'unit') <- units

      ## Quality check
      if ( (min(x,na.rm=TRUE) < -999) | (max(x,na.rm=TRUE)>2000) ) {
        print("Detected suspect data")
        print(range(x,na.rm=TRUE))
        xc <- coredata(x); xc[xc < -999] <- NA
        xc[xc > 2000] <- NA; coredata(x) <- as.matrix(xc)
      }
    write2ncdf4(x,fname,it=seq(as.Date('1893-01-01'),as.Date('2017-12-1'),by=1),
                append=append,verbose=FALSE,stid_unlim=TRUE)
    print('added to netCDF file')      
    } else {
        print('Failed to get data from GHCN:')
        print(param)
        print(stids[iii])
    }
  }