0

I have a list of rasters in a folder that I want to project, rescale values and clip. My current code looks like this

library(rgdal)
library(maptools)
library(raster

setwd("D:/working")

# define projections
utm33s <- CRS('+proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')

# read in and stack files
r.list <- list.files(pattern="*.tif", full.names=TRUE)
r.stack <- stack(r.list)

# project stack to utm 33s
r.stack.prj = projectRaster(r.stack, crs = utm33s)

# function to calculate float values
fun <- fun <- function(x) {
  y = (x - 50) / 200
  return(y)
}

# apply function on the stack
calc(r.stack.prj, fun=fun)

# read in polygon file for use in clipping
aoi = readOGR(dsn=".", layer="aoi_shape")

# clip raster using polygon
r.stack.prj <- crop(r.stack.prj, aoi, snap="out") 

# save as GeoTIFF
writeRaster(r.stack.prj, 'clean_stack.tif', overwrite=TRUE)

The script crashes at 'projectRaster' What I'm I doing doing wrong?

  • What is the exact error you are getting? It appears that the initial rasters might not have a coordinate system assigned to start with in the first place. – Vijay Ramesh Jun 5 '17 at 18:57
  • @Vijay-ramesh The rasters have this coordinate system +proj=sinu +lon_0=20 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs assigned. – jmutua Jun 5 '17 at 20:03
  • As already commented, difficult to help unless you report also which specific error you are getting. – lbusett Jun 5 '17 at 21:03
  • can you provide details of one file, i.e. the printed output of 'raster(r.list[1])' - my feeling is this sinu projection is probably of far greater extent than the valid range for this UTM zone. Generally reprojection of rasters is ill-advised, usually it's better to calculate in the native projection of the gridded data, and transform other vector data to that space for extractions etc. Also UTM is a very poor choice, in the general case - it's probably unnecessary and possibly a bad choice here. – mdsumner Jun 6 '17 at 9:26
0

I suggest doing this instead:

r.stack <- stack(r.list)
calc.stack <- calc(r.stack.prj, fun=fun)
aoi = readOGR(dsn=".", layer="aoi_shape")
r.stack.crop <- crop(calc.stack, spTransform(aoi, projection(calc.stack), snap="out") 

Then if really necessary, projectRaster the r.stack.crop - don't do it before the other operations, it's merely wasteful and destructive (in the technical senses).

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.