Search code examples
rfunctionsystemgdal

R function stops after system() call

I've written a very easy wrapper around GDAL in R. It utilises a prewritten statement which is passed to system, creating an output, which I then want to read into the R environment again.

It works by creates a temporary directory in the working directory, printing out an ESRI shape file of our area of interest, and then cuts a raster by this, with some preset information.

My problem: after successfully running the system() call and creating the output file, the function stops. It doesn't execute the next call and read the output into the R environment.

gdalwarpwarp <- function(source_file, source_srs, newfilename, reread=TRUE, clean=TRUE, cutline){

  #Create tempfolder if it doesn't exist in the working directory.
  if (!dir.exists("./tempfolder")){
    dir.create("./tempfolder")
  }

  #Write temporary shape file
  terra::writeVector(cutline, filename = "./tempfolder/outline_AOI.shp" , filetype='ESRI Shapefile',overwrite=TRUE)

  #Warp!

  if(reread==FALSE){
    system(paste0("gdalwarp -cutline ./tempfolder/outline_AOI.shp -dstnodata -9999 -s_srs EPSG:3006 ",source_file, " ",paste0("./tempfolder/",newfilename)))
    message('warp warped TRUE')

  } else if(reread==TRUE){
    system(paste0("gdalwarp -cutline ./tempfolder/outline_AOI.shp -dstnodata -9999 -s_srs EPSG:3006  ",source_file, " ",paste0("./tempfolder/",newfilename)))
    newfilename <- terra::rast(paste0("./tempfolder/",newfilename))
  }

}

This doesn't run:

newfilename <- terra::rast(paste0("./tempfolder/",newfilename))
Solution

  • The function did not return anything. Here is a somewhat improved version of your function. If you want to keep the output it would make more sense to provide a full path, rather then saving it to a temp folder. I also note that you are not using the argument source_srs

    gdalwarpwarp <- function(source_file, source_srs, newfilename, reread=TRUE, clean=TRUE, cutline){

  #Write temporary shape file
  shpf <- file.path(tempdir(), "aoi.shp")
  terra::writeVector(cutline, filename = shpf, filetype='ESRI Shapefile',overwrite=TRUE)

   outf <- file.path(tempdir(), newfilename)
    system(paste0("gdalwarp -cutline shpf -dstnodata -9999 -s_srs EPSG:3006 ",source_file, " ", outf)
   
    if (reread) {
        terra::rast(outf)
    } else {
       message('warp warped TRUE')
       invisible(filename)
    }
}

    I wonder why you don't use terra::resample or terra::project; perhaps preceded or followed by mask (I may not understand the benefit of using cutline.