TerraceM-2 – MATLAB software for the analysis und modelling of coastal and lake terraces

Posted on Updated on

Todays blog post is written by my colleage Julius Jara-Muñoz (University of Potsdam). Julius is the author of TerraceM, a software for the analysis and modelling of marine and lacustrine terraces.

TerraceM-2 is a Matlab® graphic-user interface to analyze and model marine and lacustrine terraces. TerraceM was originally designed to obtain accurate morphometric assessments of terrace morphology. The new TerraceM-2 version includes novel routines to map the elevation and spatial distribution of terraces, to model their formation and evolution, and to estimate fault-slip rates from terrace deformation patterns. TerraceM-2 uses nested Topotoolbox functions (Schwanghart and Scherler, 2014) to perform basic analyses, such as swath profile extraction, calculation between topographic parameters, and DEM visualization, significantly improving the processing speed and mapping capabilities of TerraceM-2.

TerraceM-2 comprises three modules, a) Mapping tools (MAP) b) Landscape evolution model (LEM), and c) Dislocation Model (DIM), designed to be used both as independent and complementary GUI’s. The MAP module is designed for mapping marine and lacustrine terraces using swath profiles and surface classification models. The LEM module is designed to model terraces formation simulating the effect of coastal erosion and vertical deformation, allowing comparing modelled and measured terraces in order to determine uplift and erosion rates. The DIM module uses the spatial distribution of terrace elevations or surface uplift, comparing them with fault dislocation models to determine fault slip rates.

1. Here is an example of marine terrace mapping using the TerraceM-2 MAP GUI

2. Here is an example of shoreline angle mapping using the TerraceM_staircase function.

%Add path to topotoolbox folder, TerraceM-2 Maptools, and dem and shapefile.
addpath /Users/TerraceM_2018/Stations/HUEVO
addpath /Users/TerraceM_2018/topotoolbox-master
addpath /Users/TerraceM_2018/TerraceM_Maptools

Box=shaperead('Huevo_clip.shp'); %Shapefile containing the profiles boxes 
DEM=GRIDobj('DEM_huevo.tif'); %DEM generated using Topotoolbox
profnum=2; %Profile number to process

%Create swath function

%Display swath profile 
xlabel('Distance along profile (m)'); ylabel('Elevation (m)');
legend('Min elevation','Mean elevation','Max elevation')

%Mapping a single shoreline angle
Swath profile generated using the TerraceM_makeswath function and mapped using the TerraceM_staircase function.

3. Here is another example of how to create a LEM movie using the TerraceM_anderson function.

%add path to TerraceM-2 LEM and the sealevel folders
addpath /Users/TerraceM_2018/TerraceM_LEM
addpath /Users//TerraceM_2018/TerraceM_LEM/sealevel
sl=load('Bintanja_total.mat');%load a sea-level curve
%parameters in structure format
param.tmax=500;%max time (ka)
param.tmin=0;%min time (ka)
param.dt=100;%bin size of time intervals (yrs)
param.dx=5;%bin size of x axis (m)
param.sill=150;%m, avoid sill by adding an extremely deep value
param.slope_shelf=10;%initial slope model (deg)
param.uplift_rate=1.8; %m/ka
param.initial_erosion=0.3;% initial erosion rate or E0 (m/yr)
param.cliff_diffusion=0.02;% cliff slope difussion (m^2/yr)
param.wave_height=5;%wave height (m)
param.vxx=0;%manually extend x axis if model run out of bounds
param.smth=40;%smooth factor of sea-level curve (yr)
param.video=1;%video mode activated
%optional input parameters, ignore for quick video display
param.cumula=0;%cummulative plot deactivated
param.snap=90;%number of video frames
%create axes for video
%call LEM function

The integration of TerraceM-2 modules and routines, in addition to user-friendly GUIs, optimizes the time required for mapping and modelling marine and lacustrine terraces for inexperienced users in Matlab® scripting. Furthermore, the TerraceM-2 architecture comprises independent functions allowing experienced users creating custom scripts and workflows beyond the GUI environment.

TerraceM-2 is available at www.terracem.com, there you can also find tutorials and examples, updates, and a list of recent publications using TerraceM.


Schwanghart, W., and Scherler, D. (2014). Topotoolbox 2–MATLAB-based software for topographic analysis and modeling in earth surface sciences. Earth Surf. Dyn. 2, 1–7. doi: 10.5194/esurf-2-1-2014

Jara-Muñoz, J., Melnick, D., Pedoja, K., and Strecker, M.R. (2019). TerraceM-2: A Matlab® Interface for Mapping and Modeling Marine and Lacustrine Terraces. Frontiers in earth Science 7(255). doi: 10.3389/feart.2019.00255.


