Android17
Android17

Reputation: 105

Computing weights from latitude to apply to raster stack values

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!

Upvotes: 0

Views: 127

Answers (1)

Brian Fisher
Brian Fisher

Reputation: 1367

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://stackoverflow.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 

enter image description here

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. 

Upvotes: 1

Related Questions