Visualization

Multiple colormaps in one axis

Posted on Updated on

I recently mentioned that MATLAB now lets you easily use different colormaps in one figure. The trick is to provide the axis handle as first input argument to the colormap function call.

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
ax(1) = subplot(2,1,1);
imagesc(DEM);
colorbar;

ax(2) = subplot(2,1,2);
imagesc(DEM);
% Use the axis handle as first input argument
% to colormap
colormap(ax(2),landcolor);
colorbar
```

However, how do you use multiple colormaps in one axes? That’s a little more tricky. Below, you find commented code that shows how to plot a colored and hillshaded DEM together with a stream network colored with ksn-values.

```% Some data (DEM and ksn of river network)
DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD  = FLOWobj(DEM);
A   = flowacc(FD);
S   = STREAMobj(FD,'minarea',1000);
k   = smooth(S,ksn(S,DEM,A),'K',1000);

% Create a new figure
hFig = figure;
% Call imageschs which has some parameter/value pairs
% to control colormap and colorbar appearance
imageschs(DEM,[],'colormap','landcolor',...
'colorbar',true,...
'colorbarylabel','Elevation');

% Get the current axes
h_ax = gca;
% Create a new axes in the same position as the
% first one, overlaid on top.
h_ax_stream = axes('position', get(h_ax, 'position'),...
'Color','none');
% Plot the stream network
plotc(S,k,'LineWidth',2);
% Make axis invisible
set(h_ax_stream,'Color','none');
% Create a colorbar
hc = colorbar(h_ax_stream,'Location','southoutside');
% ... and label it
hc.Label.String = 'ksn';
caxis(h_ax_stream,[0 200])
% Perfectly align both axes
set(h_ax_stream,'position', get(h_ax, 'position'),...
'XLim',xlim(h_ax),...
'YLim',ylim(h_ax),...
'DataAspectRatio',[1 1 1]...
);
% Link both axes so that you can zoom in and out

% Resizing the figure may screw it all. Using the figure
% property ResizeFcn makes sure both axes remain perfectly
% aligned
hFig.ResizeFcn = @(varargin) ...
set(h_ax_stream,'Position',get(h_ax, 'position'));
```

The trick is to use additional axes. The difficulty lies in perfectly aligning them and to make alignment robust against resizing the figure.

Overview of scientific colormaps

Posted on

Visualization is a vital part of scientific presentation and communication. This is why I have recently added some utilities that support making plots in TopoToolbox (see also my previous blogs “Jet is dead” and “Better colormaps with TopoToolbox”). One of the major additions has been Fabio Crameri’s scientific colormaps which are accessible through the function ttscm that you’ll find in the folder colormaps.

Fabio recently amended his compilation by a number of color schemes that I have now added. Please see an overview of the available colormaps in the figure below which I created by following lines of code.

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
cmaps = ttscm;
for r = 1:numel(cmaps);
subplot(6,4,r);
imageschs(DEM,[],'colormap', ttscm(cmaps{r}),...
'colorbar',false,'ticklabels','none');
title(cmaps{r});
end
```

Posted on

Jet is dead, long live parula! The criticism of MATLAB’s former default colormap is all around. Jet is perceptually non-uniform. Our brains tend to detect changes between particular kinds of colors better than the transition between others. The colormap jet is thus misguiding, it accentuates parts of your data and thus misrepresents it.

I have to admit that I used jet for quite a while. And it was the default colormap of imageschs, one of the core visualization tools in TopoToolbox. These days are over now. The default colormap is now parula which you’ll find in the folder colormaps. In addition, you’ll find other colormaps:

– landcolor
– flowcolor
– magmacolor
– and a whole bunch of terrain coloring maps in ttcmap

