Reputation: 105
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
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
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