Open PhD position in landscape evolution modelling

Posted on Updated on

Numerical models of landscape evolution (LEMs) are now increasingly sophisticated and able to model the interplay of various geomorphic processes and their controls. Parametrizing and testing these models using empirical data, however, remains challenging.

Markus Fuchs, professor for physical geography at the University of Giessen, Germany, now offers a PhD position. Equipped with a vast array of empirical data, the prospective candidate will aim at tackling some of the pressing methodological issues on model-data comparison that pertain validation and parametrization of LEMs. I have agreed to co-supervise the PhD student. If you are about to finish your Master degree and have skills in numerical computing and statistical analysis, we cordially invite you to apply for this position. If you are teacher or supervise students, then it’ll be great if you forward this open position.

The official job announcement in german and english language can be found here.


Visualizing trends in sea ice extent

Posted on

Your ability to speak a foreign language will improve with practice. Your ability to code MATLAB, too. Thus, I sometimes write some code for fun and do some visualizations. This has the advantage that I stay up to date with new developments in the graphics engine of MATLAB.

My latest code is on visualization of ice extents and their temporal dynamics and I decided to post this code here. What you need is two more function that you can download from the file exchange: gif by Chad Greene and lbmap by Robert Bemis. lbmap provides a bipolor colormap and gif allows us to write an animated gif image. The data that we visualize is from the National Snow and Ice Data Center and can be downloaded here.

data = readtable('Sea_Ice_Index_Monthly_Data_by_Year_G02135_v3.0.xlsx');
data = removevars(data,{'Var14' 'Annual'});
year = data.Var1;
month = 1:12;
monthstr = data.Properties.VariableNames(2:end);

[year,month] = ndgrid(year,month);
A    = table2array(data(:,2:end));

year = year';
month = month';
A    = A';

monthmean = mean(A,2,'omitnan');
monthdev  = (A - monthmean)./monthmean * 100;

year = year(:);
month = month(:);
A    = A(:);
t    = datetime(year,month,repmat(15,size(year)));

% Remove leading and trailing nans
I    = ~isnan(A);
I(find(I,1,'first'):find(I,1,'last')) = true;

year = year(I);
month = month(I);
A    = A(I);
t    = t(I);
monthdev = monthdev(I);

theta = (month)/12*2*pi;

% Animated figure

ncol     = 500;
cmap     = lbmap(ncol,'BrownBlue');

monthdev = monthdev(:);
monthdev = min(monthdev,20);
monthdev = max(monthdev,-20);
minmonthdev = min(monthdev(:));
maxmonthdev = max(monthdev(:));
colix    = interp1(linspace(minmonthdev,maxmonthdev,ncol),1:ncol,monthdev,'nearest');

paxanim = polaraxes('ThetaDir','clockwise',...
    'Rlim',[floor(min(A)) ceil(max(A))],...
    'Clim',[minmonthdev maxmonthdev]);

h = colorbar;
h.Label.String = 'Rel. deviation from average [%]';
h.Location = 'southoutside';

% Some figure properties
paxanim.Color = [.5 .5 .5];
fax = gcf;
fax.Color = [1 1 1];

hold on
makegif   = true;
inigif    = true;

