Applying weights using cosine of latitude to Raster Stack values

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view

Applying weights using cosine of latitude to Raster Stack values

R-sig-geo mailing list
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:
    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"
Why would the above error emerge? Is this also the correct approach to apply the weights to landCO2?

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

Thanks, and any help would be greatly appreciated!
        [[alternative HTML version deleted]]

R-sig-Geo mailing list
[hidden email]