Since I originally posted my upslope toolbox to MATLAB Central back in August, I have heard from some experts about an issue related to NaNs in the DEM data. Specifically, some data sets record elevation only for DEM pixels within a particular catchment basin. Pixels outside the catchment basin typically are NaN. Some of my code doesn't handle NaN values in the DEM.
After some e-mail correspondence, I settled on this approach:
- In the pixel flow direction computation, replace groups of NaN pixels that touch the border by a high value. (You can read about the technique for find such pixels in this post.) This replacement makes the flow direction on the watershed pixels point inward, toward the catchment basin.
- When solving the flow matrix for upslope area, use 0 instead 1 on the right-hand side for exterior NaN pixels. This means that such pixels don't contribute any flow to the system.
Here's a simple example:
[x,y] = meshgrid(-3:3); E = hypot(x,y); E(abs(x) + abs(y) > 3) = NaN
E =
NaN NaN NaN 3.0000 NaN NaN NaN
NaN NaN 2.2361 2.0000 2.2361 NaN NaN
NaN 2.2361 1.4142 1.0000 1.4142 2.2361 NaN
3.0000 2.0000 1.0000 0 1.0000 2.0000 3.0000
NaN 2.2361 1.4142 1.0000 1.4142 2.2361 NaN
NaN NaN 2.2361 2.0000 2.2361 NaN NaN
NaN NaN NaN 3.0000 NaN NaN NaN
R = dem_flow(E); T = flow_matrix(E, R); A = upslope_area(E, T)
A =
0 0 0 1.0000 0 0 0
0 0 1.0000 2.0000 1.0000 0 0
0 1.0000 1.8112 4.1888 1.8112 1.0000 0
1.0000 2.0000 4.1888 25.0000 4.1888 2.0000 1.0000
0 1.0000 1.8112 4.1888 1.8112 1.0000 0
0 0 1.0000 2.0000 1.0000 0 0
0 0 0 1.0000 0 0 0
Notice that the upslope area at the lowest pixel is 25. This is the number of non-NaN pixels in the DEM.
I uploaded version 1.2 of the upslope toolbox, with these NaN-related changes, to the MATLAB Central File Exchange in early October. If you have data sets that use NaNs outside the catchment basin of interest, you might want to download the new version.
Get
the MATLAB code
Published with MATLAB® 7.5
Hi Steve.
Thanks for the great upslope area blog. I hope, geospatial analysis will always be an integral part in Matlab. In case you need another idea for a geo-related blog category, what about one about TOPOGRID, a method to grid elevation data from irregularly distributed elevation data in a hydrologically correct way. You’ll find out more about it in Hutchinson, M. F. (1989), J. of. Hydrology, 211-232. Like the upslope topic, it could be nicely handled in combination with the image processing toolbox.
Best regards,
Wolfgang
Wolfgang—Thanks for the suggestion!