Home > ERSST, GHCN, trb > trb-0.09: Planet Ocean

trb-0.09: Planet Ocean

2010 July 25

Introduction

https://i2.wp.com/www.livingoceanproductions.com/Living_Ocean_Productions/Ocean_Conservation_files/Wave%20for%20Ocean%20Cons%20Page.jpg
Source: http://www.livingoceanproductions.com/Living_Ocean_Productions/Ocean_Conservation.html

How inappropriate to call this planet Earth when it is quite clearly Ocean.
-Arthur C. Clarke

Notes

I pretty much just duplicated sections of the previous land-based worked to run ocean data in parallel, leaving the subroutine refactor to reduce redundant code for another time.

The ocean data is ERSSTv3b data that I processed in the end of May. Discovered that I had incorrectly assigned lat/lons to the data, shifting everything 180 degrees to the West. Data visualization is a wonderful thing!

For overlaps, I currently just take the mean of the land and sea temp data.

I’m really enjoying working with the data as raster files. But it reminds that there is a whole second paradigm in GIS work – vector data.

Code

###############################################################################
# trb 0.09 : 2010 07 24 : Planet Ocean
# trb 0.08 : 2010 07 23 : Raster / Pretty Pictures / Fractional Area Weights
# trb 0.07 : 2010 07 18 : Performance Refactor (22 seconds with data read)
# trb 0.06 : 2010 07 17 : Performance Refactor (3x faster, still 8hrs for all)
# 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.09"

read <- FALSE # set to FALSE if this is your first time through
read <- TRUE # set to TRUE if you have previous data to read

#flag_latarea <- TRUE # set to TRUE if you want equal area lat bands
flag_latarea <- FALSE # set to FALSE if you want equal step lat bands

flag_land <- FALSE # set to FALSE to exclude land data
flag_land <- TRUE # set to TRUE to include land data

flag_ocean <- FALSE # set to FALSE to exclude ocean data
flag_ocean <- TRUE # set to TRUE to include ocean data

###############################################################################
# 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
# ??? url for ersst
###############################################################################

# this flag is to indicate if we must read for the first time
# or if we can read the data from previous runs

