GRIDobj

Fill missing values

Posted on Updated on

DEMs often have missing values. Reasons include lack of coherence in radar imagery or lack of backscattered signals in LiDAR data. Sometimes, it is just a pixel that is missing, but frequently there are larger contiguous regions with no data. These missing values can be annoying because they inhibit flow routing.

An approach to get rid of these pixels is to fill them with elevation data obtained by another DEM. For example, missing data in some of the first versions of the SRTM were often filled with 1 km resolution GTOPO3 data. Alternatively, one can use techniques that estimate pixel values from the surrounding pixels that have data. Those of you, who have worked with ArcGIS might know the function nibble, which uses nearest neighbor interpolation to estimate missing values.

TopoToolbox also has a function to fill regions with missing values. The function name is inpaintnans because inpainting is a common technique in image processing and restoration. The function takes either another DEM to fill missing values, or uses different techniques to estimate from the surrounding pixels. The option ‘nearest’ is the one that is equal to the tool nibble in ArcGIS.

However, while nearest neighbor interpolation is suited for nominal or ordinal scaled data, it might produce unwanted steps in continuous data such as elevation. Thus, the recommended (and default) approach in TopoToolbox is laplacian interpolation. This type of interpolation is explained nicely by Steve Eddins here and implemented in the function regionfill that is used by the TopoToolbox function inpaintnans. In short, this approach interpolates inward by solving a Dirichlet boundary value problem and thus derives a smooth surface.

The Taiwan dataset that is among the new example files on the DEM repository has some missing values. Here is how to apply inpainting:

DEM = readexample('taiwan');
axislim = [0.1970 0.2384 2.5183 2.5951]*10^6;

subplot(2,2,1);
imageschs(DEM)
hold on
plot([axislim(1) axislim(2) axislim(2) axislim(1) axislim(1)],...
[axislim(3) axislim(3) axislim(4) axislim(4) axislim(3)],...
'r','LineWidth',2)

subplot(2,2,2);
imageschs(DEM,[],'caxis',[0100])
axis(axislim);

subplot(2,2,3);
DEMin = inpaintnans(DEM);
imageschs(DEMin);
hold on
plot([axislim(1) axislim(2) axislim(2) axislim(1) axislim(1)],...
[axislim(3) axislim(3) axislim(4) axislim(4) axislim(3)],...
'r','LineWidth',2)

subplot(2,2,4);
imageschs(DEMin,[],'caxis',[0 100]);
axis(axislim);
The SRTM-3 of Taiwan. Top left: There are a regions with missing values. Top right: Zoom to this region. Bottom left: Missing values filled using inpaintnans. Bottom right: Zoom to region with DEM values filled.

It seems that some regions were not inpainted. However, these regions have a connection to the sea. To be filled, regions with no data must be entirely surrounded by pixels with data.

Inpainting can deal with missing values in DEMs. However, there are also other applications. In our recent paper by Amelie Stolle we use laplacian interpolation to reconstruct a model of the undissicted surface of Pokhara Formation in Nepal which allows us to estimate incision depths and rates.

Reference

Stolle, A., Schwanghart, W., Andermann, C., Bernhardt, A., Fort, M., Jansen, J.D., Wittmann, H., Merchel, S., Rugel, G., Adhikari, B.R., Korup, O., 2019. Protracted river response to medieval earthquakes. Earth Surface Processes and Landforms, 44, 331-341. [DOI: 10.1002/esp.4517]

Advertisements

mappingapp – a tiny tool to map points along river networks

Posted on Updated on

TopoToolbox has a new tool. Ok, it’s a small tool with limited functionality, but it might be helpful if your aim is to map points along river networks. “Mapping points can be easily done in any GIS!”, you might think. True, but as I am currently working on knickpoints in river profiles (see also the automated knickpointfinder), I wrote this tool to quickly map knickpoints both in planform and profile view of the river network.

So, here is mappingapp. Give it a try.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1000);
mappingapp(DEM,S)
DrZMINIW4AEMnH6.jpg large.jpg
mapping-app. A tiny tool to map knickpoints along river profiles.

Currently, the GUI is limited to a few basic tools. You can map a single point which automatically snaps to the river network S. The zoom tools allow you to navigate the DEM. You can add a new point using the + button. This will make the previous point permanent and add a new row to the table. Finally, the table can be exported to the workspace. The table contains the coordinates and z-values as well as a column IXgrid. IXgrid contains the linear indices into the DEM.

Controls of mappingapp

Jet is dead

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);
G = gradient8(DEM);
imageschs(DEM,G,'colormap',flipud(magmacolor(255)),'ticklabels','nice')
title('Gradient - magmacolor')
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')

 

colormaps.png

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)

gmtrelief.png
gmtrelief

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

mby.png
mby

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

gmtglobe.png
gmtglobe

I downloaded the colormaps from cpt-city, a great resource for colormaps.

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)
%
% See also: contourf, hillshade
%
% 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)

tanaka1.png

And now for the Big Tujunga DEM:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
tanakacontour(DEM,5)

tanaka2.png
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?

One tool to download them all

Posted on

Some changes and additions to TopoToolbox are subtle and may not be visible at first glance. Hence, I report them here. One such addition is part of version 2.2 (which Dirk and I are about to finalize shortly) and is found in the folder IOtools (input/output tools): readopentopo. The function uses opentopography’s RESTful Web service for global raster data access for seamless download of global digital elevation models. readopentopo offers access to SRTM-3, SRTM-1, and ALOS 3D World data at your finger tips. By default, the downloaded *.tif-file is stored in your system’s temporary folder and is deleted as soon as it is read into MATLAB. If you want to permanently save the file on your hard drive make sure to define a file path and to set the ‘deletefile’-option to false. The tif-file comes in geographic coordinates. Hence, don’t forget to reproject the DEM to UTM (see reproject2utm) before you start working with it.

Many thanks to members of the opentopography-team who make the data and service available. Please use these services reasonably and don’t excessively download data!

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]