Raster “Cover” function

September 11, 2010

(This article was first published on Steven Mosher's Blog, and kindly contributed to R-bloggers)

The last bits of raster are coming together and everything is working so we are getting really close here. Today i will cover the raster function ‘cover’ which robert has just updated to handle a brick function. Lets review the situation with a toy example: First, we read in the Land data into a Brick. A brick is object with layers of rasters. One layer per month. lat*lon*time. A single raster layer will look like this. A matrix of NA and temps. The NA can be missing reports OR cells that are not Land

3 3 3

NA NA 3

NA NA NA

The SST layers in the ocean brick also follow the same rules. The SST Layer could look like this

NA NA 4

4 4 NA

4 4 4

So generally the SST layer will have values where the land does not and vice versa. EXCEPT for coastal cells. See the upper right hand cell? the land has 3 and the ocean has 4.

So we are going to combine two operations. First, we add the ocean to the land. This will give the following:

NA, NA, 7

NA NA NA

Then we will “cover” the Land with the ocean and the coast. That will copy filled cells to NA cells. If the land has a value and the ocean has NA, then the value is copied. In the end we have all three components in one layer. It look like this in code:

Coast<-Land+Ocean

Globe<-cover(Coast,Land,Ocean)

The three bricks of coast, land and ocean are all combined with the cover function. In the end the code looks like this. Please note I have to integrate “pointsToRaster” yet.

source(“Global.R”)

Land<-getMask(file=LandWaterMask_File) # get the mask file of land percentage 0-1

OceanPercentMask <- 1 – Land # the percent of Ocean is 1-land

Inv<-getGhcnInventory() # get the stations ID

Inv<-as.GriddedInventory(Inv,Land) # This line goes away

Anom<-loadAnomalies(getwd()) # Load the anomalies for all stations

Data<-intersect.InvAnomalies(Inv,Anom) # rectify the inventory list with the stations in the anomaly period

Anom<-Data$Anomalies # the anomalies for the stations

Inv<-Data$GridInventory # the stations

CellMatrix<-as.GridCells(Anom,Inv) # MY version of points to raster. this goes away

m<-matrix(nrow=ncell(Land),ncol=ncol(CellMatrix))

b<-brick(Land) # declare a brick with a layer of format Land

clearValues(b) # clear the decks

m[as.numeric(row.names(CellMatrix)),]<-as.matrix(CellMatrix) # build a mtrix of temps

b<-setValues(b,m) # create a brick of land temps

Weights <- area(b,na.rm=T,weight=T) # area weighted temps

Temps<-Weights*b*Land # correct for land coverage percent

SST <- loadSST(getwd()) # get SST anomalies

Weights <- area(SST,na.rm=T,weight=T) # weight the cells by area

Ocean <- SST*Weights*OceanPercentMask # weight by actual area of water in the cell

Coastal <-Ocean+Temps # add the landcoast to the seacoast.

all = cover(Coastal, Ocean, Temps) # mask NAs and combine all three

x<-cellStats(all,sum) # Monthly temps

Below find HadCrut in red, Moshtemp in black. Note that differences are due to slightly different data between HadCrut and moshtemp. We use GHCN raw. And yes, I show more warming in recent months