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.
You have
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.
That is not a good approach (unless all else fails), and explains the odd extent. What you should try first is
Or better use the terra package (the replacement of raster)
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).