MATLAB Code for Earth and Environmental Science | Last updated 18 November 2009
These are functions or blocks of scripts that I've written or hacked over the course of my research. I'm making these available on the off-chance others might find them useful. If you do, please drop me a line and let me know.
PIN - Preindustrial correction for carbon isotope tree-ring series (d13C) for pC02
PIN - PreINdustrial correction for carbon isotope tree-ring series
A user-written Matlab function that performs a correction on carbon
isotope series from tree-rings to account for changes in atmospheric
pCO2 during the last ~150 years.
The function is based on the physiological mechanisms and logical
constraints discussed in detail in McCarroll et al. [2009].
RWL2MAT - a function to read Tucson decadal format files into Matlab
RWL2MAT a function to read Tucson decadal format files into Matlab
This function reads in the selected file from Tucson decadal format
and makes it available in the Matlab workspace. It is simpler, but
less powerful, than David Meko's RWL2TSM.
MAT2RWL - a function to write Tucson decadal format files from Matlab
MAT2RWL a function to write Tucson decadal format files from Matlab
mat2rwl(data,years,seriesnames,filename)
This function writes a matrix into a Tucson decadal format ASCII file.
Missing values should be 'NaN'.
EBISUZAKI - a method for estimating significant correlation values for autocorrelated time series
EBISUZAKI a method for estimating significant correlation values for
autocorrelated time series
[F] = ebisuzaki(x,y,sig,nsim)
This function creates 'nsim' random time series that have the same power
spectrum as the original time series but with random phases.
Input
x : vector of (real) numbers with the first time series
y : vector of (real) numbers with the second time series
sig : [optional] significance level for critical value estimation [0.05]
nsim : [optional] the number of simulations [10000]
Output
F : Fraction of artificial time series with higher correlation coefficents
References
----------
Ebisuzaki, W, 1997: A method to estimate the statistical
significance of a correlation when the data are serially correlated.
J. of Climate, 10, 2147-2153.
modified by Kevin Anchukaitis from a Matlab function for creating
random phase time series by Vincent Moron and from the original C code by W. Ebisuzaki
Ebisuzaki's C code: ftp://ftp.cpc.ncep.noaa.gov/wd51we/random_phase
GSBTEST - a method for estimating the significance of moving correlation analyses
GSBTEST - Monte Carlo Test of Moving Correlation Significance Levels following Gershunov et al. 2001
Implementation of a Monte Carlo type simulation of moving correlation
functions between two weakly correlated white noise series for the purposes of constructing
significance test for moving correlation functions based on their standard deviation (Gershunov et al. 2001).
Inputs:
sL = length of the time series used in the running correlation analysis
sN = number of Monte Carlo iterations to perform (should be >= 1,000)
window = width (in units of the time step for the data) of the correlation window
dR = desired overall correlation of the simulated time series (simulations will be with 1 of this value)
method [optional] = Which method to use to create correlated random series [default = 1], see below
Correlation Methods
1 = Cholsky decomposition of the covariance matrix [default]
2 = Matrix square root
NOTE: Although the solution should converge for sN >= 1000, small differences (due to the correlation between
the two artificial time series generated in the script), may still result, especially for a small number (<<10000)
of simulations; so it is recommended that you verify convergence with several runs and be cautious of values too
close to your chosen test statistic (particularly at 90 or 95 CI). You can compare the above result for 'sd95'
to Table 1 from Gershunov et al. 2001 (interseries r = 0.30, window length = 11, 95 confidence level = 0.36).
References
----------
Alexander Gershunov, Niklas Schneider, and Tim Barnett, 2001. Low-frequency modulation of the ENSO-Indian
monsoon rainfall relationship: Signal or noise?, Journal of Climate, 14: 2486-2492.
Updated 18 November 2009 to use pairwise.m
Tips and Tricks for Scientific Computing with OSX
This is mostly a collection of LaTeX, Matlab, FORTRAN, or NCL code, assorted computing tricks, and other random bits of information that I'm trying to remember and keep track of, but hopefully you'll find them useful too. If you do, drop me a line and let me know.
Description of the Problem: Despite having the correct libraries installed, running NCL version 5.1, OPENDaP-enabled for Intel Mac OSX, results in the following error:
dyld: Library not loaded: /Users/haley/tmp/5.1.0/lib/libnc-dap.3.dylib
The solution is eventually described by piecing together clues from this thread at the NCAR NCL discussion forum. If you have installed the correct libraries (libdap, libnc-dap, and ssl) via Macports, for instance, you should include the following line in one of your .bashrc or .profile files:
export DYLD_LIBRARY_PATH=/opt/local/lib
You'll also need to install gfortran, but if you get it from HPC, NCL will find it.
UPDATE:However, setting DYLD_LIBRARY_PATH seems to mess up Apple's native X11 (which Matlab uses. Doh.)
UPDATE 2: The trick is to set DYLD_FALLBACK_LIBRARY_PATH, which will allow X11 to see the default Apple path, but also let NCL find the libraries it wants in /op/local/lib. So:
[5] Installing mexnc and OPeNDAP for Matlab on Mac OSX
Description of the Problem: This is mostly just for me, since I've now done this half a dozen times, each time just long enough after the previous time to have forgotten the little details. But maybe you'll also find this helpful.
[1] Download the precompiled binaries for mexnc here. Make sure the correct files can be found on your matlabpath.
[2] I still use Chuck Dedham's object-oriented netcdf toolbox, which can be found here.
[3] You can setup OPeNDAP for Matlab using the new Matlab OPeNDAP Ocean Toolbox. Previously, you had to install libdap and libxml2, in addition to the Matlab OPeNDAP structure tools. Now it appears the Ocean Toolbox takes care of everything (although you might need to install libxml2 or libdap still?). The Ocean Toolbox appears to also put the proper files in the proper directory.
[4] There are some bugs in the loaddap and whodap functions used in Matlab to interface with OPeNDAP. First, you'll have to edit the function whodap.m, and remove the './' from in front of the system call to writedap. Second, the URL for the testing files is wrong in the documentation. The proper testing URL command is:
[5] Rob Hetland at Texas A&M University has notes on building your own netCDF for OSX here.
Tags: netcdf, Matlab, mexnc, OPeNDAP
[4] Using netCDF operators (NCO) to concatenate files without a record dimension
Description of the Problem: Imagine you have a set of netcdf files, each with a single month of a multiyear dataset. In some cases, the time dimension will not be the record dimension. This means you can't use the netcdf (NCO) operator ncrcat to concatenate the individual months along the time dimension. But, there is a work around:
Solution: run the following commands in a bash shell script [you can get mine here: makeSingleNC.sh]
list=`ls *.nc`
for i in $list;
do
ncecat -O $i $i
ncpdq -O -a time,record $i $i
ncwa -O -a record $i $i
done
ncrcat -h *.nc out.nc
Tags: netcdf, nco, .nc, no record dimension, ncrcat, ncecat
[3] Matlab m_map contours not filled
Description of the Problem: Filled contour maps using m_contourf from the m_map mapping package are not filled
Solution: Insert the following command after the call to m_grid
[2] OSX (Tiger) Updates and Security Updates Break Kerberos, SSH, and (sometimes) Keychain
Description of the Problem: Following Software Update installation of OSX Security and Software Updates, Kerberos and SSH are broken. One symptom is that ssh attempts result in a 'segmentation fault'.
Solution: The current (best? only?) description of the solution for this problem is from Eric Salathe, in this thread from the Apple Discussion Forums.
There is a discussion and work-around of the issue of Matlab mishandling permil notation here.
However, I found that this doesn't specifically address the problem as it is manifesting itself on my Mac Intel, running R2007b. Here is a description of my solution:
Here is what I entered:
=======================
ylabel(['\delta^{18}O (',char(8240),' VSMOW)'])
Here is how matlab interprets the output figure with the permil:
================================================================
%%IncludeResource: font Helvetica
/Helvetica /ISOLatin1Encoding 120 FMSR
516 1779 mt -90 rotate
(O \(\32 VSMOW\)) s
90 rotate
Here is how matlab should have interpreted it in the output:
============================================================
%%IncludeResource: font Helvetica
/Helvetica /StandardEncoding 120 FMSR
516 1779 mt -90 rotate
(O \(\275 VSMOW\)) s
90 rotate
The key here is to change the encoding to '/StandardEncoding' and the character code itself from '\32' to '\275'.
