1

I have a few .tif rasters (ELEV and POP) projected in long/lat, but I would like the CRS and extent to be the same as a shapefile I have (STUDYAREA), for example, my study area lies mostly in UTM16 so I want that to be in the crs for all objects. Also, I would like them all to have a 1000*1000m resolution.

I first make sure that the extent is the same as the shapefile, as I've learned that working with objects with a different extent will give an error. Here is how I try to reach my goal:

create object to adapt other objects extent to (needs to be equal to extent of STUDYAREA)
studyarea <- spTransform(studyarea,CRSobj= "+proj=utm +zone=16 +datum=WGS84 +units=m + no_defs")
areaextent <- extent(studyarea)
projection <- "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"

elev <- raster("Gtopo30/Gtopo30.TIF")
newelev <- setExtent(elev,areaextent)
elev <- projectRaster(newelev,res=1000,crs=projection)


pop <- raster("griddedpop/gpw.tif")
pop <- setExtent(elev,areaextent,keepres=FALSE)
pop <- projectRaster(pop,crs=projection)

But when using the ProjectRaster function, I run into very similar errors for both raster projections:

elev <- projectRaster(newelev,res=1000,crs=projection)
Error in projectExtent(from, projto) : cannot do this transformation
In addition: Warning message:
In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
704 projected point(s) not finite
Improve this question

1 Answer 1

Reset to default
1

Here is an improved version of your code. As you do not provide example data I cannot say if it works better.

library(raster)
projection <- "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"
studyarea <- spTransform(studyarea, CRSobj=projection)

# make a template RasterLayer
template <- raster(studyarea, res=1000)

# use the template 
elev <- raster("Gtopo30/Gtopo30.TIF")
elevutm <- projectRaster(elev, template)

You can also use projectRaster without a template like this:

elevutm2 <- projectRaster(elev, crs=projection, res=1000)

But I do not recommend that, as you have less control that way.

Below the same workflow, but with example data and using the terra package (the replacement of raster)

library(terra)
prj <- "+proj=utm +zone=32 +units=m"
v <- vect(system.file("ex/lux.shp", package="terra"))

sa <- project(v, prj)
tmp <- rast(sa, resolution=1000)
# to get a nicer extent 
ext(tmp) <- floor(ext(r)+5000)

elevation <- rast(system.file("ex/elev.tif", package="terra"))
elevprj <- project(elevation, tmp)

plot(elevprj)
lines(sa)

enter image description here

Improve this answer
7
  • Thank you so much Robert! It works, yetI am quite surprised that elevutm2 <- projectRaster(elev, crs=projection, res=1000) works, after all, isn't this exactly what I did, except for that the name of the assigned object is not a pre-existing object? My line being elev <- projectRaster(newelev,res=1000,crs=projection) with 'elev' being an already existing object to be replaced by this new code. Why would this make a difference and not have me run into the error?
    – Beardedant
    Commented Jul 7, 2021 at 23:17
  • I think your mistake is newelev <- setExtent(elev, areaextent) That would seem to assign a UTM extent to a lonlat raster; creating invalid coordinates (latitude >> 90 or << -90)
    – Robert Hijmans
    Commented Jul 7, 2021 at 23:26
  • That makes sense! I think in general my trouble is with the following: When working with projectRaster I have in the past gotten an error when trying to project raster 1 to raster 2 when these are of a different extent. So I try to equal the extent of my rasters before using projectRaster, but this leads to trouble when the rasters are in a different CRS (specifically a UTM vs longlat). So I guess what I'm looking for is the most straightforward way to equalize a raster to another one in terms of CRS, extent and resolution?
    – Beardedant
    Commented Jul 8, 2021 at 0:04
  • When projecting from r1 to r2, the extents are (almost) always different, so that was probably not the real issue. But you can try projectExtent and first project an empty raster (with no values) to see where it ends up and adjust that.
    – Robert Hijmans
    Commented Jul 8, 2021 at 0:38
  • 1
    project or projecRaster will only read what it needs to read, so the size of the input should not matter. Use terra for much better performance.
    – Robert Hijmans
    Commented Jul 8, 2021 at 1:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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