for r = 2:numel(A)
    if ~isnan(A(r))
        line(paxanim,theta([r-1 r]),A([r-1 r]),'color',cmap(colix(r),:));
        poi = line(paxanim,theta(r),A(r),'Marker','o',...
        if inigif && makegif
            inigif = false;
        elseif ~inigif && makegif && mod(r,4)==0
        if r == numel(A)
            poi = line(paxanim,theta(r),A(r),'Marker','o',...
            for r2 = 1:50

And here are the results. Have fun doing your own visualizations and animations.

Sea ice extent on the Northern hemisphere between 1978 and today. Values along the radius have units of 1 Mio. sqkm. Line colors indicate deviations [%] from the long-term (1978-2019) monthly averages.

CDT – The Climate Data Toolbox for MATLAB

Posted on

CDT logo illustration by Adam S. Nelsen (http://www.adamnelsen.com)

Today’s post is about CDT, the Climate Data Toolbox for MATLAB. The development of this new toolbox has been speerheaded by Chad A. Greene, an active long-time contributor to the MATLAB community and scientist at the Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. CDT has been released and an accompanying paper has just been published (Greene et al., 2019).

CDT is special in many respects. First of all, it contains a versatile set of MATLAB functions for analyzing and displaying climate data, be it gridded, spatiotemporal or time series data. The functions are computationally efficient and quite easy to use. Second, CDT comes with many tutorials that describe not only how to use CDT functions, but also detail mathematical backgrounds, and offer guidance on how to interpret the results in the context of Earth science processes.

As geomorphologist, I am often working with climate data. In fact, my diploma thesis was on dust mobilization in the Sahara and its large-scale meteorological controls (Schwanghart and Schütt, 2008). In retrospective, CDT would have helped me a lot in learning the computational tools required for the analysis. Too bad, CDT was not available at this time.

Give CDT a try. Chad Greene’s website has the documentation online. If you install CDT (download from github or the MATLAB file exchange), this documentation will be available in the MATLAB help browser, too.


Greene, C.A., Thirumalai, K., Kearney, K.A., Delgado, J.M., Schwanghart, W., Wolfenbarger, N.S., Thyng, K.M., Gwyther, D.E., Gardner, A.S., Blankenship, D.D., 2019. The Climate Data Toolbox for MATLAB. Geochemistry, Geophysics, Geosystems, accepted [DOI: 10.1029/2019GC008392]

Schwanghart, W., Schütt, B. (2008): Meteorological causes of Harmattan dust in West Africa. Geomorphology, 95, 3-4, 412-428. [DOI: 10.1016/j.geomorph.2007.07.002]

TopoToolbox dependencies

Posted on

TopoToolbox is a MATLAB software. MATLAB consists of the core product MATLAB and can be extended by numerous toolboxes that have specific purposes, and that cost extra. TopoToolbox also comes as a toolbox, and is free. When typing ver in the command line, TopoToolbox should be listed together with the other available toolboxes.

TopoToolbox relies mainly on core MATLAB functions, but unfortunately it also requires a few toolboxes. Many universities and academic institutions may have licences for these toolboxes, but if you are working for a private company or some governmental administration, access to these toolboxes may be limited, often due to monetary reasons. Accordingly, usability of TopoToolbox may be limited for many.

Dirk and I are well aware of this problem. Thus, we have designed TopoToolbox in a way that keeps toolbox dependencies at a minimum. However, this is not completely feasible. This is because, first, we only have limited time available for our TopoToolbox work. Writing some functions of other Toolboxes ourself would certainly be possible, but it is likely that they would not run with the same efficiency and speed. It would also mean reinventing the wheel. Second, we could make use of external and free third-party toolboxes and libraries that are available on the MATLAB file exchange and github. However, that would also mean maintaining these dependencies which can be difficult if there are major changes. Of course, MATLAB toolboxes are also changed frequently, but these changes are usually well documented in the release notes. Thus, to secure a positive user-experience we largely refrain from using third-party toolboxes.

What are toolbox dependencies of TopoToolbox?

DEM analysis heavily relies on image processing algorithms. In fact, gridded DEMs are just georeferenced gray-scale images. Thus, TopoToolbox requires the Image Processing (IP) Toolbox which provides access to numerous filtering and morphological algorithms (e.g. image erosion and dilation, gray-weighted distance transforms). If you don’t have the IP toolbox, you will quickly encounter errors that state that a particular function is missing. If you don’t have the IP Toolbox but you want to still use TopoToolbox, I’d recommend using dipimage, an excellent image processing toolbox. But this would mean making numerous changes to the code. I won’t say it is impossible, but it’ll be quite some work.

Then, you should have the Mapping Toolbox. If you don’t have it, reading and writing geocoded data (geotiffs, shapefiles) is a bit tedious. Moreover, some of the functions (e.g. STREAMobj2latlon) would not work because they rely on the Mapping Toolbox functions that transform coordinates from projected to geographic systems and vice versa. TopoToolbox works well without georeferencing, but having the Mapping Toolbox is much more convenient.

Some new functions (crs, smooth (only if output is positive only), quantcarve) need the Optimization Toolbox because they rely on large-scale sparse linear and quadratic programming algorithms. If you don’t have the Optimization Toolbox, you won’t be able to run these functions. There are some third-party optimization products available, but I haven’t looked at alternatives so far. If you know good alternative toolboxes that are free and implement large-scale sparse optimization techniques, please let me know. mnoptim uses the bayesopt function that come with the Statistics and Machine Learning Toolbox.

Finally, some of the functions are coded so that they make use of the Parallel Computing Toolbox. However, they run without this toolbox, too.

All in all, to run TopoToolbox you should definitely have the Image Processing Toolbox and the Mapping Toolbox. Some functions require other toolboxes, but perhaps you don’t really need these functions. Hope I didn’t forget any dependencies. If I did, please let me know.

Meet us at the EGU 2019

Posted on Updated on

Point pattern and kernel density estimate on the stream network of Big Tujunga.

Dirk and I will attend the EGU this year. And as luck would have it, our presentations on point processes on river networks (PICO4.7 | EGU2019-18489) and on identification and ordering of divides (PICO4.6 | EGU2019-14979) are in the same PICO session (Tue, 09 Apr, 10:55–10:59, PICO spot 4). Interactive presentation time is afterwards from 11:15-12:30.

Stream (a) and divide (b) network of the Big Tujunga catchment.

If you want to chat and discuss with us, please come to our interactive presentations. We will highlight some of the upcoming additions to TopoToolbox, and we are curious what you think of them.

iceTEA – Tools for Exposure Ages from ice margins

Posted on


My guest blogger today is Richard Selwyn Jones from Durham University. He developed iceTEA, a set of MATLAB tools to calculate exposure ages. He kindly accepted my invitation to write here about his software. Thanks!

iceTEA (Tools for Exposure Ages) is a new toolbox for plotting and analysing cosmogenic-nuclide surface-exposure data. The tools were originally intended for data that are used to reconstruct past glacier and ice sheet geometries, but they can potentially be used for many other geochronology and geomorphology applications of cosmogenic nuclides.

The iceTEA tools are available via an online interface (http://ice-tea.org) and as MATLAB code with a ready-to-use front-end script for each tool (https://github.com/iceTEA-code). While the online version performs all primary analysis and plotting functionality, the code provides the user with greater flexibility to apply the tools for specific needs and also includes some additional options.

There are eight tools, which provide the following functionality: 1) calculate exposure ages from 10Be and 26Al data, 2) plot exposure ages as kernel density estimates and as a horizontal or vertical transect, 3) identify and remove outliers within a dataset, 4) plot nuclide concentrations on a two-isotope diagram and as a function of depth, 5) correct exposure ages for cover of the rock surface, 6) correct ages for changes in elevation through time, and estimate 7) average and 8) continuous rates of change (e.g. ice margin retreat or thinning).

