Computing weights from latitude to apply to raster stack values

Question

I am currently trying to apply weights to the original precipitation values of my Raster Stack object, called "landCO2". These weights are to take into account the differences in area between the equator and poles. However, I am not sure how to approach applying these to the existing values of the Raster Stack. The idea would be to apply the weights to the original values of each of the 138 raster layers, and then eventually plot these values on a global map with wrld_simpl. Here is what was done so far:

library(raster)
library(maps)
library(maptools)
library(rasterVis)


data("wrld_simpl")
b <- wrld_simpl

landCO2 <- mask(RCP1pctCO2Median, b)
CO2new <- rasterToPoints(landCO2)
weightCO2 <- cos(CO2new[,"y"]*(pi/180))
CO2new[,3:ncol(CO2new)] = apply(CO2new[,3:ncol(CO2new)], 2, function(x) x * weightCO2)
avgCO2 <- colSums(CO2new[,3:ncol(CO2new)])/sum(weightCO2)

However, this approach obtains an average across all grid cells per layer, effectively creating 138 averages, which is fine. That said, I would like to similarly apply the weights to the values of landCO2 instead, so that the values across the 8192 grid cells for each of the 138 raster layers are appropriately transformed using the approach above. Evidently, this would be done prior to rasterToPoints being applied.

To do what I would like, I tried the following:

landCO2<-mask(RCP1pctCO2Median,b)
weightCO2 <- cos(landCO2[,"y"]*(pi/180)) #Notice that I skipped the "rasterToPoints" stage and   replaced CO2new with landCO2 to directly work with landCO2 for this

However, this results in the following error:

Error in landCO2[, "y"] : object of type 'S4' is not subsettable

landCO2 looks like this:

class       : RasterStack 
dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
resolution  : 2.8125, 2.789327  (x, y)
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,   
layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ... 
min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730, 0.53099146, 0.49757186,
0.45697752, 0.41382199, 0.46082401, 0.45516687, 0.51857087, 0.41005131, 0.45956529, 0.47497867, ... 
max values  :   96.30350,  104.08584,   88.92751,   97.49373,   89.57201,   90.58570,   86.67651,   
88.33519,   96.94720,  101.58247,   96.07792,   93.21948,   99.59785,   94.26218,   90.62138, ...

Why would the above error emerge? Is this also the correct approach to apply the weights to landCO2?

Here is from:

dput(landCO2)

