[GIS] stack rasters with different extent and resolution in a loop with R package raster

coordinate systemrrasterstack

I would like to stack rasters with the R package "raster". Layers are in ascii format. I want to crop them to have the same extent and project them to have the same resolution. My code gives raster layers with empty values. 

#The layers must have a similar extent and projection. We will first crop them all to the same extent.
#Read files found in the working directory containing *.asc
raster_data <- list.files(pattern='*.asc', full.names=TRUE)
#get desired extent (smallest raster extent)
e <- extent(-12, 112, 23, 72)
#get template layer with desired extent and resolution. Transform ascii template layer to raster format
d <- raster("WC_BIO1.asc", ext = e, crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0", ncol=15120, nrow=6000)
w <- writeRaster(d, "WC_BIO1_template.grd", format = "raster")
#load ascii layers, transform to raster format with extent and resolution obtained from template layer
for (i in 1:length(raster_data)){
  r <- raster(raster_data[i], crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  rc <- crop(r, e)
  rp <- projectRaster(from = rc, to = w,
                      crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
                      filename = file.path ("./crop", raster_data[i]),
                      method = "bilinear",
                      format = "raster",
                      overwrite = TRUE)
  rw <- writeRaster(rp,
                    filename = file.path ("./crop", raster_data[i]),
                    crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
                    overwrite = TRUE,
                    format = "raster",
                    na.value=-9999)
  }

Best Answer

Here are some hints. It is difficult to answer without knowing more about your data. Please show us the result of:

raster(raster_data[1])

  • You pass arguments to the raster function that are ignored when you create an object from file. You do not need to use crop on this layer. Certainly not before you project the data to lon/lat.

  • You never need projectRaster AND resample. If your layers are indeed lon/lat use resample or, if possible, just crop and (dis)aggregate.

  • My impression is that they are not lon/lat. This may be the problem and the reason that your layers are empty. In that case you must assign the correct current CRS to them, before using projectRaster. In projectRaster do not use the crs argument (it is ignored anyway, as the CRS from w is used).

  • Do not use writeRaster, as you are already saving the file with projectRaster (or could do so with resample). Also writeRaster does not have arguments crs and na.value, so these are ignored.

Here is the fixed code:

raster_data <- list.files(pattern='\\.asc$', full.names=TRUE)
w <- raster("WC_BIO1.asc", crs="+proj=longlat +datum=WGS84 +ellps=WGS84")

for (i in 1:length(raster_data)){
    r <- raster(raster_data[i], crs="???????")
    rp <- projectRaster(from = rc, to = w,
                      filename = file.path ("./crop", raster_data[i]),
                      method = "bilinear",
                      format = "raster",
                      overwrite = TRUE)
}
Related Question