Home > Uncategorized > Step by Step through Land1990-Present

Step by Step through Land1990-Present

Ok, for those of you who have downloaded Moshtemp3. I will walk you through the code, line by line. Ready?

source(“Global.R”)

Every Moshtemp program starts with loading the “global” scripts. This script, call others scripts. The scripts basically have all the functions ( well not all, but the tested ones, there is WAY WAY more) necessary to run the analysis. The subscripts are divided into different categories, like graphics, and output. Helps me keep organized ( and suggest that somebody else do graphics ) Things to watch out for is getting the loading order of the scripts wrong, as some scripts refer to others. basically in the constants area.

Land<-getLandMask()

Get the land mask. This is a file that we downloaded before. Its a 1 degree BASE map. By setting the cell size you can load it as a 2,3,5,6 degree mask. The default is 3 degrees. When we integrate the ocean, this may have to change. I could even go to a finer resolution mask and build say a 2.5deg mask. The Land Structure is a ‘raster’ object. Very cool. Since raster is a s4 object I could even write my own subclass and write methods on that class.hella cool.

Inv<-loadGridInventory(getwd())

we load the Station inventory. Now I’ve prebuilt an inventory for you. using GHCN inventory as a base. Then, we loaded the mask and added data to the inventory. Every station has a lat/lon, metadata, AND it has a ‘cell’ number, “Area” and “land” variable. WHY calculate this stuff every time you run? so we have an inventory that is ‘gridded” station X has a cell number (in a 3degree world) and that cell has an area, and a land area. No trig in the code. I hate trig.  if you like you can point at a specific inventory by directory and name

Anom<-loadAnomalies(getwd())

Load up the 1961-1990 normals. These have been calculated for all the stations in GHCN. The data is 1900-2009. You can go rebuild this baseline if you like and change the CAM criteria or the base period, whatever. I’ve picked mine.

Data<-intersect.InvAnomalies(Inv,Anom)

Now you have to make sure that the inventory you loaded and the anomalies match. Let me give you an example. Say you started with the whole inventory and then wanted Only High altitude stations or Northern hemisphere. You might load the inventory ( the complete one) take a subset (high altitude) and then load your anomalies. Problem? The anomalies are for all stations. So you need to intersect the two of the structures and take ONLY those anomalies you want. or your Anomaly object may have fewer stations than your inventory. Either way you need to put them on common footing.

Anom<-Data$Anomalies

Inv<-Data$GridInventory

Assign the Values reyurned by the “intersect” function to the respective objects. The intersection function returns TWO objects. An inventory and a set of anomalies. They have been reconciled to contain the same stations.

CellMatrix<-as.GridCells(Anom,Inv)

The next step is to take those two elements and merge them into a “cell matrix” A cell matrix has COLUMNS of Time, and ROWS of cells. The cells are Grid cells. lat/lon cells. We use the raster package to assign stations to cells.. Station at X,Y is in cell 52. cell 52 is a 3×3 Cell. So in “as.Grid” we take all the stations, assign them to a cell and take the average of the anomalies within that cell. So we go from 5000 or so Station series, to about 1500 “cell series” Time series for each cell. MONTHLY. The cell matrix has some data which we will collect later with the cell stats function.

AreaAdjusted<-CellMatrix*gridWeight(CellMatrix,Land,LAND.MASK)

Now, we adjust the anomalies in each cell by that cell’s weight. This is area weighting. gridWeight looks at each month. It sums the total AREA covered by stations in that month. It then calculates the weight for each cell. So if you have 3 cells each with an area of 50, each gets a .33 weight. if you have one cell with an area of 100000, well its weight is 1. FOR THAT MONTH. obviously the higher the total area the better. So I track the total area sampled each month ( more fun with that down the road)

AreaAdjusted<-zoo(t(AreaAdjusted),order.by=timeline,frequency=12)

I should probably “caste” the return object as.zoo before returning. Todo list.

MonthlyTemp<-monthlyAnomaly(AreaAdjusted,timeline)

Basically, you just take the AreaAdusted series and do a columnSum! I put it in a function.

YearlyTemp<-yearlyAnomaly(MonthlyTemp,timeline)

Same thing here. You just agregrate the monthly data into years.

AnnualArea<-annualAreaAnomaly(CellMatrix,timeline)

If I want to do world maps, I do the line above. That takes a monthly Cell Matrix. and turns it into a yearly matrix. CELLS in rows, Time in columns. I use the row names and cell names to carry that data. so column number 1 has the “name”  “1900”  That way I dont have to carry around extra baggage. being able to name the array dimensions is kinda cool. Note these anomalies are not area weighted. There we are DONE with computing.

