rgeospatial

Calculating the cosine of latitude as weights for gridded data


I am currently trying to calculate "weighted" spatial annual global values for precipitation using the object "RCP1pctCO2Mean", which is a raster brick. This is due to the size of the area being represented differently towards the poles relative to lower-latitude regions. I need to do this for each of the 138 years (i.e. 138 layers) of global precipitation data that I have. The idea would be to somehow apply weights to each grid cell value for each year by using the cosine of its latitude (which means grid cells at the equator would have a weight of 1 (i.e. the cosine of 0 degrees is 1), and the poles would have a value of 1 (as the cosine of 90 is 1)). The idea would also then be to make a time series of these values, from Year 1 to Year 138, after all newly derived grid cell precipitation values are averaged for each year, effectively creating 138 (weighted) averages)).

The object "RCP1pctCO2Mean" 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, ... 

Here is what I have done so far:

newraster <- rasterToPoints(RCP1pctCO2Mean) #This creates a dataframe of the raster stack, but I am not sure if this is necessary?

I then proceeded to do assign weights like this:

weight <- cos(newraster*(pi/180))

However, this yields strange but identical precipitation values (i.e. all values are 0.97 to 0.99, which is odd) across each grid cell for each layer. I am not sure what I am doing incorrectly (if there is anything incorrect) - could it be that the "pi/180" is not necessary? Also, once this is done, how to revert to a raster stack with the new values?

I also saw another function, called "getWeights", but I am not sure how relevant this is. I tried this but received strange identical precipitation values across my grid cells:

weight <- getWeights(cos(newraster))

head(weight)

[1] 0.9998492 0.9998492 0.9998492 0.9998492 0.9998492 0.9998492

Would this approach apply the weights appropriately?

Thanks, and I would greatly appreciate any help with this!!!!


Solution

  • the biggest error in your code is that once you apply rasterToPoints the returned object is a matrix and not a new raster. The way to extract information from a matrix is different to that of a raster:

    ##I am calling this prec_ts, it is a matrix not a new raster
    prec_ts <- rasterToPoints(RCP1pctCO2Mean)
    
    ##calculate weighting based on latitude
    weight <- cos(prec_ts[,"y"]*(pi/180))
    
    ##plot weight to see that it is changing as expected, i.e. 1 at equator, near 0 at poles, cosine.
    plot(prec_ts[,"y"],weight)
    

    enter image description here

    ##now apply weighting to the precipitation columns
    prec_ts[,3:ncol(prec_ts)] = apply(prec_ts[,3:ncol(prec_ts)], 2, function(x) x * weight)
    
    ##calculate averages
    averages = colSums(prec_ts[,3:ncol(prec_ts)])/sum(weight)