Home > GHCN, trb > trb-0.05: Gridding by Latitude and Longitude

trb-0.05: Gridding by Latitude and Longitude

2010 July 14

Introduction

Little boxes on the hillside,
Little boxes made of ticky tacky,
Little boxes on the hillside,
Little boxes all the same.

Equal Area Latitudinal Zones

So, we know that there are many more records for the United States than for Mongolia, more records for North America than for Asia. How do we even out the difference in the distribution within a single latitude band? And one answer, of course is to grid the band and find the mean value of each grid before calculating a latitude average. In the following, I use the same number of grids for each latitude band. And since each band is equal area (in this version), so the grids are equal area.

Code

###############################################################################
# trb 0.05 : 2010 07 14 : Gridding by Latitude and Longitude
# trb 0.04 : 2010 07 13 : Proportional Latitudinal Bands
# trb 0.03 : 2010 07 12 : Weighted Latitudinal Bands
# trb 0.02 : 2010 07 11 : Simple Latitudinal Bands
# trb 0.01 : 2010 07 10 : GHCN Simple Mean Average
# https://rhinohide.wordpress.com
###############################################################################
ver <- "0.04"

###############################################################################
# Read Data Files
###############################################################################
# ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/v2.mean
# ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/v2.temperature.inv
###############################################################################

#------------------------------------------------------------------------------
# Temperature File (GHCN formatted monthly means)
#------------------------------------------------------------------------------

# mm for 'monthly means'
mmwids=c(11,1,4,5,5,5,5,5,5,5,5,5,5,5,5)
mmcols=c("id","dupe","year","jan","feb","mar","apr","may","jun",
                            "jul","aug","sep","oct","nov","dec")

# read the file
mm <- read.fwf("../data/v2.mean", widths=mmwids, col.names=mmcols,
                header=FALSE, na.strings="-9999")

#------------------------------------------------------------------------------
# Inventory File
#------------------------------------------------------------------------------

# there are five occurrances of the pound hash in v2.temperature.inv
# R scripts chokes on reading those, remove them first

invcols <- c("id","name","lat","lon","alt","elev","ru","pop","terr","arid",
             "coast","cdist","airport","adist","vegetation","abc","bright")
invwids <- c(11,32,6,8,5,5,1,5,2,2,2,2,1,2,16,1,5)
inv <- read.fwf("../data/v2.temperature.inv",widths=invwids,col.names=invcols)

###############################################################################
# Set up arrays for latitudinal banding
###############################################################################

hindex=new.env(hash=TRUE)
hlat=new.env(hash=TRUE)
hlon=new.env(hash=TRUE)
for (i in 1:nrow(inv)) {
    hindex[[as.character(inv[i,1])]] <- i
    hlat[[as.character(inv[i,1])]] <- inv[i,3]
    hlon[[as.character(inv[i,1])]] <- inv[i,4]
}

mmlat <- rep(NA,nrow(mm))
mmlon <- rep(NA,nrow(mm))

for (i in 1:nrow(mm)) {
     mmlat[i] <- hlat[[as.character(mm[i,1])]]
     mmlon[i] <- hlon[[as.character(mm[i,1])]]
}

# one array to rule them all and ...
mm1 <- cbind(mm,mmlat,mmlon, rowMeans(mm[,4-15],na.rm=T))

# for everything there is a season
mm <- NULL
inv <- NULL
mmlat <- NULL
mmlon <- NULL
hlat <- NULL
hlon <- NULL

# general formula for a set number of latitude bands of equal area
nbands <- 36                # roughly equivalent to a 5x5 grid
bands <- rep(NA,nbands)
b = 90                      # start with santa
bands[1]=b
for (i in 1:nbands) {
    a = asin(sin(b*pi/180)-2/nbands)
    b = a * 180 / pi
    bands[i+1] = b
}
bands

# 'bars' are the longitudinal edges of the grid
# calculate the bars long now so we can just lookup later

nbars <- nbands * 2        # this could be set to any number
bars <- matrix(NA,nbands+1,nbars+1)
bars[,1] = -180            # start in Howland Island and move east ;-)
for (i in 1:(nbands+1)) {
  for (j in 1:nbars) {
    bars[i,j+1] = 360*(j/nbars)-180
  }
}
bars[1,]
bars[nbands+1,]

###############################################################################
# Calculate Annual Mean Average
###############################################################################