# TO-DO split into land read and ocean read ?
if (!read) {

  if (flag_land) {

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

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

    # read the file
    lmm <- read.fwf("../data/v2.mean", widths=lmmwids, col.names=lmmcols,
                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

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

  if (flag_ocean) {
  #----------------------------------------------------------------------------
  # Ocean Temperature File (GHCN formatted monthly means)
  #----------------------------------------------------------------------------
  # I'm cheating a little here, will need to handle import of native format
  # but for now, I am using the ERSST data I formatted as GHCN earlier

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

    # read the file
    omm <- read.fwf("../data/my.ersst.mean.all",widths=ommwids,
                    col.names=ommcols, header=FALSE, na.strings="-9999")

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

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

  } # end read ocean

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

# TO-DO make into subroutine
  hindex=new.env(hash=TRUE)
  hlat=new.env(hash=TRUE)
  hlon=new.env(hash=TRUE)
  if (flag_land) {
    for (i in 1:nrow(linv)) {
        hindex[[as.character(linv[i,1])]] <- i
        hlat[[as.character(linv[i,1])]] <- linv[i,3]
        hlon[[as.character(linv[i,1])]] <- linv[i,4]
    }

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

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

    # add inv lat lon
    lmm1 <- cbind(lmm,mmlat,mmlon)
 }

 if (flag_ocean) {
    for (i in 1:nrow(oinv)) {
        hindex[[as.character(oinv[i,1])]] <- i
        hlat[[as.character(oinv[i,1])]] <- oinv[i,3]
        hlon[[as.character(oinv[i,1])]] <- oinv[i,4]
    }

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

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

    # add inv lat lon
    omm1 <- cbind(omm,mmlat,mmlon)

  }

# ---------------------------------
# grid geometry
# ---------------------------------
# general formula for a set number of latitude bands of equal area
  nbands <- 36
  bands <- rep(NA,nbands)
  b = 90                      # start with santa
  bands[1]=b
  bands[nbands+1]=-90         # rounding errors at nbands ~ 32
  for (i in 1:(nbands-1)) {
    if (flag_latarea) {
      # equal area lat bands
      a = asin(sin(b*pi/180)-2/nbands)
      b = a * 180 / pi
      bands[i+1] = b
    } else {
      # equal lat step size
      bands[i+1] = 90 - (180/nbands) * i
    }
  }
  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
    }
  }

## weights for each cell
  weights <- matrix(NA,nbands,nbars)
  for (i in 1:(nbands)) {
    for (j in 1:nbars) {
     # assumes equal cells within a lat band
	 weights[i,j]<-abs((sin(bands[i]*pi/180)-sin(bands[i+1]*pi/180)))/(2*nbars)
    }
  }
  sum(weights[,]) # should equal 1

#------------------------------------------------------------------------------
# assign each record to a grid point
#------------------------------------------------------------------------------
# lat, lon, year, stn annual mean
# TO-DO: make into subroutine

    bar1 <- bars[1,]
    fm1 <- function(i) { sum(bands>i) }
    # sensitive to end points (-180,180)
    #fm2 <- function(j) { floor( (j + 180) * nbars / 360 ) +1 }
    fm2 <- function(j) { ceiling( (j + 180) * nbars / 360 ) }

    mmlat <- lapply(omm1[,16],fm1)
    mmlon <- lapply(omm1[,17],fm2)
    mmlon <- as.numeric(mmlon)

  if (flag_land) {
    mmlat <- rep(NA,nrow(lmm1))
    mmlon <- rep(NA,nrow(lmm1))

    mmlat <- as.integer(lapply(lmm1[,16],fm1))
    mmlon <- as.integer(lapply(lmm1[,17],fm2))

    lmm2 <- cbind(mmlat,mmlon,lmm1[,3],rowMeans(lmm1[,4:15],na.rm=T)/10)
    colnames(lmm2) <- c("ilat","ilon","yr","mean")
  }

  if (flag_ocean) {
    mmlat <- rep(NA,nrow(omm1))
    mmlon <- rep(NA,nrow(omm1))

    mmlat <- as.integer(lapply(omm1[,16],fm1))
    mmlon <- as.integer(lapply(omm1[,17],fm2))

    # my ersst data is in 100ths C not 10ths C
    omm2 <- cbind(mmlat,mmlon,omm1[,3],rowMeans(omm1[,4:15],na.rm=T)/100)
    colnames(omm2) <- c("ilat","ilon","yr","mean")
  }

  #------------------------------------------------------------------------------
  # set up the range of years to calculate
  #------------------------------------------------------------------------------
  if ( flag_ocean && flag_land ) {
    # find common overlap
    year0 <- max(min(lmm1[,3]),min(omm1[,3]))
    year1 <- min(max(lmm1[,3]),max(omm1[,3]))
  } else if (flag_land) {
    year0 <- min(lmm1[,3])
    year1 <- max(lmm1[,3])
  } else if (flag_ocean) {
    year0 <- min(omm1[,3])
    year1 <- max(omm1[,3])
  }

  # or override the above
  year0 <- 1880
  year1 <- 2009

  # calculate the range of years
  range <- year1 - year0 + 1

#------------------------------------------------------------------------------
# trim the matrix to the span of years
#------------------------------------------------------------------------------
# TO-DO make into subroutine
  if (flag_land) {
    idx <- lmm2[,3]>= year0 & lmm2[,3]<=year1
    lmm2 <- lmm2[idx,]
    lmm2[1,]
  }

  if (flag_ocean) {
    idx <- omm2[,3]>= year0 & omm2[,3]<=year1
    omm2 <- omm2[idx,]
    omm2[1,]
  }

#------------------------------------------------------------------------------
# clean up
#------------------------------------------------------------------------------
  hlat <- NULL
  hlon <- NULL

  if (flag_land) {
    lmm <- NULL
    lmm1 <- NULL
    lmmlat <- NULL
    lmmlon <- NULL
    linv <- NULL
  }

  if (flag_ocean) {
    omm <- NULL
    omm1 <- NULL
    ommlat <- NULL
    ommlon <- NULL
    oinv <- NULL
  }

#------------------------------------------------------------------------------
# write the data
#------------------------------------------------------------------------------

  if (flag_land) { write.table(lmm2,".Rlmm2",row.names=F,col.names=T) }
  if (flag_ocean) { write.table(omm2,".Romm2",row.names=F,col.names=T) }
  write.table(bands,".Rbands",row.names=F,col.names=F)
  write.table(bars,".Rbars",row.names=F,col.names=F)
  write.table(weights,".Rweights",row.names=F,col.names=F)
  write.table(c(year0,year1),".Ryears",row.names=F,col.names=F)

} else {

#------------------------------------------------------------------------------
# read the data
#------------------------------------------------------------------------------

  if (flag_land) { lmm2 <- read.table(".Rlmm2",header=T) }
  if (flag_ocean) { omm2 <- read.table(".Romm2",header=T) }

  bands <- read.table(".Rbands")[,1]
  nbands <- length(bands) -1
  bars <- read.table(".Rbars")
  nbars <- ncol(bars) -1
  weights <- read.table(".Rweights")
  yrs <- read.table(".Ryears")[,1]
  year0 <- yrs[1]
  year1 <- yrs[2]
  range <- year1 - year0 + 1
} # end if for data reads

# this is a trick to ensure that some value (NA) exists for every lat band
# required to ensure that the tapply mean factors every lat band
a1 <- matrix(NA,nbands,4)
a1[,1] <- rep(1:nbands)
a1[,2] <- 1
a1[,3] <- year0
colnames(a1) <- c("ilat", "ilon", "yr", "mean")
if (flag_land) { lmm2 <- rbind(lmm2,a1) }
if (flag_ocean) { omm2 <- rbind(omm2,a1) }
a1 <- NULL

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

# convenient list of years for later use
lyear <- seq(year0:year1) + year0 -1

#------------------------------------------------------------------------------
# grid the temperature data
#------------------------------------------------------------------------------
if (flag_land) {
  lgmean <- array(NA, dim=c(nbands,nbars,range)) # holds the mean of each cell
  lgmean <- tapply(lmm2[,4], list(lmm2[,1],lmm2[,2],lmm2[,3]), mean, na.rm=T)
}

if (flag_ocean) {
  ogmean <- array(NA, dim=c(nbands,nbars,range))  # holds the mean of each cell
  ogmean <- tapply(omm2[,4], list(omm2[,1],omm2[,2],omm2[,3]), mean, na.rm=T)
}

#------------------------------------------------------------------------------
# calculate the global land temperature
#------------------------------------------------------------------------------
# Saving these values until I calculate them from the raster map
#if (flag_land) { mean(lamean, na.rm=T) } # 12.00907
#if (flag_ocean) { mean(oamean, na.rm=T) } # 13.13318

###############################################################################
# Raster Display
###############################################################################

library(sp)
library(rgdal)
library(raster)
library(maps)

# take 2 rasters and return the mean
mean_raster_values <- function(r1,r2) {
  # TO-DO error check to ensure both are same size
  rr <- matrix(nrow=nrow(r1), ncol=ncol(r1))
  for (r in 1:nrow(rr)) {
   rr[r,] <- colMeans(rbind(getValues(r1,r),getValues(r2,r)), na.rm=T)
  }
  return(rr)
}

# take 2 rasters and return the mean
plot_raster_values <- function(r1,title,subt) {
    plot(r1,col=rev(rainbow(256)), main=title, sub=subt)
    grid(nx=nbars,ny=nbands,col="lightgray")
    grid(nx=nbars/2,ny=nbands/2,col="darkgray",lty="solid")
    map("world",add=T)
}

