I always try to continue working on updates of TopoToolbox. Updates, however, will mainly comprise tools that I need for my work. If you are developing tools by yourself and you are willing to share them with others, please don’t hesitate to contribute to TopoToolbox. How does it work? Just send me an email with your function, I’ll review it and include it in the next update.

Functions should fulfill following criteria:

  • No scripts! If you have a good script that explains a type of analysis you are doing, we could consider including it in the documentation section.
  • The function should have a proper function name that doesn’t clash with built in functions of Matlab and common toolboxes.
  • H1-line, so that the function can be found by lookfor.
  • A comprehensive documentation structured as other functions in the toolbox.
  • The code should be commented as much as possible.
  • The code should be optimized using vectorization techniques and execute fast and reliably.
  • At best, the code should only depend on the Image Processing Toolbox. If other products are required, please notify.

Do you think TopoToolbox is missing important tools, but you are not willing or cannot programme them by yourself, just send me a request. Either I’ll try to work it out by myself or I’ll post it here hoping someone else will do it.


20 thoughts on “Contribute

    Jonathan Phillips said:
    January 12, 2016 at 11:43 am

    Howdy! I appreciate your including my blog on your Blogroll. Don’t ask me why I clicked on my own blog, but I discovered you’ve got my old URL. The new one is I don’t know why it had to be changed; it is difficult for me to get a straight answer from our IT people.



      wschwanghart responded:
      January 15, 2016 at 8:41 am

      Thanks for the note, Jonathan. Corrected. Cheers, Wolfgang

    Denis said:
    May 8, 2017 at 12:00 pm

    Hi Wolfgang,

    Thank you for this tool.
    After having used the ‘fillsinks’ filling option on my DEM, I got the following error when either trying to crop the filled DEM, or when e.g. subtracting the original DEM from it (I would like to have an idea about the residuals):

    Error using GRIDobj (line 169)
    The resolution in x- and y-direction must be the same

    Error in GRIDobj/crop (line 176)
    DEMc = GRIDobj(x,y,demc);

    My question is: is the ‘fillsinks’ option altering the resolution differently along X and Y?


      wschwanghart responded:
      May 8, 2017 at 12:20 pm

      Hi Denis, thanks for your comment. This is not a problem with the function fillsinks, but with the function crop. The problem has been noted before (see here: ) and is due to very small differences in cell sizes in x and y direction. These differences don’t play a role at first. Yet, after recalculating the georeferencing after cropping these differences blow up due to some numerical issues. I will need to fix this as soon as I have time. Until then, check the work arounds in above link.
      Best, Wolfgang

        wschwanghart responded:
        May 8, 2017 at 2:16 pm

        Note that I have updated the function GRIDobj/crop and the problem should not persist any longer in the latest version available on github (

        Denis said:
        May 9, 2017 at 8:24 am

        Many thanks Wolfgang,
        The ‘crop’ function now works perfectly – no resolution conflicts with the parent matrix anymore.

        Denis said:
        May 9, 2017 at 9:11 am

        I have another question regarding the ‘crop’ function.
        I would like to re-use the mask obtained from the ‘interactive’ option in order to crop several DEMs according to the exactly same extent and location. How can I do this?
        I looped into the ‘crop.m’ script, and it seems that there should be a way using the ‘sub2coord’ subfunction.

    wschwanghart responded:
    May 9, 2017 at 9:28 am

    Hi Denis, in fact crop has some undocumented second output argument, which is the MASK.
    [DEMc,MASK] = crop(DEM,…);
    In case you have different GRIDobjs A, B, C… with exactly the same extent as DEM folllowing code should work
    Ac = crop(A,MASK);
    Bc = crop(B,MASK);

    Hope this works,

      Denis said:
      May 9, 2017 at 9:37 am

      Thanks again!
      This is a very nice tool.

    Shantamoy said:
    October 13, 2017 at 1:56 pm

    Hi Wolfgang,
    I was wondering that if we can get maximum gradient and corresponding relief values for the definite grids. It might be a slight variant of the EXCESSTOPOGRAPHY function. Is there any way to get these two variables with the existing functions in TopoToolbox?
    Thank you in advance.

      wschwanghart responded:
      October 13, 2017 at 3:15 pm

      Hi Shantomoy, thanks for your comment. I am not sure whether I understand correctly what you mean by maximum gradient and corresponding relief. Do you mean the maximum gradient measured for a specific pixel in the grid, and the relief (e.g. the topographic range within a certain radius) at this specific pixel location?
      Cheers, Wolfgang

        Shantamoy said:
        October 13, 2017 at 7:56 pm

        Hello Wolfgang,
        Many thanks for your reply and sorry for not being clear. Let me put it in this way.
        Say, I have divided the DEM into several grids (let the size of the square grid is 5 km). This is not like traditional focal statistics, rather all the grids have to be discrete and will have no overlap. Now suppose we have divided the whole DEM into 100 grids. (So the size of DEM is 50 km by 50km). Therefore, we have 100 relief values from each grid.
        The second part is the maximum slope. What I mean by maximum slope is one slope value from each grid. The relief is the difference between highest and lowest elevation value. Now if the locations of these highest and lowest points can be traced then the distance can be calculated. So
        Max_slope = tan_inverse(Relief/Distance).
        One problem will arise when there will be more than one max or min elevation value in one grid. But that can be solved by some other method. (Note that I have changed word gradient into slope)

    wschwanghart responded:
    October 13, 2017 at 10:40 pm

    Hi Shantamoy, ok, thanks for clarifying. Now that is what I’d refer to as block statistics. So far there is no function to do this, and I thank you for bringing this topic to my attention. I will implement block statistics definitely soon.
    So how could this issue be solved. I’ll start with calculating the topographic range to show how to approach the problem.

    DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
    R = GRIDobj(DEM);
    R.Z = blockproc(DEM.Z,[50 50],@(x) repmat(max(,x.blockSize));

    The function blockproc is probably best to do calculations on individual blocks. Here, it tiles the DEM in equal blocks of 50 by 50 pixels. You need to change the vector [50 50] to do the calculation on different block sizes. In addition, blockproc takes an anonymous function that determines what to do with the block. Note, that x is a structure array that is described at the bottom of the function’s documentation. Here I simply calculate the difference between the maximum and the minimum, and then replicate the elements to obtain the original block size as output. So, basically, this is a one-liner.

    For calculating the maximum gradient, we need to define another function that is a bit more bulky. Hence, we have to write a function that resides on MATLAB’s search path and that takes the a block structure as input. The function below is a first try and doesn’t solve for the issue of replicate maxima or minima.

    function g = maxgradient(blockstruct)
    % calculate minimum and maximum and their locations
    [zmax,ixmax] = max(;
    [zmin,ixmin] = min(;
    % calculate row and column indices of maxima and minima
    [rmax,cmax] = ind2sub(blockstruct.blockSize,ixmax);
    [rmin,cmin] = ind2sub(blockstruct.blockSize,ixmin);
    % calculate the distance
    d = hypot(cmax-cmin,rmax-rmin);
    % calculate the gradient
    g = (zmax-zmin)/d;
    % replicate elements
    g = repmat(g,blockstruct.blockSize);

    Now you can again call the function blockproc with following command.

    R.Z = blockproc(DEM.Z,[50 50],@maxgradient);

    Finally, this gives the gradient for a unit pixel spacing. Hence, calculate R as

    R = R/R.cellsize;

    And plot it:


    Hope this does the trick.

    Best regards, Wolfgang

      Shantamoy said:
      October 14, 2017 at 8:53 am

      Hello Wolfgang,
      That worked really nice. Thank you very much for making the things easier. However, there are still some pitfalls that I want to address. The maxgradient function needs to take care of the duplicate data. There might be more than one pixel which carries minimum or maximum value within a grid. Say for instance, if there is a flat plateau and/or a flat valley bottom within a single grid. It means that there might be multiple paths that connect max and min values (as there will be multiple values of max and min). So the minimum distance can be taken for the calculation of slope. Hope this will be taken care in the future.
      This tool is a great contribution towards the field of quantitative geomorphology. Thank you for the constant effort.

        wschwanghart responded:
        October 14, 2017 at 11:47 am

        Hi Shantamoy, for duplicate data, here is a modified function maxgradient. The trick is to use the function bwdist.

        function g = maxgradient(blockstruct)
        % calculate minimum and maximum and their locations
        zmax = max(;
        zmin = min(;
        % zero gradients, no further action required
        if zmax == zmin
            g = zeros(blockstruct.blockSize);
        IMAX = == zmax;
        IMIN = == zmin;
        D = bwdist(IMAX,'euclidean');
        d = min(D(IMIN));
        % calculate the gradient
        g = (zmax-zmin)/d;
        % replicate elements
        g = repmat(g,blockstruct.blockSize);

        Cheers, Wolfgang

    Richard Ott said:
    May 23, 2018 at 12:47 pm

    the crop function currently only is able to crop out rectangles. Would it be easy to implement a polygon cropping into the function? I looked into it but wasn’t able to modify the function and decided to use a slow workaround. I think such an option would be very useful (at least for me).

    wschwanghart responded:
    May 23, 2018 at 2:10 pm

    Hi Richard, you are right. crop currently supports clipping rectangles. However, if you supply a logical GRIDobj as second input argument, you can choose nan as fillval. Hence, my workaround would be according to following example:

    DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
    FD = FLOWobj(DEM);
    S  = STREAMobj(FD,'minarea',1000);
    S  = klargestconncomps(S);
    D  = drainagebasins(FD,S);
    MS = GRIDobj2polygon(D);
    I  = polygon2GRIDobj(DEM,MS);
    DEMc = crop(DEM,I,nan);

    Clearly, the step from D -> MS -> I is not necessary here but should illustrate how to combine polygon2GRIDobj and crop.

    What was your workaround?

    Cheers, Wolfgang

    betsy said:
    June 14, 2019 at 5:54 pm

    Hi Wolfgang, I sent you a message through mathworks but thought I should also post here. Thanks for sharing your codes. Why do X and Y have to be the same size in your quantile regression code?

      wschwanghart responded:
      June 15, 2019 at 10:03 am

      Hi Betsy, this is because y is a function of x. Hence, each observation must have a value of the independent and dependent variable.
      Best, Wolfgang

        Betsy said:
        June 15, 2019 at 9:11 pm

        Right, so the lengths have to be the same, but I was thinking you should be able to use multiple explanatory variables. Thanks for the quick reply.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.