1

I'm attempting to reproject a multiband raster image in R to Orthographic for export https://www.naturalearthdata.com/downloads/50m-raster-data/50m-cross-blend-hypso/

However when I project the raster I get the following error, and the result although reprojected is not the bands that I started out with

Error in if (value[1] != nrow(x) | value[2] != ncol(x)) { : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In `dim<-`(`*tmp*`, value = c(nr, nc)) :
  NAs introduced by coercion to integer range

This is how I load the raster

library(raster)

# Create RasterStack object
r <- stack('HYP_50M_SR/HYP_50M_SR.tif')

In the console this returns

coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 


        HYP_50M_SR.1 HYP_50M_SR.2 HYP_50M_SR.3
Min.              64           94           88
1st Qu.          199          199          193
Median           250          250          250
3rd Qu.          251          251          251
Max.             255          255          255
NA's               0            0            0

Then I reproject

# Define the Proj.4 spatial reference 
sr <- "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"

# Project Raster
projected_raster <- projectRaster(r, crs = sr)

In the console this returns

coord. ref. : +proj=ortho +lat_0=-10 +lon_0=-60 +ellps=WGS84 

      HYP_50M_SR
Min.    7.105387e+01
1st Qu. 2.027389e+02
Median  2.502865e+02
3rd Qu. 2.507040e+02
Max.    2.550000e+02
NA's    4.880675e+06

How can I reproject a multiband raster while keeping the band information?

Ive been using this for reference https://datacarpentry.org/r-raster-vector-geospatial/05-raster-multi-band-in-r/

2

Firstly, you must be looking at a projected_raster created earlier, because on that error R will not create an object, and what's more the projection strings don't match. That explains why your output looks wrong (because its not the output from reprojecting the raster) but not the error.

So we read the help for projectRaster:

     projectRaster(from, to, res, crs, method="bilinear", 
                  alignOnly=FALSE, over=FALSE, filename="", ...) 

You are doing projectRaster(r, sr), so its matching r as from and sr as to. What are those?

    from: Raster* object
    to: Raster* object with the parameters to which 'from' should be projected

But your to is a CRS string, not a raster object. You only use this unnamed parameter form of rasterProject when you have a second raster and you want to project the first one to that one. Instead you have a CRS string so you need name the crs parameter to get the matching right: projectRaster(r, crs=sr). You should also specify a resolution or you might end up with an output grid of the world in metres which might be a bit big.

Here's a worked example (but I'm cropping the source raster down a bit since I want this to run before I have to go to work)

Crop the large source to a 1x1 degree chunk

> rc = crop(r,extent(45,46,23,24))

Try to project and fail because we didn't name the argument:

> rp = projectRaster(rc,sr)
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘res’ for signature ‘"character"’
In addition: Warning message:
In min(dim(to)[1:2]) : no non-missing arguments to min; returning Inf

and notice that rp does not exist:

> rp
Error: object 'rp' not found
> 

Project it to a CRS with a named argument and get an RGB stack back:

> rp = projectRaster(rc,crs=sr,res=1000)
> rp
class      : RasterBrick 
dimensions : 153, 165, 25245, 3  (nrow, ncol, ncell, nlayers)
resolution : 1000, 1000  (x, y)
extent     : 8023238, 8188238, -5279.056, 147720.9  (xmin, xmax, ymin, ymax)
crs        : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs 
source     : memory
names      : HYP_50M_SR.1, HYP_50M_SR.2, HYP_50M_SR.3 
min values :     187.1538,     177.0546,     160.1354 
max values :     219.6638,     208.5521,     191.2164 

Plot them:

> plotRGB(rc)
> plotRGB(rp)
> 
  • I'd been thrown off by the error I was getting but it was really quite simple! Thank you for this thorough answer, you've really helped me out. Setting the resolution is a great tip too – user3821345 Jul 16 at 10:17

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.