And the good thing: Now you can have multiple colormaps in one figure…

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
subplot(2,2,1);
imageschs(DEM,[],'colormap','landcolor','ticklabels','nice')
title('DEM - landcolor')
subplot(2,2,2);
imageschs(DEM,G,'colormap',flipud(magmacolor(255)),'ticklabels','nice')
subplot(2,2,3);
FD = FLOWobj(DEM);
A  = flowacc(FD);
A  = dilate(A,ones(5));
imageschs(DEM,sqrt(A),'colormap','flowcolor','ticklabels','nice')
title('Flow accumulation - flowcolor')
subplot(2,2,4)
[clr,lims] = ttcmap(DEM,'cmap','gmtglobe');
imageschs(DEM,[],'colormap',clr,'caxis',lims,'ticklabels','nice')
title('DEM - gmtglobe')
``` Of course, you can use all kinds of colormaps with imageschs such as hot, spring, etc. If you use custom colormaps, remember that

• imageschs requires colormaps with less than or equal 256 colors or 255 if there are nans in the data.
• you can easily flip colormaps, e.g. flipud(gray(255)).
• you can easily build your own colormaps using colormapeditor.

Do you have a favourite colormap that you want to share with others and that should be included in a future release of TopoToolbox. TopoToolbox could need a bipolar colormap for displaying curvature, for example.

Better maps with TopoToolbox

Posted on

I am increasingly using MATLAB to create maps. I simply find the workflow cumbersome that combines TopoToolbox, ArcGIS or QGIS. This is among the reasons why I implemented a new function in TopoToolbox that helps you coloring DEMs, particularly those that include topography and bathymetry. The function is ttcmap and you’ll find it in the colormaps folder.

Here is an example application of this function. I’ll use a DEM that covers one of the world’s steepest topographic gradients from the Andes down to the Peru-Chile deep-sea trench.

Before you run the function, choose a colormap. Calling ttcmap without input arguments displays a table that shows the different colormaps and their respective elevation ranges.

```ttcmap

Colormap_Name    Min_Elevation    Max_Elevation
_____________    _____________    _____________

'gmtrelief'       -8000           7000
'france'          -5000           4000
'mby'             -8000           4000
'gmtglobe'       -10000           9500

```

Now here are three colormaps that span the elevation range in our DEM. Note that ttcmap outputs the array zlimits which is required as input to imageschs so that the land-sea boundary is exactly at zero. ttcmap is similar to the mapping toolbox function demcmap but more precisely sets the land-sea transition to the zero value.

```[cmap,zlimits] = ttcmap(DEM,'cmap','gmtrelief');
imageschs(DEM,[],'caxis',zlimits,'colormap',cmap)
```
```[cmap,zlimits] = ttcmap(DEM,'cmap','mby');
imageschs(DEM,[],'caxis',zlimits,'colormap',cmap)
```
```[cmap,zlimits] = ttcmap(DEM,'cmap','gmtglobe');
imageschs(DEM,[],'caxis',zlimits,'colormap',cmap)
```

Tanaka contours – a MATLAB solution

Posted on Updated on

Today, I read about the Tanaka method, an approach to depicting relief. The method relies on filled contours and coloring contour lines according to their direction and light origin. I have seen solutions in QGIS (here and here) which made me thinking whether the Tanaka method can also be implemented in MATLAB.

One difficulty with Tanaka’s method is that coloring of contour lines requires knowledge about whether higher values are found to the left or the right of a contour line. GDAL’s contouring algorithm has consistent behavior in this respect. However, MATLAB’s algorithm doesn’t. My trick is to calculate surface normals and their angle to the light source, and then interpolate the gridded directions to the vertex locations of the contour lines.

Now here is my approach which uses some undocumented features of the HG2 graphics. An excellent resource for undocumented features in MATLAB is Yair Altman’s blog.

```function tanakacontour(DEM,nlevels)

%TANAKACONTOUR Relief depiction using Tanaka contours
%
% Syntax
%
%     tanakacontour(DEM,nlevels)
%
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 8. December, 2017

if nargin == 1;
nlevels = 10;
end

% calculate solar vector
azimuth = 315;
azid    = azimuth-90;
azisource = azid/180*pi;
[sx,sy,~] = sph2cart(azisource,0.1,1);

% get coordinate matrices and z values
[Z,X,Y] = GRIDobj2mat(DEM);

