Home > Uncategorized > RghcnV3 version 1.1

RghcnV3 version 1.1

I’ve just uploaded version 1.1 of  the package RghcnV3 to Cran. I’ve made a few changes that should make it easier for some folks to use. First I removed the requirement for rgdal. At the present time “rgdal” is not required. On the MAC installing it can be a little trouble, but if you RTFM you can get it done. I’ll look for ways to provide capability without “rgdal” but it will involve reformating several spatial datasets and that’s not always a simple task. Further I like to be able to allow the user to build everything from the data sources themselves. That’s not always possible, especially since some datasets are behind registration walls. That too can be handled using “RCurl”, however “RCurl” has also been a hassle for me on the MAC. That hassle is entirely my own screw ups in having outdated curl libs all over the place. The other change I made was to change the R requirements from 2.13 back to 2.11. There was a bug in the demo code ( how’d that happen?) that has been fixed.

New features in 1.1

windowV3: a function to select windows of years from the V3 data

downloadSST, readSST: functions to download and read HADSST2 data.

intersectInvAnomalies: a minor change to the output of the function to keep anomalies as a zoo object.

rasterizeAnomaly: a function to “grid” anomalies into a 3D raster object.

some other tidbits here and there, more demo code as well.

The most important addition is the “rasterizeAnomaly” function. This function hides all the nastiness of projecting a zoo object of anomalies into a raster brick. I’ll show you the code and explain what the nasty bits are.

rasterizeAnomaly<-function(inventory,anomaly,land=GLOBE5){
  require("zoo")
  require("raster")
  # check the various calling parameters

  # check inventory as inventory data.frame
  #check anomaly as a zoo
  # stop if they dont match
  # proceeding otherwise
  if(!is.zoo(anomaly)){
    print(cat("anomaly is a ",class(anomaly),"\n" ))
    stop("anomaly must be a zoo object")
  }
  if(!is.data.frame(inventory)){
    print(cat("inventory is a ",class(inventory),"\n" ))
    stop("inventory must be a data.frame")
  }
  if(names(inventory)[1]!="Id"){
    print("inventory column 1 must be Id")
    stop(cat("column 1 name is ",names(inventory)[1],"\n"))
  }
  if(names(inventory)[2]!="Lat"){
    print("inventory column 2 must be Lat")
    stop(cat("column 2 name is ",names(inventory)[2],"\n"))
  }
  if(names(inventory)[3]!="Lon"){
    print("inventory column 3 must be Lon")
    stop(cat("column 3 name is ",names(inventory)[3],"\n"))
  }

  matching <- identical(getStations(inventory),getStations(anomaly))
  if(matching==FALSE){
    print(" mismatching stations between inventory and anomalies")
    stop("use the function intersectInvAnomalies to rectify")

  }
  if(!class(land)[1] == "RasterLayer") stop("land must be a raster object")
  resolution <- res(land)
  if(resolution[1]!=resolution[2])stop(" must use a land mask with square cells")

  ###########  ready to process
  #  get the lon lats of stations
  #  suck the time out of the zoo object
  #  transform the zoo object into a matrix with right orientation
  #  rasterize with mean
  #test <- rasterize(xy,r,field=vals,fun=mean)
  xy <- getLonLat(inventory)
  TIME <- time(anomaly)
  # core data wacks the row names
  data <- coredata(anomaly)
  dimnames(data)[[1]]<-TIME
  tdata<-t(data)
  gridded <- rasterize(xy,land,field=tdata,fun = mean,na.rm=TRUE)
  return(gridded)

}

As you can see we do some checking to insure that  the inputs are as required. Essentially, if you feed the out put of intersectInvAnomalies() to the rasterizeAnomaly() function, then everything will be in place.

Then we do the following.

we get the stations Lon/lat data from the inventory.

we extract the time from the anomaly zoo object. This is going to be used as layer names in the final raster brick. We have to extract it before doing any other operations on the zoo object because those operations will drop the data we need.

we pull the coredata from the zoo object. When you do this, you lose the time data. at this point data is a matrix

then we  put the time data back into the row names of the data matrix. Thus, we have a matrix where the row names are TIME and the column names are station Ids.  To use the rasterize() function in raster, this has to change. rasterize works like this:

You have a 2 column data structure of Lon/Lat.  Every Row is a new pair of points.  If you had 1000 stations, you would have 1000 rows of data. The matrix, on the other hand, has 1000 columns. A column per station. And the rows represent time. The solution is just to transpose the matrix t(). after that you have stations in rows and columns represent times. When this is fed to rasterize, the points that lie in the same grid square are going to have a function applied: mean. That is, if 3 stations lie with a grid square, the final value of that grid will be the mean of those three. After rasterizing you have a raster brick. The brick will have layers for every month.  The layerName of the layer will be the date in fractional years. 100 years = 1200 layers.  That brick is returned to the caller. You can then plot the layers of that brick, aggregate them into years or decades, take the mean of the layer, etc. More on that in the next release where will look at area averaging and working with SST data.  I’ll also add some other convenience functions for doing things in zoo series. I have requests in for different styles of summarizing years, and seasons. That’s all pretty straightforward in zoo, but getting fluent in zoo takes a little time. so I may just wrap some zoo processing into simple functions.  That code is done ( from long ago). I’d rather not do too much of that, my preference would be to extend zoo and extend raster as packages, but that is going to require some OOP work.  Maybe later in 2.0. The big issue there is one package is S3 OOP and the other is S4. That’s not an issue, just a learning curve.

