Error calculating area of raster with lat/lon projection

Question

I have a global raster stack (of three rasters) whose pixel values are the percent of a land use for that pixel. Here's the raster metadata:

class      : RasterBrick 
dimensions : 3600, 7200, 25920000, 3  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 7200, 0, 3600  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : grass_baseline.tif 
names      : grass_2020, grass_2040, grass_2100

I'm trying to calculate the total area of land use in each pixel by multiplying the pixel value by the area of the raster, using the area() function in the raster package.

When I do that, I get the following error:

Warning message:
In .couldBeLonLat(x, warnings = warnings) :
  raster has a longitude/latitude CRS, but coordinates do not match that

Here's the metadata for the area raster:

class      : RasterLayer 
dimensions : 3600, 7200, 25920000  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 7200, 0, 3600  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : -710.0924, 2211922  (min, max)

Does anyone have any insight into what might be going on?

In case it's relevant, I assembled this raster stack from a few .nc files that I read into R with the ncdf4 package and converted to rasters with the following line of code:

raster(first_nc, xmn=0, xmx=7200, ymn=0, ymx=3600, crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0")

I then combined several of these rasters together as a stack and exported using the stars package (to preserve the names of each raster):

stack <- stack(first_nc,second_nc,third_nc)
names(stack) <- c('first_nc','second_nc','third_nc')
stars::write_stars(stars::st_as_stars(stack), "stack.tif")

I then read the .tif into a separate script, which is where I'm trying to calculate the area.

Robert Hijmans · Answer 1 · 2022-05-24T17:43:40.530

You have

#extent     : 0, 7200, 0, 3600  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs

That is, a latitude between 0 and 3600 degrees. That makes no sense as you cannot go beyond 90 degrees N and it is thus not possible to compute area for these cells. And the specified longitude is not likely to be correct either, unless your data really covers the globe 20 times.

I assembled this raster stack from a few .nc files that I read into R with the ncdf4 package and converted to rasters

That is not a good approach (unless all else fails), and explains the odd extent. What you should try first is

library(raster)    
s <- stack(filenames)

Or better use the terra package (the replacement of raster)

library(terra)    
s <- rast(filenames)

It should not be necessary, but if you are going to set the extent yourself, more plausible values would be (-180, 180, -90, 90), or (0, 360, -90, 90).

Error calculating area of raster with lat/lon projection

1 Answers1