This is a simple script that demonstrates how to open netcdf files (a typical format of file used for storing large amounts of data, and often used to display output from 3D Earth system models). This example uses a version of the marine reservoir ages output from Butzin et al. 2017.
In [13]:
# load required packages
library(ncdf4)
library(maps)
In [14]:
#Open netcdf file
nc <- nc_open( "mra14_intcal13_pd_21kcalBP.nc")
In [15]:
#list names of variables in netcdf file
names(nc$var)
In [16]:
#list names of dimensions in netcdf file
names(nc$dim)
Lets say that you want to extract the marine reservoir ages at a specific location for which you have the longitude, latitude and depth (e.g. -55oN, -70oE, 3000m). The following is the code for how to extract such data.
In [17]:
#Lets have a look at the matrix containing the variable of interest. In this case marine reservoir ages ("MRA)
nc_var <- ncvar_get( nc, varid="MRA" )
#list how many dimensions the matrix has
dim(nc_var)
In [19]:
#Note how this compares to the number of entries in each dimension
length(nc$dim$lon$vals)
length(nc$dim$lat$vals)
length(nc$dim$depth$vals)
In [41]:
#The matrix containing our variable is made up of the following [lon,lat,depth]
#Lets try to extract the data for our core site
#First define the variables
input_lat = -55.1
input_lon = -70
input_depth = 1250
#Now find the location in the matrix that corresponds to our data.
#It is likely that your core is not at the exact location as each data point in the cdf file,
#so you will have to find the nearest grid point
#Find correct colours for Interpolated_masterfile
nc_lat<-nc$dim$lat$vals
nc_lon<-nc$dim$lon$vals
nc_depth<-nc$dim$depth$vals
#This is written in a loop to make it easier when you have more than one site
index_vals=NULL
for(i in 1:length(input_lat)){
lat_index<-which(abs(nc_lat-input_lat[i])==min(abs(nc_lat-input_lat[i])))
#longitudes may need correcting from -180 to 180 ---> 0 to 360
input_lon2<-ifelse(input_lon[i]< min(nc$dim$lon$vals), input_lon[i]+360, input_lon[i])
lon_index<-which(abs(nc_lon-input_lon2)==min(abs(nc_lon-input_lon2)))
depth_index<-which(abs(nc_depth-input_depth[i])==min(abs(nc_depth-input_depth[i])))
a<-data.frame(lat_index=lat_index,lon_index=lon_index,depth_index=depth_index)
index_vals<-rbind(index_vals,a)
}
index_vals
#Now use these index values to find the RMA at the core site by using the variable matrix
MRA_output=NULL
for(i in 1:nrow(index_vals)){
MRA<-nc_var[index_vals$lon_index[i],index_vals$lat_index[i],index_vals$depth_index[i]]
MRA_output<-rbind(MRA_output,MRA)
}
MRA_output
The same principle can be used to extract marine reservoir ages from multiple sites. You can import in an excel sheet containing all of your longitudes, latitudes and depths, and output the data as a csv file.
In [44]:
#This package lets you read in excel docs
require(gdata)
input_cores<- read.xls("All_chilean_margin_cores.xlsx", sheet=1, header=TRUE)
#See the top few lines of your excel file
head(input_cores)
In [96]:
#Define the input variables
input_lat = input_cores$Lat..oN.
input_lon = input_cores$Long..oE.
input_depth = input_cores$WD..m.
#Again apply the function to find the index locations of each site within the variable matrix
index_vals=data.frame(lat_index=numeric(0),lon_index=numeric(0),depth_index=numeric(0))
for(i in 1:length(input_lat)){
lat_index<-which(abs(nc_lat-input_lat[i])==min(abs(nc_lat-input_lat[i])))
#longitudes may need correcting from -180 to 180 ---> 0 to 360
input_lon2<-ifelse(input_lon[i]< min(nc$dim$lon$vals), input_lon[i]+360, input_lon[i])
lon_index<-which(abs(nc_lon-input_lon2)==min(abs(nc_lon-input_lon2)))
depth_index<-which(abs(nc_depth-input_depth[i])==min(abs(nc_depth-input_depth[i])))
a <- data.frame(lat_index=lat_index,lon_index=lon_index,depth_index=depth_index)
index_vals<-rbind(index_vals,a)
}
#Now use these index values to find the RMA at the core site by using the variable matrix
MRA_output=NULL
for(i in 1:nrow(index_vals)){
MRA<-nc_var[index_vals$lon_ind[i],index_vals$lat_ind[i],index_vals$depth_ind[i]]
MRA_output<-rbind(MRA_output,MRA)
}
#Add the output to the original datafile
input_cores["MRA"]<- MRA_output
#Display datafile
input_cores
#Output as csv file
write.csv(input_cores,"input_cores_with MRA.csv")