The package should hit CRAN in a couple days ( source first, binaries later)

Categories: Uncategorized
  1. June 28, 2011 at 3:48 AM

    Hi,

    Perhaps the recent changes in raster are useful for you. Now the Raster objects include a slot named “z” (a list) to store information as a time index. This slot can be set with “setZ” and retrieved with “getZ”. Besides, there is a function named “zApply” to apply a function on a Raster object over subsets of the values of the “z” slot (as the aggregate.zoo).

    Last, I am working in a new package for enhanced visualization of Raster objects. Still in development but already usable: https://r-forge.r-project.org/projects/rastervis/

    Cheers,

    Oscar.

  2. Steven Mosher
    June 28, 2011 at 3:53 AM

    Thanks oscar, I was just looking at your project this very second on Rforge as I’m probably going to host there. I was talking to Robert yesterday about the setZ and getZ, something I’ve been waiting for since I started working on this. I’ll probably add some functions like “window” for the brick object to wrap dropLayers. Early on I used to get raster drops daily because I was testing out the new stuff robert was doing, but then I got distracted by looking at OOP so I’ve got to get back up to speed on the raster code base.

  3. June 28, 2011 at 2:04 PM

    You don’t need to wrap your calls to cat() in print(), and you don’t need cat() in the stop() calls. The various checks don’t follow the same format; was that intended?

    • Steven Mosher
      June 28, 2011 at 10:02 PM

      Thanks, It turned out that way merely because I had made them all prints to begin with
      and then realized that I needed cat.. It’ll get cleaned up.

  4. July 16, 2011 at 8:35 AM

    Hi Steven,
    Finally getting to work on it. Trouble with downloadMask() (tho I may not use it). The message was:
    px/publish/stevenmosher/land_percent_qd.asc Publish/land_percent_qd.asc’
    Error in download.file(url = url, dest, mode = “wb”) :
    cannot open URL ‘http://www.drivehq.com/file/df.aspx/publish/stevenmosher/land_percent_qd.asc Publish/land_percent_qd.asc’
    In addition: Warning message:
    In download.file(url = url, dest, mode = “wb”) :
    cannot open: HTTP status was ‘400 Bad Request’

    • Steven Mosher
      July 16, 2011 at 10:59 AM

      Thats weird. Its the second time that’s failed. I’m putting the mask into the package in the future. You should be able to go to the ftp and download by hand

    • Steven Mosher
      July 16, 2011 at 11:18 AM

      ha. I forgot to pay my ftp bill!

  5. Roger Bivand
    August 16, 2011 at 9:49 PM

    You install rgdal on intel OSX by setRepositories(ind=1:2); install.packages(“rgdal”), unless you need special drivers, or are on PPC. Simple, no?

    • steven mosher
      August 17, 2011 at 1:36 AM

      Hi roger glad you stopped by. In the past I did all my work on the MAC. maybe I should have been clearer:

      “First I removed the requirement for rgdal. At the present time “rgdal” is not required. On the MAC installing it can be a little trouble, but if you RTFM you can get it done.”

      On the PC rgdal is a breeze. On the MAC it was also a breeze once you helped me and pointed me at the KingChaos stuff. So what I am going to do is partition my packges so that the only one that needs to has a depends on Rgdal. That way I can help those mac people who need to go through the little deal with installing from source. So Rgdal is in the flow, I’ve just re arranged by packages so that only one will depend on it. Its just a cleaner organization.

      basically I’ll have 5 packages

      1. GHCNV3 : which will only deal with getting non geometrical data: works on all OS
      2. CHGN; same as above
      3. DailyGHCN: same as Above
      4. stats Processing tools: same as above plus raster
      5. GEO Specific metadata stuff : will depend on Rgdal and Rcurl. and I’ll help those mac
      people who want to use this package. That package is maybe a couple weeks away

      basically 4 will take inputs from 5. If people want to redo everything from source data
      then they need 5 and they absolutely need rgdal. So the metadata package uses Rgdal and raster to create data that can be used by 4. At first I was going to put everything in 1.
      But then people started asking for other stuff in 1 and It wasnt looking as clean as the organization above. which is better for me because I get to use rgdal and isolate support
      questions to only those people who want to use package 5

  6. January 23, 2014 at 11:03 PM

    Attractive section of content. I just stumbled upon your site and in
    accession capital to assert that I acquire actually enjoyed
    account your blog posts. Any way I’ll be subscribing to your
    augment and even I achievement you access consistently quickly.

  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: