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
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:
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.
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
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
What about the differences?
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.26200 -0.25080 -0.04142 -0.04890 0.15190 1.56000