根据纬度计算权重以应用于栅格堆栈值

Question

我目前正在尝试将权重应用于栅格堆栈对象的原始降水值，称为“ landCO2”。这些权重应考虑到赤道和两极之间的面积差异。但是，我不确定如何将其应用于栅格堆栈的现有值。想法是将权重应用于138个栅格图层中每个图层的原始值，然后最终使用wrld_simpl将这些值绘制在全局地图上。这是到目前为止所做的：

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)

但是，此方法获得了每层所有网格单元的平均值，有效地创建了138个平均值。也就是说，我想类似地将权重应用于landCO2的值，以便使用上述方法适当地转换138个栅格层中每个栅格层的8192个栅格像元的值。显然，这将在应用rasterToPoints之前完成。

要做我想做的，我尝试了以下方法：

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

但是，这导致以下错误：

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

landCO2看起来像这样：

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, ... 

为什么会出现上述错误？这也是将权重应用于土地二氧化碳的正确方法吗？

这里来自：

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())

谢谢您，任何帮助将不胜感激！

Answer 1

我认为这种方法应该有效。

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 

这看起来很合理，在赤道附近的权重约为1，并且在您接近两极时急剧下降。现在，您可以使用栅格代数将权重应用于RasterBrick。\

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

Answer 2

这里是如何设置示例数据

library(raster)
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)

要获得重量，您可以这样做

a <- area(b)