Looping through multiple NetCDF files, rotate, crop and update their crs-CodePudding

I have multiple NetCDF files that contain fire weather data. The files are rotated and also, other R packages (e.g., stars) and programs (e.g., cdo) do not recognize their crs but I can see the crs as WGS84 when I opened them with the raster package.

I'm looking into rotating the data to the common/regular lonlat structure, update the crs, crop the data to the extent of my study area, and writing out each NetCDF file (containing 365/366 layers) just as they were loaded. Example data here: https://www.filemail.com/t/JtIAB1zC

Here's the code that I'm working with:

 all_nc = list.files("C:/file_path/", 
           pattern = ".nc$", recursive = F, full.names = T)
    
    for(i in seq_along(all_nc)){ 
          
        r = stack(all_nc[i])
        r2 = raster::rotate(brick(r))
        
        projectRaster(r2, crs = crs(" proj=longlat  datum=WGS84  ellps=WGS84  towgs84=0,0,0")) 
          
        r2 <- terra::crop(r2, terra::extent(-141.0069, -123.7893, 60, 69.64794))
          
          writeRaster(r, file.path('D:/path/', names(r)), 
                      force_v4 = T, overwrite=TRUE, format="CDF", compression = 7)
          
          }

Error in if (filename == "") { : the condition has length > 1

I can get individual layers with all the changes that I want by tweaking the code but that's not what I'm hoping to achieve. I don't use loops often. So, I suspect that I'm screwing up the indexing somehow.

CodePudding user response：

You can do that like this

fin <- list.files("C:/file_path/", pattern = ".nc$", full.names=TRUE)
fout <- gsub("C:/file_path/", "D:/path/", fin)

library(terra)    
e <- ext(-141.0069, -123.7893, 60, 69.64794)

for(i in seq_along(fin)){     
    r <- rast(fin[i])
    crs(r) <- " proj=longlat  datum=WGS84"
    r <- rotate(r)   
    r <- crop(r, e)
    writeCDF(r, filename=fout[i], compression=7)
}

Instead of

r <- crop(r, e)
writeCDF(r, filename=fout[i], compression=7)

writeCDF is a better netCDF writer than writeRaster, but you could also use

r <- crop(r, e, filename=fout[i], gdal=c(COMPRESS="DEFLATE", ZLEVEL="7"))
## or 
r <- crop(r, e)
writeRaster(r, filename=fout[i], gdal=c(COMPRESS="DEFLATE", ZLEVEL="7"))

Note that you are mixing up a lot of "raster" and "terra" code. In the above I only use "terra". Also could you point to or share one of the files. I would like to find out why GDAL does not detect the crs.