datanotation = "FLT4S", byteorder = c(value = "little"), 
nodatavalue = -3.4e+38, NAchanged = FALSE, nbands = 138L, 
bandorder = c(value = "BIL"), offset = 0L, toptobottom = TRUE, 
blockrows = 0L, blockcols = 0L, driver = "raster", open = FALSE), 
data = new(".MultipleRasterData", values = structure(logical(0), .Dim = c(0L, 
0L)), offset = 0, gain = 1, inmemory = FALSE, fromdisk = TRUE, 
    nlayers = 138L, dropped = NULL, isfactor = FALSE, attributes = list(), 
    haveminmax = TRUE, min = c(3.51323445746594, 2.92551451275358, 
    3.76379186167912, 3.3067384979661, 3.34126201329632, 
    3.08386133801832, 3.01129632333484, 3.0647601880446, 
    3.48453623055105, 3.03053231936106, 3.19536433344149, 
    3.52855495856725, 2.90104612158757, 3.56892638125178, 
    3.47426237221953, 3.36702133980949, 3.48315423643726, 
    3.92361587221995, 3.90687036742463, 4.01543228703232, 
    3.41304727244007, 3.49068809707965, 3.2989754752411, 
    4.08699701583828, 2.85441317825034, 3.97032106120605, 
    3.684871323053, 3.44885465128755, 3.41068716698003, 3.71867905821301, 
    3.63124173755822, 3.36029075859188, 3.51076291417485, 
    3.60606616807582, 3.85469010434708, 3.62363212746417, 
    4.01465470183506, 3.90460196846042, 3.84035711667821, 
    4.47918876746572, 3.55786634099786, 3.93571036946891, 
    3.46479833997364, 4.15353306379117, 4.6771440494922, 
    3.76401468863884, 3.68025763686699, 3.98433894706614, 
    3.91096955447224, 3.43559741885861, 4.00445031471232, 
    4.48086544565188, 3.87482981677749, 4.08332665795915, 
    4.31136149250093, 3.70382263043114, 3.82149420265705, 
    4.0976117612541, 4.72182125382113, 3.84831596125847, 
    4.20844965732288, 3.79496924067602, 3.82120403737645, 
    3.96484623698029, 4.79580486673734, 4.03934290043253, 
    3.78843141554031, 4.04144783063136, 3.88475518193445, 
    3.77580790845968, 4.22916891094863, 3.37891057351953, 
    3.76609718198318, 3.66466763443896, 4.35687029283038, 
    4.15186590607846, 3.5213297389571, 3.81840347326943, 
    4.40712918335291, 4.49926237706677, 4.18081811171593, 
    4.23534950568547, 4.03825622270233, 4.19584640962629, 
    3.9942848069586, 3.60301350528971, 3.66142267849031, 
    4.41059042526361, 4.02642647615039, 3.96002336497629, 
    4.33060831293324, 4.56909277397865, 4.15728244422731, 
    4.95208776068967, 4.54139118465677, 4.25662558990673, 
    4.61052828421067, 4.15774228504737, 4.49014597435756, 
    4.48774327753287, 5.04747852185112, 4.41152644270915, 
    3.90896636477584, 3.87957397457568, 5.03593959365389, 
    4.27441673938052, 4.75848660621224, 4.3730213657783, 
    4.11031913867645, 5.53989297335251, 4.86124819051654, 
    4.11291051062372, 4.64780621813198, 4.67344226009189, 
    4.97821775717528, 4.87561446464774, 4.235786708869, 5.21463406723742, 
    4.91386599564275, 4.5040326461617, 4.39755890397611, 
    4.1122798656679, 5.62129607988027, 4.06041000097652, 
    5.78230509522252, 4.52243176940702, 4.33450278497182, 
    6.33154314118987, 4.84194052422235, 4.36239409992561, 
    5.06536910154253, 5.36775565871968, 4.58622692133259, 
    3.97402715566745, 4.38483710357093, 5.63010312542135, 
    5.19084796017202, 4.1911775173503), max = c(187.826931949416, 
    156.6457731009, 172.60496802628, 150.653561843192, 163.464290379554, 
    153.660276870287, 184.438549579016, 175.015738903312, 
    162.894831787447, 166.379381638683, 157.147247755231, 
    175.804302409214, 168.281236826442, 162.953071545666, 
    177.617207366887, 176.066608043038, 176.134624167252, 
    183.170432230683, 180.123080536513, 165.758984489366, 
    153.273621544173, 166.31783216726, 191.256207353269, 
    181.834164471366, 174.520181090338, 184.23551817623, 
    171.416536999497, 165.081042726422, 193.315459441979, 
    179.3304157447, 159.511279217247, 172.87064976477, 177.989077667051, 
    171.338298245905, 197.234144158504, 179.725063566265, 
    180.454647608295, 176.538985130977, 167.484739904838, 
    169.59549226423, 178.845183911962, 219.227785511576, 
    179.302375539991, 189.70609536145, 176.671747164801, 
    177.727204008252, 177.685635778916, 171.556028157277, 
    171.660911495565, 182.604148474644, 178.880778257735, 
    175.510428522397, 206.251924111661, 182.791122960718, 
    173.389614745975, 172.956023947333, 195.719172262432, 
    179.013389209727, 177.013185397468, 190.490072804889, 
    173.58294588578, 208.348677327084, 182.333424891621, 
    192.639120889362, 199.178892716937, 206.56841288783, 
    187.576048950698, 209.337462111751, 232.501635693939, 
    236.489306009982, 190.710977550501, 189.469236022643, 
    192.935159127228, 207.684714429639, 180.671009264941, 
    180.381004577515, 169.700639328308, 203.160405447161, 
    205.799407036166, 202.479083081546, 169.195944012218, 
    213.597443973355, 188.696902641166, 181.865464907605, 
    210.356417980395, 186.55099964459, 176.471593682072, 
    167.473492395052, 188.925691105779, 182.676410237466, 
    192.553436647052, 177.382951924105, 187.143067376759, 
    217.646413810112, 201.400167227448, 226.508857402951, 
    176.276089741395, 210.412824002824, 189.672038419686, 
    189.197220375707, 182.771473729409, 198.749990872168, 
    194.103429280842, 165.76628240291, 188.975499918648, 
    185.64837067388, 197.874000910364, 202.24484545222, 198.701451325542, 
    196.977756579872, 220.808595710649, 194.134417466921, 
    190.265171994542, 221.453181095372, 181.838130578399, 
    216.538949125262, 212.407913302794, 212.82087062114, 
    197.914417285938, 221.26754031219, 184.472694969736, 
    204.124747824096, 194.51746976059, 209.452072845306, 
    197.70029677602, 189.389852035546, 238.026551504564, 
    213.272616759493, 209.350743405636, 192.209439970673, 
    226.128907106621, 185.23781793192, 204.761643599548, 
    214.069989413456, 208.494819342333, 200.178394743342, 
    223.090127686495, 228.673297329609), unit = "", names = c("layer.1", 
    "layer.2", "layer.3", "layer.4", "layer.5", "layer.6", 
    "layer.7", "layer.8", "layer.9", "layer.10", "layer.11", 
    "layer.12", "layer.13", "layer.14", "layer.15", "layer.16", 
    "layer.17", "layer.18", "layer.19", "layer.20", "layer.21", 
    "layer.22", "layer.23", "layer.24", "layer.25", "layer.26", 
    "layer.27", "layer.28", "layer.29", "layer.30", "layer.31", 
    "layer.32", "layer.33", "layer.34", "layer.35", "layer.36", 
    "layer.37", "layer.38", "layer.39", "layer.40", "layer.41", 
    "layer.42", "layer.43", "layer.44", "layer.45", "layer.46", 
    "layer.47", "layer.48", "layer.49", "layer.50", "layer.51", 
    "layer.52", "layer.53", "layer.54", "layer.55", "layer.56", 
    "layer.57", "layer.58", "layer.59", "layer.60", "layer.61", 
    "layer.62", "layer.63", "layer.64", "layer.65", "layer.66", 
    "layer.67", "layer.68", "layer.69", "layer.70", "layer.71", 
    "layer.72", "layer.73", "layer.74", "layer.75", "layer.76", 
    "layer.77", "layer.78", "layer.79", "layer.80", "layer.81", 
    "layer.82", "layer.83", "layer.84", "layer.85", "layer.86", 
    "layer.87", "layer.88", "layer.89", "layer.90", "layer.91", 
    "layer.92", "layer.93", "layer.94", "layer.95", "layer.96", 
    "layer.97", "layer.98", "layer.99", "layer.100", "layer.101", 
    "layer.102", "layer.103", "layer.104", "layer.105", "layer.106", 
    "layer.107", "layer.108", "layer.109", "layer.110", "layer.111", 
    "layer.112", "layer.113", "layer.114", "layer.115", "layer.116", 
    "layer.117", "layer.118", "layer.119", "layer.120", "layer.121", 
    "layer.122", "layer.123", "layer.124", "layer.125", "layer.126", 
    "layer.127", "layer.128", "layer.129", "layer.130", "layer.131", 
    "layer.132", "layer.133", "layer.134", "layer.135", "layer.136", 
    "layer.137", "layer.138")), legend = new(".RasterLegend", 
    type = character(0), values = logical(0), color = logical(0), 
    names = logical(0), colortable = logical(0)), title = character(0), 