if ( flag_land ) {

    lr <- raster(nrows=nbands,ncols=nbars,
            xmn=bars[1,1],xmx=bars[1,nbars+1],
            ymn=bands[nbands+1],ymx=bands[1],
            crs="+proj=longlat +datum=WGS84")
    # set a year to display
    y <-130
    lr <- setValues(lr,as.vector(t(lgmean[,,y])),layer=1)

    img=paste("trb-",ver,"-land-grid-",lyear[y],".png", sep="")
    tit=paste("Global Gridded Mean Temperature (",lyear[y],
              ")\n Land Data: GHCN unadjusted",sep="")
    subt=paste("https://rhinohide.wordpress.com (trb ver ",ver,")",sep="")
    png(img, width=600, height=350)

    plot_raster_values(lr,tit,subt)
}

if ( flag_ocean ) {

    or <- raster(nrows=nbands,ncols=nbars,
            xmn=bars[1,1],xmx=bars[1,nbars+1],
            ymn=bands[nbands+1],ymx=bands[1],
            crs="+proj=longlat +datum=WGS84")
    # set a year to display
    y <-130
    or <- setValues(or,as.vector(t(ogmean[,,y])),layer=1)

    img=paste("trb-",ver,"-ocean-grid-",lyear[y],".png", sep="")
    tit=paste("Global Gridded Mean Temperature (",lyear[y],
              ")\n Ocean Data: ERSSTv3",sep="")
    subt=paste("https://rhinohide.wordpress.com (trb ver ",ver,")",sep="")
    png(img, width=600, height=350)

    plot_raster_values(or,tit,subt)
}

# plot both land and ocean
if ( flag_land && flag_ocean ) {

    lor <- raster(nrows=nbands,ncols=nbars,
            xmn=bars[1,1],xmx=bars[1,nbars+1],
            ymn=bands[nbands+1],ymx=bands[1],
            crs="+proj=longlat +datum=WGS84")
    # set a year to display
    y <-130
    lor <- setValues(lor,as.vector(t(mean_raster_values(lr,or))),layer=1)
    img=paste("trb-",ver,"-land-ocean-grid-",lyear[y],".png", sep="")
    tit=paste("Global Gridded Mean Temperature (",lyear[y],
              ")\n Land Data: GHCN unadjusted / Ocean Data: ERSSTv3",sep="")
    subt=paste("https://rhinohide.wordpress.com (trb ver ",ver,")",sep="")
    png(img, width=600, height=350)

    plot_raster_values(lor,tit,subt)
}

# plot land and ocean overlap
if ( flag_land && flag_ocean ) {

    mask1 <- ! is.na(getValues(lr))
    mask2 <- ! is.na(getValues(or))
    mask <- mask1 & mask2
    aav <- as.vector(getValues(lor)) * mask
    aav[!mask] <- NA
    lor2 <- raster(nrows=nbands,ncols=nbars,
            xmn=bars[1,1],xmx=bars[1,nbars+1],
            ymn=bands[nbands+1],ymx=bands[1],
            crs="+proj=longlat +datum=WGS84")
    lor2 <- setValues(lor,aav,layer=1)
    img=paste("trb-",ver,"-land-ocean-grid-overlap",lyear[y],".png", sep="")
    tit=paste("Global Gridded Mean Temperature Overlap (",lyear[y],
              ")\n Land Data: GHCN unadjusted / Ocean Data: ERSSTv3",sep="")
    subt=paste("https://rhinohide.wordpress.com (trb ver ",ver,")",sep="")
    png(img, width=600, height=350)

    plot_raster_values(lor2,tit,subt)

}

trb-0.09.R

Results

trb-0.09-land-grid-2009.png

trb-0.09-ocean-grid-2009.png

trb-0.09-land-ocean-grid-2009.png

trb-0.09-land-ocean-grid-overlap-2009.png

Advertisements
  1. 2010 July 26 at 11:59 am

    How are you dealing with cells that have both land and ocean stations? Is the result the land-fraction weighted average of the two?

  2. 2010 July 26 at 8:09 pm

    Nope. In this version, it is a simple average of the two: (land + ocean)/2
    To get to fractional weighting for the land -v- ocean in a cell, I will need to bring in a coastal raster map and compute the fraction from it. Pretty simple conceptionally, just haven’t got there yet.

  3. steven Mosher
    2010 August 4 at 7:44 am

    Raster is very cool

    looking forward to your native import of the SST data ( lm lazy)

  1. No trackbacks yet.
Comments are closed.