% calculate surface normals
[Nx,Ny] = surfnorm(Z);
N = hypot(Nx,Ny);
Nx = Nx./N;
Ny = Ny./N;

H = GRIDobj(DEM);
H.Z = Nx*sx + Ny*sy;

% Get coordinate matrices and Z
[~,h] = contourf(X,Y,Z,nlevels);
axis image
drawnow

% colormap
cmap = gray(255)*255;
cval = linspace(0,1,size(cmap,1))';

% the undocumented part
for r = 1:numel(h.EdgePrims)
% the EdgePrims property contains vertex data
x = h.EdgePrims(r).VertexData(1,:);
y = h.EdgePrims(r).VertexData(2,:);
% interpolate hillshade grayscale values to contour vertices
c = interp(H,double(x(:)),double(y(:)));
% interpolate to colormap
clr = interp1(cval,cmap,c);
% adjust color data to conform with EdgePrims ColorData property
clr = clr';
clr = uint8([clr;repmat(255,1,size(clr,2))]);
% set all properties at once
set(h.EdgePrims(r),'ColorBinding','interpolated','ColorData',clr,'LineWidth',1.2);
end
```

Ok, now let’s try this on some data. In this code, I have used GRIDobj, the class for gridded, georeferenced data.

```[X,Y,Z] = peaks(500);
tanakacontour(GRIDobj(X,Y,Z),10)
``` And now for the Big Tujunga DEM:

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
tanakacontour(DEM,5)
``` Now that looks ok, but I am not totally convinced by these unrealistic terraces, to be honest. Moreover, my approach so far doesn’t include changing line widths. Contour lines should actually be thicker for slopes that face both towards and away from the light source. However, I didn’t figure out how to do this so far. Any ideas?

An alternative to along-river swath profiles

Posted on Updated on

Today, I want to show another application of TopoToolbox’s new smoothing suite that Dirk and I have recently released together with our discussion paper in ESURF. I will demonstrate that using different methods associated with STREAMobj enable you to create plots that are quite similar to swath profiles calculated along rivers. Let’s see how it works.

Swath profiles require a straight or curved path defined by a number of nodes. Each node has a certain orientation defined by one or two of its neighboring nodes. Creating swath profiles then involves mapping values from locations that are perpendicular to that orientation and within a specified distance. We will use a simplified version of mapping values lateral to the stream network which is implemented in the function maplateral. maplateral uses the function bwdist that returns a distance transform from all pixels (the target pixels) to a number of seed pixels in a binary image. In our case, seed pixels are the pixels of the stream network and bwdist identifies the  target pixels from which it calculates the euclidean distance to the stream pixels. Usually, there are several target pixels for each seed pixel so that we require an aggregation function that calculates a scalar based on the vector of values that are associated with each seed pixel. Unlike swath profiles, however, this approach entails that some seed pixels may not have target pixels  up to the maximum distance. We can thus be pretty sure that some pixels may be missing some of the information that we want to extract. While we cannot avoid this problem using this approach, we can at least reduce its impact using the nonparametric quantile regression implemented in the crs function.

Here is the approach. We’ll load a DEM, derive flow directions and a stream network which is in this case just one trunk river. We use the function maplateral to map elevations in a distance of 2 km to the stream network. Our mapping uses the maximum function to aggregate the values. Hence, we will have for each pixel along the stream network the maximum elevation in its 2 km nearest-neighborhood. The swath is plotted using the function imageschs and the second output of maplateral. Then we plot the maximum values along with a profile of the stream network.

```DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1e6,'unit','map');
S = klargestconncomps(S);
S = trunk(S);

z = getnal(S,DEM);

subplot(2,1,1)
'colorbar',false,...
'ticklabels','nice');

subplot(2,1,2)
plotdz(S,z)
hold on
plotdz(S,zmax)
legend('River profile','maximum heights in 2 km distance')
``` This looks quite ok, but the topographic profile has a lot of scatter and a large number of values seem not to be representative. Most likely, those are the pixels whose closest pixels fail to extend to the maximum distance of 2 km. Rather, I’d expect a profile that runs along the maximum values of the zigzag line. How can we obtain this line? Well, the crs function could handle this. Though it was originally intended to be applied to longitudinal river profiles, it can also handle other data measured or calculated along stream networks. The only thing we need to “switch off” is the downstream minimum gradient that the function uses by default to create smooth profiles that monotoneously decrease in downstream direction. This is simply done by setting the ‘mingradient’-option to nan. I use a smoothing parameter value K=5 and set tau=0.99 which forces the profile to run along the 99th percentile of the data. You can experiment with different values of K and tau, if you like.

```figure
plotdz(S,z);
hold on
hold off
``` This looks much better and gives an impression on the steepness of the terrain adjacent to the stream network. Let’s finally compare this to a SWATHobj as obtained by the function STREAMobj2SWATHobj.

```figure
SW = STREAMobj2SWATHobj(S,DEM,'width',4000);
plotdz(SW)
xlabel('Distance upstream [m]');
ylabel('Elevation [m]');
``` Our previous results should be similar to the maximum line shown in the SWATHobj derived profile. And I think they are. In fact, the SWATHobj derived profiles also shows some scatter that is likely due to the sharp changes of orientation of the swath centerline.  While we can remove some of this scatter by smoothing the SWATHobj centerline, I think that the combination of maplateral and the CRS algorithm provides a convenient approach to sketch along-river swath profiles.

Orientation of stream segments

Posted on Updated on

Today’s guest post comes from Philippe Steer, assistant professor at the University of Rennes. Philippe contributed two new functions and I asked him to describe them here to show potential applications. Thanks, Philippe, for this contribution!

I recently needed to compute river network orientation. Indeed river network orientation can reveal distributed tectonics deformation at the vicinity of strike-slip faults (e.g. Castelltort et al., 2012) or can enable to identify developing vs paleo river networks.  One way to obtain orientation is first to extract the river segments that are linking confluences to outlets or confluences to channel heads and then to compute the orientation of the resulting segments. The streampoi function allows to extract these feature points of the network such as channel heads, confluences, b-confluences (the points just upstream any confluence) and outlets. To connect these feature points into segment, one can use the propriety that each segment belongs to a unique watershed that is given by the function drainagebasins using as outlets the true river outlets and the b-confluences. By definition, each segment also corresponds to a unique value of strahler order.

Once the river segments are identified, segment orientation is easily obtained. What is interesting is that defining these segments also paves the way to compute several other geometric metrics averaged over the river segment such as segment direct or flowpath length, sinuosity, direct or flowpath slope, mean drainage area and even steepness index (if assuming a value of concavity).

To cleanly do the job, I have written two MATLAB functions:

Here is a piece of code showing you the use of these functions (assuming you have already loaded the DEM of interest):

```%% Flow direction and drainage area
DEMf = fillsinks(DEM); FD= FLOWobj(DEMf); A = flowacc(FD);

%% Extract rivers
% Set the critical drainage area
A_river=1e7/(res^2); % m^2
W = A>A_river;
% Extract the wanted network
S = STREAMobj(FD,W);

%% River segments
% Extract segment and compute their orientations
mnratio=0.5;
[segment]=networksegment(S,FD,DEM,A,mnratio);
% Plot results
plotsegmentgeometry(S,segment)
```

To illustrate the results, I thought the South Island of New-Zealand was a good example. Indeed, as first suggested by Castelltort et al. (2012) the river catchments at the south of the Alpine fault are tilted, possibly suggesting that they have recorded a distributed shear deformation. This tilting is also quite obvious when looking at the river segments in the following figure. Another interesting feature is that the river segments north of the Alpine fault, facing the tilted rivers, display a higher steepness than their south neighbours. This could be due to a long list of forcings, including active tectonics, orographic precipitation gradient or the transition from glacial to fluvial erosion.

References:
Castelltort, S et al. (2012). River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain. Nature Geoscience, 5(10), 744-748. [DOI: 10.1038/ngeo1582]