Moshtemp 5.1

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

library("animation")

source("Global.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"))

dir.create(dir)

ani.options(nmax=nlayers(world),ani.width=1080,ani.height=540,

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

filename = "world.html",ani.dev=png,title="Global Anomaly",

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

ani.start()

aninmateWorldMap(world)

ani.stop()

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

source("Global.R")

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(!is.na(Anom))

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(!is.na(Anom)))

# 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

png("CoolStations/Trends.png")

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

dev.off()

# generate the Charts

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

series <- zoo(Anom[,Inv$Index[j]],order.by=TIMELINE)

plotStationSeries("CoolStations",series,Inv[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

writeKml(Inv,outfile="CoolStations/UrbanCool")

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