extent = new("Extent", xmin = -181.40625, xmax = 178.59375, 
    ymin = -89.258464857103, ymax = 89.258464857103), rotated = FALSE, 
rotation = new(".Rotation", geotrans = numeric(0), transfun = function () 
NULL), ncols = 128L, nrows = 64L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"), 
history = list(), z = list())

Thank you, and any assistance would be greatly appreciated!

Can you see any "y" in your landCO2? It's very likely that you still need rasterToPoints where there is a "y". — jay.sf
Can you create a weights raster with the same resolution and extent as your landCO2 stack and use raster math to multiply the weights by the RasterStack? Alternatively, have you considered reprojecting your RasterStack a projection that better matches your needs, maybe an equal area projection? I don't work with worldwide data, so don't have suggestions off hand. The second part of your question may be better suited for the GIS stack exchange site: gis.stackexchange.com — Brian Fisher
@jay.sf Unless because this is a Raster Stack, there is a y component (i.e. the stacking of raster layers creates a y component)? Otherwise, it is essentially composed of 8192 values per raster layer, so ideally, I would like each of these values to have the weights applied to them. — Android17
@BrianFisher I am not sure if it is possible to create another raster based off of the original landCO2 raster stack object. Maybe it is possible to re-project landCO2 to a suitable projection that takes into account the differences in area between the area and the poles by making them equal. Do you happen to know of a suitable projection to do this? Thanks — Android17
Using dput is a not a good way to provide example data. Rather do something like b <- brick(ncol=128, nrow=64, nl=2, xmn=-180, xmx=180, ymn=-90, ymx=90, crs="+proj=longlat +datum=WGS84"); values(b) <- c(1:ncell(b), ncell(b):1) — Robert Hijmans

Brian Fisher Brian Fisher · Accepted Answer · 2019-11-26T02:45:21

I think this approach should work.

library(raster)
r <- landCO2[[1]] # Extracts one layer of the raster brick to pull info from. Might not be necessary

## Create raster layer that has the latitude values for each cell. using method from <https://stackguides.com/questions/22848836/r-assigning-x-or-y-coordinate-to-cells-of-a-raster-to-perform-calculations>
lats <- raster(matrix(yFromCell(r, 1:length(r)), nrow = 64, byrow = TRUE),  #get latitude from existing raster
               xmn = -181.4062, xmx=  178.5938, ymn = -89.25846, ymx = 89.25846 ,  # Set the extents of the raster
               crs = crs(r))  # Using projection of existing raster
weights <- cos(lats*(pi/180))  #Feed this raster into your weighting formula.
plot(weights) # to see what we're looking at

This looks reasonable with weights around 1 near the equator and falling off sharply as you approach the poles.
You can now use raster algebra to apply the weights to your RasterBrick.\

CO2weighted <- landCO2 * weights # this should be a RasterBrick with 138 layers.

Computing weights from latitude to apply to raster stack values

1 Answers