DEM to volumetiric

converts a Digital Elevation Model (or 2.5D data) to a volumetric raster (or true 3D data)

Description

In a DEM height information is strored as h(x,y). This scripts converts h(x,y) to I(x,y,z) based on

(a) if z < h(x,y) then I(x,y,z) = 1
(b I(x,y,z) = 0 otherwise.

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
boxcount

• MATLAB release: MATLAB 7.7 (R2008b)

[Wolfgang Schwanghart]

While this function has a reasonable help, it lacks various features that I'd wish to see in the future to make this function valuable for others. First, provide a h1 line. It's the first line in your help block and can be found by lookfor. Currently, your h1 line is empty and thus the function can never be found.

You allocate a 3d-array of zeros which you subsequently fill with ones in 3 nested while loops. If there are only ones and zeros, why don't you use a logical array. It requires much less memory than using double precision.

In line 9 you than call the function zeros in following manner:
XYZ = zeros(size(XY,1), size(XY,2), max(max(XY)));

This may produce a warning (actually I was wondering that it not even returns an error) for DEMs with a maximum elevation that is not an integer value.
Warning: Size vector should be a row vector with integer elements.

As such, you implicitly assume that everyone uses integer representations of DEMs, that may often be not the case. Moreover, wouln't it be better to give the user more control on how your function handles pixel to voxel conversion? I think, you should allow the user to define a base level (now it will always be the sea level or another fixed reference value) and vertical resolution of the generated 3d-array. Moreover, I'd like to see that the function returns a coordinate system for the 3d array (for example by calling meshgrid(x,y,z) or so.

Currently the function is heavily looped and I am sure there are ways to make it computationally more efficient.