MonthStats<-monthlyStatistics(CellMatrix,Land)

Now, we collect stats from the objects. You can look at the ascii versions of these files to see what I output. If you like, just dont call this function.

CellStats<-getCellStats(CellMatrix,Land,Inv)

Now, I collect stats from the cells, again check the ascii output for questions. Its all labelled and formated so you can pull it into excell or some other systen

StationCount<-stationCount(Anom)

Count the stations per month. How many report?

There we are done collecting stats from our process. Load the data. Process. Collect stats. Done.

#####################################################

# Now we write out the  OBJECT files we want to to our project

Directories<-directorySetUp(getwd(),”Land1900-Present”)

If you like yur projects neat. Just set up a directory structure. You name the top level, and the function builds default subdirectories. See below

PD<-Directories$Project

A place for the high level data, input data needed to reconstruct the study.

GD<-Directories$Graphics

A place for graphics written to file

AD<-Directories$Ascii

Want to share data with others in an ascii form? a place for that

MD<-Directories$Maps

A place for global maps to be written, if you like. Maps, sits below graphics.

# directory

Now, we document what we did and WHEN we did it.

params<-InputParameters

params$Title<-“BASELAND1900-2009”

params$Cellsize<-xres(Land)

params$Years<-yearRange

params$Cam<-camCriteria

params$Landmask<-LAND.MASK

params$SourceData<-GhcnV2_Zipped

params$SourceDate<-lastModified(GhcnV2_Zipped)

params$AnomalyData<-V2ANOMALIES.RDATA

params$AnomalyDate<-lastModified(V2ANOMALIES.RDATA)

params$Stations<-NULL #need to deprecate this

params$Time<-Sys.time()

params$Author<-“Steven Mosher”

Now we create a list object of all the data and outputs. We will write this blob to disk.

## opps found a line of code that does nothing.

Now we save that as an .Rdata object which we can reload whenever we like and extract the objects

OBJOutList<-list(Anom,Inv,CellMatrix,MonthlyTemp,YearlyTemp,

AnnualArea,MonthStats,CellStats,StationCount)

saveOBJ(OBJOutList,PD,”Land1900Present.Rdata”)

Then we write out the ascii versions of the data.

writeStationCount(AD,filename=STATIONCOUNT.RDATA,StationCount)

writeParameters(PD,filename=PARAMETERS.RDATA,params)

writeStations(PD,filename=STATIONS.RDATA,Inv$Id)

writeMonthStats(AD,filename=MONTHSTATS.RDATA,MonthStats)

writeCellStats(AD,filename=CELLSTATS.RDATA,CellStats)

writeInventory(AD,filename=GRID3INV.RDATA,Inv)

writeAnomalies(AD,filename=V2ANOMALIES.RDATA,Anom)

writeMonthlyTemp(AD,filename=MONTHLYTEMP.RDATA,MonthlyTemp)

writeYearyTemp(AD,filename=YEARLYTEMP.RDATA,YearlyTemp)

writeCellAnnual(AD,filename=ANNUALAREA.RDATA,AnnualArea)

Lastly, we plot some graphics to disk, dividing by 10 to get to the correct units.

plotAnomalySeries(GD,”YearlyTemps”,YearlyTemp/10,Smooth=ipcc,Title=”Annual”)

plotAnomalySeries(GD,”MonthlyTemps”,MonthlyTemp/10,Smooth=tri48,Title=”Monthly”)

And finally we draw global maps for every year 1900-2009

drawGlobalArea(MD,pname=”World”,AnnualArea/10)

That’s it

###############UPDATE

This fix will come shortly or you can just make the change

ProjectOBJ<-list(Parameters=params,

Anomalies=Anom,

Inventory=Inv,

MonthlyTemp=MonthlyTemp,

YearlyTemp=YearlyTemp,

AnnualArea=AnnualArea,

MonthStats=MonthStats,

CellStats=CellStats,

StationCount=StationCount)

OBJOutList<-list(Anom,Inv,CellMatrix,MonthlyTemp,YearlyTemp,

AnnualArea,MonthStats,CellStats,StationCount)

saveOBJ(OBJOutList,PD,”Land1900Present.Rdata”)

saveOBJ(ProjectOBJ,PD,”Land1900Present.Rdata”)

Categories: Uncategorized
  1. No comments yet.
  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: