Search Results

Keyword: ‘moshtemp’

Moshtemp 5.1

September 30, 2010 Leave a comment

Time for another dump of the entire package.

Get the zip file 5.1 in the box to the right. ( shortly). unzip and run the following files if you havent already:

1. downloadall.R

2. setup.R

If you are running for the first time, You’ll note I added diagnostics. Sometimes the files get corrupted on download and these files will help me figure out where your problems are

3. If you have a old folder of “functions” from 5.0, replace it.

4. LandOcean.R should be run again to output data for the animation script.

New in this release is Animation.R



if(!file.exists("LandOcean/Globe.Rdata"))stop(" run Land Ocean to generate Input")

world <- loadGlobe("LandOcean")

layerNames(world)<- as.character(TIMELINE)

dir <- "Animation"

if(file.exists(dir))stop(cat(dir, "exists create new dir"))



ani.type ="png",outdir=dir,imgdir="IMAGES",

filename = "world.html",,title="Global Anomaly",

description ="world animation",autobrowse=TRUE,interval=.1)




print(" to generate a stand alone Movie use Imagemagick on the PNG")

When you run “landocean.r” again it will output a Globe.Rdata object. Animation runs by loading this or ANY raster brick and then creating a html animation of the layers over time. Very simple. The parameters for the animation are set with “ani.options()” then, ani.start() is called, the the actual plots of each layer are drawn in animateWorldMap(). You pass that routine a brick and a color pallete. And everything else just happens. The folder also conatins all the source images as *.jpg. You can turn this into a move by using the tool Imagemagick or on the mac I just loaded the jpeg into iPhoto and exported as a movie. You get a 1320 frame  *.mov file with a 1 second frame rate and on playback you can control the speed. The package Animation actually allows you to “call”   saveMovie() but I wasnt able to get Imagemagick to compile on Osx 10.5. And they only have binaries for 10.6. Project for a rainy day.

Next you can run “coolUrban.R” that program looks like this


yearsrequired <- 90

minMonths     <- yearsrequired*12

if(!file.exists("CoolStations")) dir.create("CoolStations")

# load the anomaly data and the total inventory

Anom <- loadAnomalies()

Inv  <- loadInventory()

# now we count the months each station (column) has actual temperatures and not

# missing data

