Search
Feeds
Archives
About
Links
BLOG
Sunday 28 July 2019 at 09:06 am.This section presents some examples of the command line capabilities of TerraceM outside the Graphic user interface.
Loading TerraceM-2 Maptools outputs
TerraceM-2 Maptools creates several Matlab® files used for data exchange between the other TerraceM-2 modules. These files can be called independently to create customized plots. In this example, we load the main output files containing the shoreline angles and swath profiles.
%Load swath profiles
load('HUEVE_swath.mat');
%Load shoreline angles
load('HUEVE_SHORELINES.mat');
The shoreline-angle file (.mat) is a structure containing several fields that define: a) method (met) used to map the shoreline angles, including the linear regression coefficients and points used for mapping, and b) the shoreline angle level (level), which includes the profile, elevation, error, map coordinates, location along profile and the station or mapping site name.
Displaying swath profiles and shorelines angles mapped in TerraceM-2
This excercise follows the previous example.
Here we present an example plotting swath profiles and shoreline angles using the output files of TerraceM-2 Maptools.
i=2; %profile number
j=1; %marine terrace level
%these values can be eventually looped to plot all mapped shoreline angles
%extract the linear regressions for platform and cliff
cliff=SHORELINE(i).met(j).REGRESSIONS.cliff;
plat=SHORELINE(i).met(j).REGRESSIONS.plat;
points=SHORELINE(i).met(j).POINTS;
%extract the shoreline angle
SH(1)=SHORELINE(i).level(j).DISTANCE_PROFILE;
SH(2)=SHORELINE(i).level(j).ELEVATION;
SH(3)=SHORELINE(i).level(j).ERROR;
%plot the result
figure
hold on; plot(swath(i).x,swath(i).z_max,'-k',swath(i).x,swath(i).z_min,...
'-k',swath(i).x,swath(i).z_mean,'-k');%plot the swath profile
plot(cliff(:,1),cliff(:,2),'--r'); plot(cliff(:,1),cliff(:,4),'--r');%cliff points
plot(cliff(:,1),cliff(:,3),'-r');plot(plat(:,1),plat(:,2),'--r');%platform points
plot(plat(:,1),plat(:,4),'--r');plot(plat(:,1),plat(:,3),'-r');%regressions
errorbar(SH(1),SH(2),SH(3),'ok'); plot(SH(1),SH(2),'ok',... %shoreline angle
'MarkerFaceColor','r','Markersize',6);
xlabel('Distance along profile (m)');ylabel('Elevation (m)');%axis labels
Creating swath profiles from shapefile polygons outside the GUI environment
Swath profiles can be generated using the TerraceM_makeswath.m function and Topotoolbox. The inputs are a DEM in geotiff (.tif) or Esri grid format (.asc) and a shapefile containing the swath profiles boxes, which may be created using the rectangle tool. The output is a structure containing the distance, minimum, mean, maximum and slope of a swath profile.
%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
swat=TerraceM_makeswath(DEM,Box,profnum);
%Display swath profile
figure
plot(swat.x,swat.z_min,'-c',swat.x,swat.z_mean,'-b',swat.x,swat.z_max,'-k');%plot the swath lines
xlabel('Distance along profile (m)'); ylabel('Elevation (m)');
legend('Min elevation','Mean elevation','Max elevation')
Mapping shoreline angles outside the GUI environment
This excercise follows the swath profile plot excercise presented before.
We can directly map shoreline angles calling the function TerraceM_staircase.m, which displays an interactive window to map the shoreline angle by selecting the paleo-platform and paleo-cliff. Outputs include the shoreline angle elevation and its associated error (shoreline), the points used to fit the linear regressions (points), and the linear regression coefficients and 95% confidence bounds (regr).
%Mapping a single shoreline angle
[shoreline,points,regr]=TerraceM_staircase(swat);
TerraceM-2 LEM scripts
The TerraceM-2 LEM module includes two functions that can be used to create custom models, export results in video formats, and compare the results to swath profile topography.
Creating a LEM movie
Two options exist to create a LEM movie: Simple or customized display. The optional input parameters allow for custom movie size, style, and speed.
clear
%add path to TerraceM-2 LEM and the sealevel folders
addpath /Users/TerraceM_2018/TerraceM_LEM
addpath /Users//TerraceM_2018/TerraceM_LEM/sealevel
%inputs
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
param.h1=subplot(3,1,1:2); param.h2=subplot(3,1,3);
%call LEM function
[LEM]=TerraceM_anderson(sl,param);
Comparing measured and synthetic LEM terraces
The topographic expression of marine terraces can be directly compared to LEM results by tying them using a single shoreline angle. This requires assigning an estimated or measured age to the shoreline angle and adjusting the input parameters to minimize the weighted root mean squared error WRMSE or the difference between model and observations. To compare modeled and measured topography, we use the TerraceM_rms_fit.m function and a swath profile .mat input file, and select a piercing point shoreline angle. The function uses two LEM’s, one including all the user-defined parameters and the second with slope diffusion equal to 0, in order to find the position of shoreline angles in the LEM. Here we present an example fitting a profile with a constant uplift rate (Fig 5SP); note that the quality of the fit decreases upwards suggesting that this area was probably affected by variable uplift rate along time.
clear
%addpath to TerraceM-2 LEM folders and input files
addpath /Users/TerraceM_2018/TerraceM_LEM
addpath /Users/TerraceM_2018/TerraceM_LEM/sealevel addpath /Users/TerraceM_2018/Stations/HUEVO/HUEVE
%First load a profile and the shoreline angle
load('HUEVE_swath.mat') %loadswath profiles
load('HUEVE_SHORELINES.mat') %load shoreline angles
i=9; %profile number
j=1; %marine terrace level
%extract the shoreline angle
SH(1)=SHORELINE(i).level(j).DISTANCE_PROFILE; SH(2)=SHORELINE(i).level(j).ELEVATION; SH(3)=SHORELINE(i).level(j).ERROR;
%extract the swath profile
SW(:,1)=swath(i).x;
SW(:,2)=swath(i).z_min;
SW(:,3)=swath(i).z_mean;
SW(:,4)=swath(i).z_max;
%create LEM’s
%LEM input parameters in structure format
sl=load('Bintanja_total.mat');%load a sea-level curve
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;%avoid sill effect of sl by adding an extremely deep value (m)param.slope_shelf=4;%initial slope model (deg)
param.uplift_rate=0.8; %uplift rate m/ka
param.initial_erosion=0.05;% initial erosion rate or E0 (m/yr)
param.cliff_diffusion=0.005;% cliff slope diffusion (m^2/yr)
param.wave_height=5;%wave height (m)
param.vxx=0;%manually extend x axis if model run out of bounds (m)
param.smth=40;%smooth factor of sea-level curve (yr)
param.video=0;%video mode deactivated
param2=param; %modify diffusion to cero for the default LEM
param2.cliff_diffusion=0;% cliff slope diffusion for the default model (m^2/yr)
%Create a LEM
[LEM]=TerraceM_anderson(sl,param);
%Create default LEM without slope diffusion to locate the shoreline angles
[LEM2]=TerraceM_anderson(sl,param2);
%LEM fitting parameters
age=125;%estimated shoreline angle age (ka)
ratio=100;%search ratio to avoid overlapping shoreline angles in LEM (yr)
cycle=15000;%mean span between glacial cycles (yr), smaller value increase detail in shoreline angle detection
%profile subset for fitting
up=9;%distance inlands from shoreline angle (km)
down=3;%distance seawards from shoreline angle (km)
mode=1;%plot fit activated
savepdf=0;%save figure deactivated
%%%% run the fit function
[Result,MOD]=TerraceM_rms_fit(LEM,LEM2,i,SW,SH,age,ratio,cycle,up,down,mode,savepdf);
TerraceM-2 DIM scripts
Creating a custom elastic dislocation model using mapped faults in a polyline shapefile as input
Elastic dislocation models can be developed using inputs as Shapefiles and the function TerraceM_okada.m; here we present an example using two faults.
clear
%set directories to TerraceM-2 okada functions and input files
addpath /Users/TerraceM_2018/stations/San_Andreas/SANDREAS/SANDREAS_EXPORT
addpath /Users/TerraceM_2018/TerraceM_DIM
%Input shapefiles created in arcgis
%faults, a polyline shapefile containing the faults
Sf=shaperead('SANDREAS_fault.shp');% faults polyline file, contains 2 faults
%box is a polygon shapefile of the area for modeling, usually is a square
%and should be bigger than the fault extent
Bx=shaperead('Sandreas_box.shp');
%Input parameters in structure format
param.d=100;%cellsize (m)
param.Dip(1)=80;%dip of fault 1 (deg)
param.Dip(2)=80;%dip of fault 2 (deg)
param.Slip(1)=5;%Slip fo fault 1 (m or m/ka)
param.Slip(2)=25;%Slip fo fault 2 (m or m/ka)
param.Updip(1)=0;%Updip of fault 1 (m)
param.Updip(2)=0;%Updip of fault 2 (m)
param.Downdip(1)=10000;%Downdip of fault 1 (m)
param.Downdip(2)=10000;%Downdip of fault 2 (m)
param.Rake(1)=90;%Rake of fault 1 (deg)
param.Rake(2)=90;%Rake of fault 2 (deg)
param.pplot=1;%plot result activated
%call the function to create dislocation models
[E3,N3,El,Nl,U]=TerraceM_okada(Sf,Bx,param);
