1

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)
  }
  • 1
    resample from the raster package has been designed for this particular purpose. Inside the loop, use rp <- projectRaster(r, w) before running resample(rp, w). As far as I can tell, the use of crop is not required. – fdetsch Dec 8 '15 at 16:19
  • @fdetsch is right on the mark. You do not want to use projectRaster. This funciton is for transforming the projection space of the raster and not subsetting or resampling. Also, wildcards using "*" do not produce the expected results in list.files and you are missing the path argument. This should be a regular expression and would looks like: list.files(path = getwd(), pattern = "asc$") – Jeffrey Evans Dec 8 '15 at 19:33
2

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)
}
  • I didn't know projectRaster includes spatial resampling when 'to' is specified. Thanks for the hint! – fdetsch Dec 9 '15 at 8:38
  • Here is the result of raster(raster_data[1]) class : RasterLayer dimensions : 199, 503, 100097 (nrow, ncol, ncell) resolution : 0.25, 0.25 (x, y) extent : -13, 112.75, 23, 72.75 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : ... names : CMSAF_SDL > As you can see my dataset do not have crs so I used: crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" ...@fdetsch solution also worked well, but I guess that your code is the right way to do it. – jl-blancopastor Dec 9 '15 at 16:43
  • No, in that case projectRaster is not the way to go, but you should indeed use resample; or a combination of disaggregate and crop – Robert Hijmans Dec 9 '15 at 22:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.