I want to crop an elevation raster to add it to a raster stack. It's easy, I did this before smoothly, adding a ecoregions raster to the same stack. But with the elevation one, just doesn't work. Now, there are several questions here in overflow adressing this issue and I tryed a lot of things...
First of all, we need this:
library(rgdal)
library(raster)
My stack is predictors2:
#Downloading the stack
predictors2_full<-getData('worldclim', var='bio', res=10)
#Cropping it, I don' need the whole world
xmin=-120; xmax=-35; ymin=-60; ymax=35
limits <- c(xmin, xmax, ymin, ymax)
predictors2 <- crop(predictors2_full,limits)
Then I've downloaded the terr_ecorregions shapefile here: http://maps.tnc.org/files/shp/terr-ecoregions-TNC.zip
setwd("~/ORCHIDACEAE/Ecologicos/w2/layers/terr-ecoregions-TNC")
ecoreg = readOGR("tnc_terr_ecoregions.shp") # I've loaded...
ecoreg2 <- crop(ecoreg,extent(predictors2)) # cropped...
ecoreg2 <- rasterize(ecoreg2, predictors2) # made the shapefile a raster
predictors4<-addLayer(predictors2,elevation,ecoreg2) # and added the raster
# to my stack
With elevation, I just can't. The Digital elevation model is based in GMTED2010, which can be downloaded here: http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/mn30_grd.zip
elevation<-raster("w001001.adf") #I've loaded
elevation<-crop(elevation,predictors2) # and cropped
But elevation gets a slightly different extent instead of predictors2's extent:
> extent(elevation)
class : Extent
xmin : -120.0001
xmax : -35.00014
ymin : -60.00014
ymax : 34.99986
>
I tried to make then equal by all means I read about in questions here... I tried to extend so elevation's ymax would meet predictors2's ymax elevation<-extend(elevation,predictors2) #didn't work, extent remains the same
I tried the opposite... making predictors2 extent meet elevation's extent... nothing either.
But then I read that
You might not want to play with setExtent() or extent() <- extent(), as you could end with wrong geographic coordinates of your rasters - @ztl, Jun 29 '15
And I tried to get the minimal common extent of my rasters, following @zlt answer in another extent question, by doing this
# Summing your rasters will only work where they are not NA
r123 = r1+r2+r3 # r123 has the minimal common extent
r1 = crop(r1, r123) # crop to that minimal extent
r2 = crop(r2, r123)
r3 = crop(r3, r123)
For that, first I had to set the resolutions:
res(elevation)<-res(predictors2) #fixing the resolutions... This one worked.
But then, r123 = r1+r2+r
didn't work:
> r123=elevation+ecoreg2+predictors2
Error in elevation + ecoreg2 : first Raster object has no values
Can anyone give me a hint on this? I really would like to add my elevation to the raster. Funny thing is, I have another stack named predictors1 with the exact same elevation's extent... And I was able to crop ecoreg and add ecoreg to both predictors1 and predictors2... Why can't I just do the same to elevation? I'm quite new to this world and runned out of ideas... I appreciate any tips.
EDIT: Solution, Thanks to @Val
I got to this:
#Getting the factor to aggregate (rasters are multiples of each other)
res(ecoreg2)/res(elevation)
[1] 20 20 #The factor is 20
elevation2<-aggregate(elevation, fact=20)
elevation2 <- crop(elevation2,extent(predictors2))
#Finally adding the layer:
predictors2_eco<-addLayer(predictors2,elevation2,ecoreg)
New problem, thought...
I can't write stack to a geotiff
writeRaster(predictors2_eco, filename="cropped_predictors2_eco.tif", options="INTERLEAVE=BAND", overwrite=TRUE)
Error in .checkLevels(levs[[j]], value[[j]]) :
new raster attributes (factor values) should be in a data.frame (inside a list)
options="INTERLEAVE=BAND"
and usebylayer = F
– aldo_tapia