monthCount <- colSums(!

monthCount <- ifelse(monthCount >= minMonths,T,F)

# Above we create a mask of true/false TRUE if the station has enough months of data

# then below, we mask off those stations that fail the test.

Anom <- Anom[ ,monthCount]

# now that we have changed the number of stations in the Anom structure we have to

# intersect the two key data structures to get the same stations in each

Data <- intersect.InvAnomalies(Inv,Anom)

Inv <- Data$Inventory

Anom <- Data$Anomalies

# transpose Anomalies so that stations are in columns and time is in rows

# recall that "intersect" will transpose the Anomaly array.

Anom  <- t(Anom )

Inv <- cbind(Inv,Months = colSums(!

# and above we add a column to our Inventory that shows how many months each station has.

results <- vector(length=nrow(Inv))

# now we are going to loop through each station and generate a slope using glm()

for ( i in 1:nrow(Inv)){

x <- (glm(coredata(Anom[, i])~time(Anom[, i])))

results[i] <- x$coefficients[2]

print( i)


# store the slopes  as columns in Inv as decadal slopes *120

# write an index to keep track of which anomalie series is associated

Inv <- cbind(Inv,Slope=results*120,Index=1:nrow(Inv))

# Now, select those rows of the inventory where Slope is less than zero

Inv <- Inv[which(Inv$Slope < 0),]

# sort it ( just for output later)

o <- order(Inv$Slope,decreasing =FALSE)

Inv <- Inv[o,]

# write the inventory sorted by slopes from coolest to zero

write.table(Inv, file="CoolStations/CoolStation.inv",row.names=FALSE,quote=F,sep=" ")

# write the histogram


hist(as.numeric(Inv$Slope), breaks=20, main="Linear slopes of Long term Stations",xlab="Decadal Trend")

# generate the Charts

for ( j in 1:nrow(Inv)){

series <- zoo(Anom[,Inv$Index[j]],



# now subset the urban sites that are cooling and generate a KML

Inv <- Inv[which(Inv$Rural=="U"),]

# then write out another inventory with just the urban sites


write.table(Inv, file="CoolStations/CoolUrbanStation.inv",row.names=FALSE,quote=F,sep=" ")

Categories: Uncategorized

MoshTemp 4.1

September 6, 2010 1 comment

Replaces 4.0, minor changes made here and there. totally replaces 4.0 as the testing move forward.

see 4.1

Run: DownloadAll.

Run: DataSetUp.

SSTTest reruns the SST analsyis and compares at HADSST2

Land1900Present: reruns the old version of program with a bench against CRU land figures.

Next version will integrate more of raster to the analysis, then we throw out a bunch of old code and use raster for everything.

First raster bit is working

















Adding area weighting and land masking






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




Categories: Uncategorized


September 4, 2010 2 comments

Now that SST has been processed with raster we can return to the land processing and put that on the same footing. To get there I’ll rebuild from scratch making some minor adjustments along the way. As we go foward the amount of code we have to write should become less and less because we are letting raster do things for us.

As much as I hate to trash old code I wrote putting everything into raster will be the right thing to do. In the end the vast majority of the code we have will be basic bookkeeping code and project organizing code. So, I spent sometime looking at new R packages ( project template) and spent time trying to bullet proof the installation and the file handling. As a result, the code will now check for installed packages and install them only of they are not installed. The other problem that remains is the “uncompress” issue with .Z files. Following Ron Broberg, SteveMc has a nifty little solution for this problem. First, what’s the problem:

Many of the files we deal with are tar and/or .Z files. .Z is a old Unix compression standard. many folks have moved on from this. But there  are some who still use the format. Ron’s solution and Steve’s after him was to put the uncompress exe in on the system path and then invoke an R system() call, calling the uncompress with  the file as a parameter. Pretty slick!  I tried to use the R package “uncompress” but I get failures as did SteveMc. I prefer the R package approach. So, I’m trying to contact the maintainer to update his package. The source is available, so any c programmer should be able to figure it out. When Time permits I’ll look at his source and see if I can find the bug. It’s a unexpected EOF type error. I’m suspecting that he only tested on Linux since he’s a linux guy and that the tests Cran runs are not sufficient. SteveMc runs windows and I use MAC. Wonder if Ron Broberg had any luck on linux.

We will start by rebaselining. Moshtemp4 Zip will go up at the drop and in the box on this page.

Unzip it. Make it your working directory and run Downloadall.R That script will take a while. It will download all the assets you need. The files, masks, inventories, and CRU results which we will eventually calibrate against.  you will then have to manually uncompress v2.mean.Z. On the Mac thats just a click and on linux uncompressed should also be built it.

next run SetUpData.R.  That script will create

1. V2Mean.Rdata: this is a processed version of v2mean with all duplicates averaged

2. V2Anomalies.Rdata: this is giant matrix of all stations. StationID is in the column. Every row is a month of anomaly in C ( not tenths) from 1900 to 2009.

3. SST.Rdata: a brick of SST anomalies from 1900 to 2009

After SetUpData completes, you can test it by doing this:


You could even loop thru all 5072 stations and plot every one [ see the routine, plotAnomalySeries() ] You can also check out the calibration data:


[1] “Years”   “Annual”  “Monthly” “Percent”
then Y$Monthly or Y$Annual to check out the numbers
Finally, Robert has sped up some features in raster (stackApply) and he’s got some new code that will help me build a raster brick directly from my matrix of V2Anomalies. So in one line we will take that V2Anomaly object and ‘grid’ it, average stations within a grid, and build a brick where every layer is a month. I end up throwing away some of the final bits of the code I wrote. So, while he and I test that stuff have at Moshtemp4.
Categories: Uncategorized

MoshTemp 3.0

August 22, 2010 4 comments

A final release for 3.0

Some substantial changes from 1.0, and 2.0. Mainly in project directories and output objects. All objects for a project are now in a project object. Individal objects can always be written out.

Next up, a  FDM module. Just an exercise to do a quick comparision.  Then a study or two..

Downlad the zip. Run the script “Land1900Present”

It creates a directory for the project runs the analysis and writes out way too much data.

Categories: Uncategorized


August 22, 2010 Leave a comment

Well Almost there. There are a couple of Graphics bugs in the World Maps. They worked before, but I made a slight change, so I’ll figure that out.

A bunch of HOUSEKEEPING code got integrated: When I fix the bug and release here is what you will see:

And so forth. The files all output, graphics all draw ( cept the damn maps)

All Zipped up this weekend

Categories: Uncategorized

MoshTemp303 CELL-a-brate

August 20, 2010 Leave a comment

Whats this button do?

A quick morning update rebuilding what I did last night, checking it all out, building from scratch.  We are one step closer. The last step we took was to get our 3 major pieces of data coordinated. We needed to ‘intersect” the Land,Anomalies and the Inventory.  Now, we turn it into Grid cells. For today, I’ll show what it looks like in code, and then you have the magnifico chart above. BEWARE this is still not area WEIGHTED.
















It should be somewhat clear. Except for the referencing part where you read in or “load” RDATA. Basically, if you load(xyz.Rdata) you get an object called…..obj. so load(GHCNINV.RDATA); invobj<-obj means load this Rdata file. That file gets stuffed into an object called “obj” Probably need to do some more reading on this, but I will bury this detail in an API. so that this will work : Inventory<-loadOBJ(GHCNINV.RDATA)

Anyway, so the Anomalies get loaded, the inventory gets loaded, we load the LandMask. Then we grid the inventory. That’s adding celldata and area data. Then we intersect all three. Then we hand it over to as.GridCell() which takes the stations, groups them in cells and averages per cell. Lastly, we take the the COLUMN means ( columns are months) and we plot. The graph is now area averaged. Not weighted yet, that will take one step. gridWeight.


Categories: Uncategorized


August 20, 2010 23 comments

I’ll probably update the code dump this weekend unless somebody screams. For now I’ll just go through some of the code that I’ve just added in and retested.

Recall where we are. We have V2Anomalies. Rows of time, columns of stations. We have a gridded inventory. The two have to work together. Here is the issue. When we process v2mean we loose stations. We apply a time window and we have CAM criteria. So V2Anomalies has around 5100 stations. The inventory has 7280. So, we should trim the inventory. That’s easy. However, we have another issue. The land mask. Do we want to weight cells by the area of the cell OR weight by the land area. If the land mask is TRUE we want to weight by land area, This means we might have to trim both v2mean and the gridded inventory. Also, we want be able to handle situations where we take the 5100 stations in V2Anomalies and trim them. Like, only select the airport stations. There are bunch of solutions to this. None pretty. I can just assume that the user will do this trimming on his own. Its easy. BUT if you forget, you get bad answers. So I want to enforce a reconciliation.

R has no global variables. The scoping rules are such that a function can’t modify a parameter you pass to it. Well, you CAN with a trick. I will avoid that. Instead I handle the problem this way:




Basically, call a function that intersects these 3 pieces of information in one place and returns TWO values, held in a list. That’s a basic R paradigm. Lets look inside that code:


if(is.null(gridInv$Cell)) stop(“Pass a valid grid inventory”)


if(mask) stations<-gridInv$Id[which(gridInv$Land>0)]



if(length(common)==0) stop(“No Stations in intersection”)

InvOut<-gridInv[which(gridInv$Id %in% common),]

anomalyOut<-anomalies[,which(anomaly.stations %in% common)]



Starightforward. We check to make sure that the caller passed in a gridded inventory. Then we select stations that have a value for Land that is greater than ZERO, iff the landmask is true. This drops a few stations on small islands and ships. Then we ‘get’ the stations from the Anomalies. the function that does this just gets the column names. Then we do some set logic to find the intersection. common<-intersect(stations,anomaly.stations) Then we trim BOTH the anomalies and the inventory so they are justified with each other. If you dont do this ( trust me) you will have troubles. Lastly we set the values to two items in a list.

The process flow looks like this: read in the stored v2Anomalies from disk. Read in the entire Inventory. Do whatever you like to slice and dice the inventory. Like select airports, rural, high altitude, whatever. Build as many different inventory options as you like. Then intersect them with the BASELINE Anomalies.

The next piece of code takes us to the next level of geographical data.


if(is.null(gridinv$Cell)) stop( ” Pass a Gridded Inventory”)



if(!(nstations==anomaly.stations)) stop( “Inventory Mismatch”)






cellmatrix<-aggregate(anomalies, by=idx, mean, na.rm=T)






What is going on here? In a nutshell. We are averaging stations WITHIN a grid cell. First I do yet ANOTHER check on the issue of justification of the inventories. I cannot prevent the user from passing in the WRONG inventory, sizewise. Then we sort the inventory and the anomalies, then we apply the function aggregate. This collects all the stations in a cell and averages them. Note also that I transpose anomalies to get the rows as columns and columns as rows. We fix the NA issues that arise when there is no data. Lastly we take the first column ( which contains cell numbers) and make row names out of it. We now return a Cell Matrix that is Gridcells by row, and Time by column.

Thats it. A couple more routines and we are done.

Categories: Uncategorized

MoshTemp301: Location,Location,Location

August 20, 2010 8 comments

MoshTemp301 is all about Location. Putting stations on the map.









Just a short explanation here, and then more before I do the final drop of MT301. If you run MT301, you get the chart above and the output below

Id                    Name    Lat   Lon Cell      Area      Land
1   10160355000                  SKIKDA  36.93  6.95 2103  88477.57  27428.05
2   10160360000                  ANNABA  36.83  7.82 2103  88477.57  27428.05
3   10160390000            DAR-EL-BEIDA  36.72  3.25 2102  88477.57  23888.94
4   10160395001              FTNATIONAL  36.52  4.18 2102  88477.57  23888.94
5   10160400001              CAP CARBON  36.80  5.10 2102  88477.57  23888.94
6   10160402000                  BEJAIA  36.72  5.07 2102  88477.57  23888.94
7   10160403000                  GUELMA  36.47  7.47 2103  88477.57  27428.05
8   10160419000            CONSTANTINE   36.28  6.62 2103  88477.57  27428.05
9   10160425000                  CHLEF   36.22  1.33 2101  88477.57  13959.79
10  10160425001            ORLEANSVILLE  36.17  1.50 2101  88477.57  13959.79
11  10160430000                MILIANA   36.30  2.23 2101  88477.57  13959.79
12  10160444000          BORDJ BOU ARR   36.07  4.77 2102  88477.57  23888.94
13  10160445000                  SETIF   36.18  5.42 2102  88477.57  23888.94
14  10160457000          MOSTAGANEM VI   35.88  0.12 2221  91909.99  91195.14

To the code. The first function refreshGhcnInv() downloads the station inventory from NOAA and then cleans up the text fields in that file ( so they print a bit nicer..more  work here) and creates and writes a Rdata object to your working directory. Then we Globe<-getLandMask() which means we load the data in the land mask. Next we read in the station inventory object:invObj<-getGhcnInventory(). The next step is critical. We take that inventory and we feed it to a function which will add the extra geographical information we need. Basically the grid cell the station lives in, the AREA of that gridcell and the actual LAND area. The GHCN inventtory has data from surface ships and some tiny islands. They will register ZERO land.  There are not many and it doesnt matter how you treat these stations in the final analysis. I will drop them. But rest assured I tested it both ways before making that decision. Here is that function. Take the ghcnInv data and the Globe and return me a gridded inventory gridObj<-as.GriddedInventory(invObj,Globe). Recall that in MoshTemp202 I added some bits ( Parameters) to the output objects to make them more encapsulated self documenting objects. The DATA I want is in the $Data element. get that: gridInv<-gridObj$Data. There now I have a GHCN Inventory file that has been “gridded.” I’ve added geographical information to it that is KEYED to the gridding I selected for my land mask: Lets see how:

> gridObj$Parameter


[1] “BASELAND1900-2009”


[1] 3


[1] NA


[1] NA


[1] NA


[1] “v2.temperature.inv”


[1] “2010-08-19 15:29:00 PDT”

So, whenever I am reading and writing these objects to file I have to make sure I update anything that changes. In this case when I created the gridded inventory I gave the source data, the DATE of that source data and since I took information that depends upon the cell size, I updated that as well. Recall the cell number has no INHERENT meaning. If the world is gridded in 3 degrees, then cell number 1 will cover a specific 3 degree by 3 degree grid (-180E,90N to -177E,87N) If I did a 5 degree grid cell one would cover a different grid. There are 7200 total grids in a 3degree world.
Next I do a simple summary to quick check that I ported by old code correctly  stations<-table(gridInv$Cell) and I make a histogram that shows a lot of cells with one station only. hist(stations,breaks=40,labels=T). The table function gives me a frequency count of stations in every cell. If a cell has 40 stations, the count is 40. 5 stations, the count is 5. The chart for that is above

lastly I’ll output a portion of the inventory and check gridInv[1:100,c(1,2,3,4,17,18,19)] I have selected the first 100 records, and output the first 4 fields and the last three fields. Those fields get appended when you call “as.GriddedInv”. They are $Cell, $Area, $Landarea.

so the first cell has a land area that is LESS than the grid cell area. WTF? ah, that means one thing, noob.
sucker must be near the coast and only PART of the 3 degree cell is land: Check for yourself:
Id    Name   Lat  Lon Altitude Grid_Alt Rural Population Topography Vegetation Coastal Cell     Area     Land
1 10160355000  SKIKDA 36.93 6.95        7       18     U        107         HI         NA      CO 2103 88477.57 27428.05

YUP! $Coastal=CO Hit Google to check that thing( later we output google earth tours wo0 hoo):

More later. I guess my inbox is filling

Categories: Uncategorized

MoshTemp 202

August 19, 2010 Leave a comment

In this installation we will ad a few bits of code here and there checking our work as we go along. I will dump a new Zip tomorrow. First, lets hit the biggest change. I’ll ship the package with a couple more files, V2Anomalies.Rdata and V2mean.Rdata. All the work from here on out uses V2anomalies.

a quick note. Every object now comes with a self documenting feature












Author=”Steven Mosher”)

This gets appended to every file. Basically giving you tracibility.

I also added some code to the geography section: We read in a land mask:There are two things here of note: the CELL.SIZE which I fixed at 3 degrees and the base resolution which is 1 degree. I could use ,25 degree data. Look at the aggregate term. If I start with 1 degree data I can select 2 degree squares, 2.5, 3 degree data or 5 degree data. factor<-cellsize/basesize is a key line. Input .25 degree data and you can aggregate differently. “factor” controls how many grids i average to make my final resolution. Execute this and you get a Land mask of size ‘cellsize’. The data is then turned into a fraction and you have the fractional area covered by land 0 to 1




land<-read.table(file,sep=” “)

world<-raster(as.matrix(land),xmn=-180, xmx=180, ymn=-90, ymx=90,

crs=”+proj=longlat +datum=WGS84″)




The I defined a bunch of helper functions. basic stuff






















And I added some analysis code.





Here is Moshtemp202. Loads your Anomalies and does a few simple tests:












Code at the drop. later cause I’m beat

Categories: Uncategorized


August 18, 2010 Leave a comment

Time to get Organized!

In MoshTemp201 we are going to spend a bunch of time organizing our resources. While not the final form the  directory structure you see is the one we will be working with for a while. The code from now on will be dropped in Zip files. Or individual files that you can drop into the structure that we set up.

Moving forward we need to organize our code and data to keep from getting overwhelmed. When  unzip moshtemp 201 you should get a file system that looks something like the structure above. Lets start, with a list of what you get and what it does.

v2.mean.Z and v2.mean. I’ve included the compressed file of GHCN v2.mean and the unzipped version. and land_water_masks_xdeg (Folder). These files will be generated AFTER you run moshtemp201 successfully. We download and unzip the file for you. Inside that folder when your done you’ll find the land mask file. Its a 1 degree file of land coverage percentages.

InstallPackages.R. If you want to know what packages we will use here the file shows you. Get started early and install these packages. You can run the file, or do it using your package manager. Let me know if there are any missing dependencies. I do need to work on this aspect of the system build and haven’t dug into managing packages.







Moshtemp10X(Folder) All the prior code for you pleasure.

V2mean.Rdata and V2Anomalies.Rdata. These two files will be created when you moshtemp201. It takes a LONG time, but we only have to do it once. Sorry about that. From here on out we just use the V2Anomalies data. They form the basis of everything we do. These are temperature Normals, in the style of CRU. The only reason for regenerating this file is to recalculate your window (1900-2009) or your anomaly period. The file contains everything you need. I may update this, by making the object more self documenting, that is by encoding the camCriteria and file dates into a larger object. Later on that, need to test it.

Global.R This file is key to everything we write from here on out. It executes many scripts located in the functions folder. All the functions we write are going to be “organized” loosely into a different scripts. For now some of the scripts are just empty dummy files. The routines will be added as we go along. By the time we are done all the scripts will have more code. Plus if you want to add something ( Like graphics) you can see where it goes. I’ve tried to modularize things so that graphics are independent of the programs you will write. The scripts we write will all take on the same form: get the data; process it; output it; graph it. I like to do things that way, because I may come back to a study I did and decide to graph other things without running the whole thing over.  The various scripts are named as below:












The next bit is the Functions folder. which contains a bunch of scripts.

The scripts include, constants, filenames,graphics,datafunctions,download,inventory functions, analysis functions, scripts for Robjects or Rdata files and urls. As I said some of these files, like Url.R are fully developed. Url for example contains the urls of most of the data I have worked with on this project. There will be more. For Moshtemp201, there are a few key ones so we will walk through those first, leaving the others for your self enjoyment, confusion or whatever.

download.R is straightforward and requires little explanation







# add code for uncompress




In Moshtemp201, we will just download the land mask and unzip it. And you see what is coming as we have other files to download. Also, I still need to report a problem on uncompress to the maintainer.

Next RObjects



Again straightforward. As we go along there will be a few more .Rdata files that we save, so I collect all the names here. I hate having quoted strings in my code. Typos R Us.

Next Output functions. I will collect all data output in this place. Grunt code. In the final version there are also routines for outputing all .Rdata files as ascii, for data exchange with non R systems. You’ll find some grunt code that starts to lay the foundation for that, but here is the workhorse function:





if(!((o==f)&(o==d))) stop(“number of objects must match file names and directores”)

for(calls in 1:length(objects)){






Just a function that takes three lists: a list of objects, a list of directories and filenames. When we get to the end, if we need to write out 5 .Rdata files, then we just call this function and it takes care of the rest. It calls the saveOBJ function which checks your file extension and writes the object. You’ll also find a function that takes a “xyxx.Rdata” file name and switches it to “xyxx.dat” for ascii output. So if you save V2mean.Rdata as.ascii() then the system will switch the names and call the right file writing method.

That’s It.

Now, for Moshtemp201







#now write the Rdata objects





This will take a long time to run as we push up to R’s memory limit at least on the MAC. But thankfully we only have to do it once. we download the land mask. It is unzipped for us. Then we run the functions that we learned about in Moshtemp101-105. Finally, we test out our output function by writing the .Rdata files to our working directory. Should work. Lemme know if you have troubles.
Finally test this at the console:
> Land<-getPercentLand()
> Land
class       : RasterLayer
filename    :
nrow        : 60
ncol        : 120
ncell       : 7200
min value   : 0
max value   : 1
projection  : +proj=longlat +datum=WGS84
xmin        : -180
xmax        : 180
ymin        : -90
ymax        : 90
xres        : 3
yres        : 3
If you see that then you have successfully downloaded the land mask. Try this:

> cellFromXY(Land,c(0,0))
[1] 3661
> cellValues(Land,cellFromXY(Land,c(100,27)))
> Area<-area(Land)
> cellValues(Area,cellFromXY(Area,c(100,27)))
The first command gets the ‘cell’ number of the longtitude/latitude pair 0,0 The world is a big long vector of cells based on the size of the grid I’ve selected in the constants.R file.
The second command gets the actual value of the FRACTION of land in the specified lon/lat.
The third command and fourth command get you a New raster of Area (sqkm) and then reports the actual area for a cell. Slick. Steve hates trig, so we will NOT mess around with any code that calculates weights of cells based on sin and cos and earth radius and degrees to radian or any of that. The package will assign stations to grid cells so NO hand written binning problems, boundary issues, issues at the pole. None of that nonsense to clutter up our code. Raster to the rescue. Cleaner code, always good.
Code at the drop in a Zip. Download. Unzip. Make it your working directory! or suffer the pain of file not found.

Categories: Uncategorized