Home > Uncategorized > RasterBrick testing

RasterBrick testing

The R community is great. If you want to learn R the help lists are an invaluable source for getting solutions quickly. I’ll suggest you following the posting guides to get the most out of it. That means reading the manuals first, reading the help archives and when you do post a question, POST CODE. post a working ‘sample problem’

I use the list on a daily basis. Although I haven’t always followed the rules, I found that when you do, you will get help. The quality of the helpers is great. Anyway, i’ve posted some noob questions about raster and as a result the package has got a few more features today. They havent been released to Cran yet, but you can pick them up from R Forge. If you are new to R I will suggest that you wait, as you have to install packages by hand. It’s straight forward, so at your own risk you can do this:

install.packages("raster", repos="http://R-Forge.R-project.org")

load the library in the standard fashion [library(“raster”)] and you are good to go. ( see ‘bricktest.R’ )

The new capabilities on the brick, are a “cellStats” feature, an ‘as.logical’ feature and some other otions for taking the area, which will get us weights.

Lets recall what we want to do. For every layer we want to take the temps, apply a weight. the weight is simply the contribution that a cell’s area makes to the entire area sampled weighted by the percentage of Ocean in that cell. So Imagine you have just one temperature report in a month. 10C. Imagine that cell is 50% water. Imagine that this cell has an area of 100sqKm. the calculation goes like this:

Areaweight= (Area_of_ cell*percent_of_Ocean_in_Cell)/Total Sum of area of all cells with temperatures

Final temp= temp*Areaweight.

That’s pretty simple and our goal is to do that all with brick math. So conceptually I have this:

BrickOfcellWeights ( the combination you get by multiplying a weighted Area Brick by the ocean percent layer)

BrickOfTemps:

The answer is the product of these two. And then you SUM the result per month. cellStats(brick,sum)

From a 3D brick to a vector of temps. To get there, I’ll start with a simple test of the new capability. One feature at time,  giving the developers feedback as I go.  Ready. What we do is test out the new cellStat()

######### GET the land/ocean  Percent of LAND

LandOcean<-getLandOceanMask(cellsize=5)

##### switch it to percent of ocean

Ocean<-1-LandOcean

#  go get the temperatures and read into a brick

seafile<-HadSST2ncdf

SST<-brick(seafile,varname=’sst’)

# set a zoo type time index (year/month). for plotting You can use ts() is you perfer. I like zoo.

timeline<-as.yearmon(1850+seq(0,length=nlayers(SST))/12)

# NOW, in one line calculate the average temperature per Month ( per LAYER, each layer is a month)

# THIS IS NOT AREA WEIGHTED. This is new rasterbrick capability I am testing

results<-cellStats(SST,mean)

That’s it. Now you dont have to loop over layers to pull out the mean temp of each layer in to a vector that is nlayers long.  brilliant. Next, I’ll show you how we multiply bricks. In this case I will multiply a brick by a single layer. Later we will have to multiply a brick by a brick. But for now, I’ll just do a simple illustration of brick*layer. I will use the raster that represents the percent of water in each gridcell. This is a constant month to month.  The point of this is not to calculate a temperature (yet) but rather just to show that the code works ( the new version of raster). One thing I wanted ( have to read the manual) was a ‘subset’ or ‘window’ like command on the brick. Like so  newbrick<-subset(bigbrick, layer=1,nlayers=36). I’ll have to search for that. Anyway, on to the code: get the cell stats(mean) of the (SST*Ocean) and stuff it in a vector that is nlayers long.

resultsPct<-cellStats(SST*Ocean,mean)

And then the graph. remember this is just a test, but we will compare the results anyways just to make sure that the underlying math is right. Applying a percentage ocean mask should make only a tiny difference, becuase we are just adjusting the coastal stations by the percent of water in them.

At the drop

Categories: Uncategorized
  1. Patagon
    August 26, 2010 at 6:16 PM

    One thing I wanted ( have to read the manual) was a ‘subset’ or ‘window’ like command on the brick. Like so newbrick<-subset(bigbrick, layer=1,nlayers=36). I’ll have to search for that.

    May be?

    unwantedtlyrs = c(1:3,7)
    newbrick = dropLayer(bigbrick, unwantedtlyrs)

    • Steven Mosher
      August 27, 2010 at 12:03 AM

      Latest drop will have a bug fix. Not sure what the build cycle is. Also a new function subset() which is a nice addition.
      There will be a new area function as well which does the “weighting” that is needed in this type of application.
      and applyLayers… can’t wait.

      intersect, union, setdiff applied to bricks? (hmm maybe that is in there already)

  2. Steven Mosher
    August 26, 2010 at 8:34 PM

    Thanks, Patagon.

    I found that last night. testing

  3. Roj
    August 27, 2010 at 2:36 AM

    Not my lucky day. ‘install.packages(“raster”, repos=”http://R-Forge.R-project.org”)’ just installs 1.3-11 again. Going to R-Forge and trying to download the latest Windows .zip tries to download ‘https://r-forge.r-project.org/bin/windows/contrib/latest/raster_1.4-7.zip’, but that directory has just raster_1.3-11.zip in it. The source files and the Mac files are both version 1.4-7, and would download. Not sure I want to get into trying to do it from the source.

    Looks like they don’t do latest builds for Windows, everything in that ‘latest’ folder for Windows is dated 14-August-2010, for the others it’s 26-August-2010. Maybe the stuff about FusionForge is relevant. Just have to sit back and watch…

    • Steven Mosher
      August 27, 2010 at 2:54 AM

      I had trouble ( MAC) trying to do a “hand” install, so I went the API route and it worked. HOWEVER the pacakge manager still finds two versions, so if you hit that problem its just trial and error with the manager. library(“raster”) may just grab the old version if it was not installed properly so I just used the GUI and package manager to do it. Robert H, has been adding capability faster than I can test it. have a look at his code, he is very very good. I think he’s commiting to Cran in a week or so, so you can always pick it there when it comes. I believe he does a daily build.

  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: