Reputation: 199
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))
Upvotes: 0
Views: 74
Reputation: 47146
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
.
Upvotes: 1