Here is an example of how you can use the code to correct data for long-term uplift:

data_name = 'input_data.xlsx';  % File name used for sample data

scaling_model = 'LSD';     % 'DE','DU','LI','ST','LM','LSD','LSDn'
% Set elevation correction
correction_type = 'rate';  % Select 'model' or 'rate'
GIA_model = [];            % If 'model', set GIA model to use - 'I5G' or 'I6G'
elev_rate = 2;             % If 'rate', set rate of elevation change (m/ka)

sample_data = get_data(data_name);  % Load sample data
% Calculate Elevation-corrected Exposure Ages
if strcmp(correction_type,'model')
    elev_input = GIA_model;
elseif strcmp(correction_type,'rate')
    elev_input = elev_rate;
corrected = elev_correct(sample_data,scaling_model,correction_type,elev_input);

The corrected exposure ages can then be plotted as kernel density estimates:

% Plot settings
feature = 1;    % Data from single feature?  yes = 1, no = 0
save_plot = 0;  % Save plot?  1 to save as .png and .eps, otherwise 0
mask = [];      % Select samples to plot (default is all)
time_lim = [];  % Optionally set x-axis limits (in ka)
weighted = [];  % Optionally select weighted (1) or unweighted (0) mean and standard deviation (default is weighted)
% Plot figure

A two-isotope diagram. Samples that plot inside the “simple exposure region” (marked by black lines) were continuously exposed, whereas samples that plot below this region have been buried in the past with a complex exposure history.
A kernel density estimate plot for exposure ages from a moraine. Individual ages are shown in light red, with the summed probability as a dark red line. One of the tools allows the user to evaluate the effects of surface cover on their dataset using different types and depths of surface cover (shown here as green summed probability lines for snow and till).
Rates of change can be estimated using linear regression in a Monte Carlo framework, or continuously using Fourier Series analysis or Bayesian penalized spline regression. The latter is shown here. The exposure age constraints (with age and sample position uncertainties) are modelled (upper panel) in order to generate rates of change (lower panel). In this example, ice sheet thinning was most rapid at approx. 8 ka.

The purpose of these tools is to allow users to explore the spatial and temporal patterns in their data in a consistent and inter-comparable way, and also to initiate discussion of further improvements in the application and analysis of surface-exposure data. We welcome suggestions for additional plotting or analysis tools.


Jones, R.S., Small, D., Cahill, N., Bentley, M.J. and Whitehouse, P.L., 2019. iceTEA: Tools for plotting and analysing cosmogenic-nuclide surface-exposure data from former ice margins. Quaternary Geochronology, 51, 72-86. [DOI: 10.1016/j.quageo.2019.01.001]