DEM grid size specification and learning a bit of TopoToolbox

I have been wanting to learn TopoToolbox for a while, and last week I had a chance to get started. Spending time at the Universität Potsdam, Institut für Erd- und Umweltwissenschaften, I have met with Wolfgang Schwanghart and Dirk Scherler occasionally. They are coauthors of TopoToolbox.

My ASU colleagues Adam Forte and Kelin Whipple have been using TopoToolbox increasingly over the last few years. Recently Kelin had a question as to why some DEMs grids were causing an error in TopoToolbox and why some were not. I dug into it and here is what I found. Projected DEMs should have EXACTLY the same x and y cell sizes and expressed in precise values (nearest meter or 10th or 100th of a meter).

DEMs (such as SRTM 30 meter) can be downloaded as geotiff formats and with geographic coordinate systems from sources such as OpenTopography: select Global Data tab.

For example, I have chosen a piece of data from Java along the Cimandiri fault (see Marliyani, et al., 2016 for example).

To be useful for most geomorphic analyses, the data should be projected to UTM so that the horizontal and vertical units are the same. Often, the projection is done in ArcGIS. And, usually it seems fine to let Arc determine the cell size automatically.

Notice how ArcGIS decided that the resolution for the cells should be: 30.8462728281739.

Here I set the cell size to exactly 30 m

Turning to TopoToolbox, as we compute drainage network properties, including contributing area, there is a check that the cell sizes are the same (inside the GRIDobj.m function):

if abs(abs(dx)-abs(dy))>1e-9;
    error('TopoToolbox:GRIDobj',...
    'The resolution in x- and y-direction must be the same');
    end

The 1e-9 is a somewhat arbitrarily small number which one would think would not cause a problem.

I looked into this problem tracking along with the TopoToolbox processing and some of the MATLAB built in tools.

1. Read the geotiff using MATLAB's geotiffread:

   [demdata, R] = geotiffread(filename);
   geotiffreaddifference = R.CellExtentInWorldX-R.CellExtentInWorldY;
   

   The R object has a number of values including the cell size of the geotiff:
   

   R.CellExtentInWorldX= 30.79776426908749800
   R.CellExtentInWorldY= 30.79776426908791000
   R.CellExtentInWorldX-R.CellExtentInWorldY= -4.12115e-13
   

2. Read the geotiff into the DEM object in TopoToolbox and do a similar check as above (refmat is similar to R):
   

   DEM = GRIDobj(filename);
   refmatdifference = abs(DEM.refmat(2,1))-abs(DEM.refmat(1,2));
   

   Again we see the same very small grid size difference:
   

   abs(DEM.refmat(2,1))-abs(DEM.refmat(1,2))= -4.12115e-13
   

3. The problem comes in a calculation that builds out the x and y vectors that includes a cumulative multiplication from the grid sizes inside a TopoToolbox function refmat2XY:
   

   nrrows = siz(1);
   nrcols = siz(2);
   
   x = [ones(nrcols,1) (1:nrcols)' ones(nrcols,1)]*R;
   x = x(:,1)';
   
   y = [(1:nrrows)' ones(nrrows,2)]*R;
   y = y(:,2);
   

   R is the refmat object. But, here we get the error:
   

   dx=x(1)-x(2)= -30.79776426916942000
   dy=y(1)-y(2)= 30.79776426777243600
   abs(dx)-abs(dy)= 1.39698e-09
   

   Which then would fail the GRIDobj test.
If I run the same set of checks on the file I projected with exactly 30 m grid cell size (recall above), I don't get the problem:


geotiffread:
R.CellExtentInWorldX= 30.00000000000000000
R.CellExtentInWorldY= 30.00000000000000000
R.CellExtentInWorldX-R.CellExtentInWorldY= 0
refmat difference as produced from GRIDobj:
abs(DEM.refmat(2,1))-abs(DEM.refmat(1,2))= 0
Difference after GRIDobj2mat:
dx=x(1)-x(2)= -30.00000000000000000
dy=y(1)-y(2)= 30.00000000000000000
abs(dx)-abs(dy)= 0

I think what is happening is a bit of roundoff error (see this link for detailed discussion: What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg): "Therefore the result of a floating-point calculation must often be rounded in order to fit back into its finite representation." MATLAB does all of its calculations by default in double precision, so we should have access to 16 decimal digits (e.g., https://en.wikipedia.org/wiki/IEEE_754-1985). So, why we loose precision up to 10^-9 at times is a little surprising to me, but it can obviously happen.

Thus there are a couple of workarounds:

1. Comment out the dx and dy check in GRIDobj:

   %if abs(abs(dx)-abs(dy))>1e-9;
   %    error('TopoToolbox:GRIDobj',...
   %    'The resolution in x- and y-direction must be the same');
   %    end
   

   Or make the threshold larger.
2. Explicitly set the grid resolution to a round number when projecting in ArcGIS or other software.
3. Use the reproject2utm.m command to project dems before using other TopoToolbox operations. It sets the resolution to be exactly the same for x and y

Here is a script to run these calculations: DEM_cell_size_Script.m

Thanks to Chris Crosby, Benjamin Gross, Wolfgang Schwanghart, Kelin Whipple, and Mike Zoldak for discussions.