Raster “Cover” function

September 11, 2010
By

(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

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

What about the differences?

summary(err)

Min.             1st Qu.   Median     Mean       3rd Qu.     Max.

-1.26200 -0.25080 -0.04142 -0.04890  0.15190  1.56000


To leave a comment for the author, please follow the link and comment on his blog: Steven Mosher's Blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags:

Comments are closed.