#------------------------------------------------------------------------------
# set up the range of years to calculate
#------------------------------------------------------------------------------
#year0 <- min(mm1[,3])
year0 <- 1880              # 1880 is an arbitrary start date
#year0 <- 2000              # use a short range during developement
year1 <- max(mm1[,3])
range <- year1 - year0 + 1

#------------------------------------------------------------------------------
# initialize some vectors to hold annual results
#------------------------------------------------------------------------------
lyear <- rep(NA, range)
lmean <- matrix(NA,range,nbands)
lcount <- matrix(NA,range,nbands)
gmean <- array(NA, dim=c(range,nbands,nbars))
gcount <- array(NA, dim=c(range,nbands,nbars))

#------------------------------------------------------------------------------
# loop through the range of years
#------------------------------------------------------------------------------
for (i in 1:range) {

    year <- i + year0 -1

    #--------------------------------------------------------------------------
    # loop through the latitude bands
    #--------------------------------------------------------------------------

    for (j in 1:nbands) {

        #----------------------------------------------------------------------
        # loop through the longitude bars
        #----------------------------------------------------------------------

        for (k in 1:nbars) {

            # select the year of this loop
            # select a latitude band (descending order)
            # select a longitude bar (ascending order)
            # note that stns lying on a grid line are counted twice - ok by me
            idx=(mm1[,3]==year
                 & mm1[,16]<=bands[j] & mm1[,16]>=bands[j+1]
                 & mm1[,17]>=bars[j,k] & mm1[,17]<=bars[j,k+1] )
            mmlt=mm1[idx,]

            # calculate the annual means
            # divide by 10 since means are stored as 10ths of T

            # fractional area for this latitudinal band
            # fa <- abs((sin(bands[j]*pi/180) - sin(bands[j+1]*pi/180)))/2

            # don't need fractional area when every band has equal area
            fa <- 1

            # only use years with data for every month
            #gmean[i,j,k]= (mean(rowMeans(mmlt[,4:15]),na.rm=T)/10) * fa

            # use every year, even those with partial data
            gmean[i,j,k]=(mean(rowMeans(mmlt[,4:15],na.rm=T),na.rm=T)/10)*fa
            gcount[i,j,k]=nrow(mmlt)
        }

        # combine the grids into lat zones
        lmean[i,j]=mean(gmean[i,j,],na.rm=T)
        lcount[i,j]=sum(gcount[i,j,])
    }
    year[i]=year
}

#------------------------------------------------------------------------------
# calculate the global land temperature
#------------------------------------------------------------------------------
# since each band is equal area, take a simple mean
ameans = rowMeans(lmean,na.rm=T)
mean(ameans, na.rm=T) #

###############################################################################
# Display
###############################################################################
# plot some stuff
# do this nband times
for (i in 1:nbands)
{
    # had some issues when trying to plot vectors with no data
    if ( length(lmean[!is.na(lmean[,i]),i]) > 0 ) {
        img=paste("trb-",ver,"-zone-",i,".png", sep="")
        tit=paste("Latitudes",bands[i],"to",bands[i+1], sep=" ")
        png(img, width=600, height=150)
        par(mar=c(2, 2, 2, 1))
        plot(lmean[,i] ~ lyear, ylab="Temp (C)", xlab="", type="l",
        main=tit, cex.main=0.8)
  }
}

# to avoid the 0 for 2010, drop last value
img=paste("trb-",ver,"-zone-all.png", sep="")
subt=paste("https://rhinohide.wordpress.com (",ver,"-0.04)",sep="")
png(img, width=600, height=150)
par(mar=c(2, 2, 2, 1))
plot (ameans[1:(length(ameans)-1)] ~ lyear[1:(length(ameans)-1)],
    type="l", xlab="", ylab="Temp (C)",
    main="GHCN Gridded Averaged Temps", sub=subt)

code available at here

Results

I have my MS Update Manager to auto-reboot and it did so during the night. So no pretty pics at this time for the ‘results’. I’ll update them when they become available. Note that the code was run successfully on a much short time span (2000-2010) so I feel comfortable posting at this time.

BWHAHAHAHA! There seems to be a rounding error in the trig functions that happen when the number of bands becomes greater than about 31. I have a new run started.

Advertisements
  1. Steven Mosher
    2010 July 17 at 12:38 pm

    the trig in R killed me too. I hate that. along with bugs in regexp.

    That said, Your approach is interesting.

    do you have GSOD in an amenable format..

  1. No trackbacks yet.
Comments are closed.