Monday, 5 December 2011

World Ocean Atlas

How to open data from World Ocean Atlas in R
 
####################################################
####################################################
# OPEN data from World Atlas NOAA 
# http://www.nodc.noaa.gov/cgi-bin/OC5/WOA05/woa05.pl
####################################################
### Tim Bover from WOA09 (how to obtain the file from internet, and download to your computer):
#   - Go to http://www.nodc.noaa.gov/OC5/SELECT/woaselect/woaselect.html
#   - select the oceanographic variable you are interested in
#   - select geographic area (or leave whole world), type of data
#    (probably climatological mean in your case), time period and depth
#    (depth is not that important.  The data files contain all standard
#     depths available)
#   - click on 'Show Figure'
#   - on the resultant page, on the left, above the figure, you will see
#    'Data Access: ASCII ArcGIS'
#   - click on the ASCII link (or left mouse button and 'save target as').
# readLines(file.choose()) # Let the user choose the local file to read (thank's to Luke Miller)
# setwd("C:\\Users\\Marta Rufino\\world")
h(readLines("data\\WOA05\\temperature.gz")) #this may take some time...
 
# note the path indicates the folder wher you have the file. 
# You can get the path, by clicking the left mouse button and properties, 
# using windows explorer, for example or changing R working directory to the folder where you have your file
# Note also that the file name will be different from the current one (it is generated automatly by WOA)
# Whatever way you have to obtain the path, note that you need to change the bars: instead of \ put \\ or /
T=read.table("data\\wfig1323085620.23584.csv.gz",
     sep=",", col.names=c("lat","lon", 
     paste("t",c(0,10,20,30,50,75,100,125,150,200,250,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500,1750,2000,2500,3000,3500,4000,4500,5000,5500), sep="")))
dim(T)# to check if everything is ok...
summary(T)
 
## If we want to work with the data as a spatial object:
require(sp)# this package is required for working on most spatial objects.
T.spa=T
coordinates(T.spa)=~lon+lat
gridded(T.spa)=TRUE
#proj4string(T.spa)=CRS("+proj=longlat +datum=WGS84") # if we wanted to define the projection... which is not exactly this one for all but 
summary(T.spa)
 
#image(T.spa, "t0") #for plotting temperature, at zero depth, for example. Note we can plot the 35 depth levels.
spplot(T.spa, c("t0","t10"), col.regions=rev(heat.colors(20)))
#spplot(T.spa, c("t0","t10"), col.regions=rev(rainbow(100)), at=seq(30,40,l=100),main="Temperature")
 
Created by Pretty R at inside-R.org


This map represents temperature at 0 and 10 m depth ;)

Tuesday, 29 November 2011

Benthic stack

The benthic d18O LR4 stack spans 5.3Ma. It is based on an average of 57 benthic d18O records (which measure global ice volume and deep ocean temperature), globally distributed sites along the Earth. Among with the ice cores (from Antartida and Greenland), the benthic stack provides an independent proxy, as a template for global clima and a stratrigraphic tool for paleoceanographers. For example, using the benthic stack we can identify the glacial cycles, climate cycles from Martrat et al. (2004) or marine isotopic stages (so called MIS).

References:
Lisiecki, L. E., and M. E. Raymo (2005), A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records,
Paleoceanography, 20, PA1003, doi:10.1029/2004PA001071

Here I show a simple code to download it directly from the internet (don't forget to cite), so we can use it in R. No extra packages needed. The plot shown is just a simple example for the first 500 ka (note also the way to write the legend, where d18O is on a correct syntax, I hope :) which I was unable to do in blogger).



# Benthic stack 
 # database created with 1000 corers around the world. 
 # http://lorraine-lisiecki.com/stack.html
 # Please cite: Lisiecki, L. E., and M. E. Raymo (2005), A Pliocene-
 # Pleistocene stack of 57 globally distributed benthic d18O records, 
 # Paleoceanography,20, PA1003, doi:10.1029/2004PA001071.  
 # readLines("http://lorraine-lisiecki.com/LR04stack.txt")
 benthic.stack=read.table("http://lorraine-lisiecki.com/LR04stack.txt", skip=5, col.names=c("age","d18O","se")); 
 benthic.stack$d18O=benthic.stack$d18O*-1
 plot(d18O~c(age), benthic.stack[benthic.stack$age<450,], typ="l", col=4, xlab="Age (ka)", ylab=expression(paste(delta^18~O," Benthic stack")))
 abline(v=seq(round(min(benthic.stack$age)),round(max(benthic.stack$age)), by=20), col=8, lty=3)
Created by Pretty R at inside-R.org



Thursday, 24 November 2011

Eccentricity, Obliquity and Precession

The Earth climate is primarily influenced by the Sun,
thus changes in the orbital variations are generally
considered in paleoceanographic studies. Thus,
my first step was to download this data.
In this code, you download the data directly from the net,
and plot it in R. No extra packages required.

#  Orbital Variations and Insolation Database: Insolation.readme file
 orbit=read.table("ftp://ftp.ncdc.noaa.gov/pub/data/paleo/insolation/orbit91", 
      skip=3, col.names=c("age", "ECC","OMEGA","OBL","PREC","I65NJul","I65SJan","I15NJul","I15SJan"))
 #summary(orbit)
 orbit$age=abs(orbit$age)
 par(mfrow=c(2,2), mar=c(4,4,1,1), cex=.7)
 plot(ECC~age, orbit[orbit$age<450,], typ="l", col=4, main="Eccentricity")
 plot(OMEGA~age, orbit[orbit$age<450,], typ="l", col=2)
 plot(OBL~age, orbit[orbit$age<450,], typ="l", col=3, main="Obliquity")
 plot(PREC~age, orbit[orbit$age<450,], typ="l", col=5, main="Precession")
Created by Pretty R at inside-R.org

 And than you get something like this: