DEM raster to 3D PLY format

Since I originally posted this 3D printing conversion method, several other options are now available to convert a raster DEM/DTM (digital elevation/terrain model) to a 3D format. Here are just a few worth checking out before I list a workflow using a GDAL/Python script method as posted by Jake from Zurich.

• Blender (free) has a plug-in to open up HiRISE DTMs. See instructions here. From Blender, I think you can export to STL or OBJ? http://www.uahirise.org/dtm/howto.php
• SpaceBlenders - We also support a variation on this tool but it requires getting your hands dirty with some configurations to get the GDAL library installed with Blender: https://github.com/USGS-Astrogeology/SpaceBlender (based on a NAU capstone) 
• GRASS: https://grass.osgeo.org/screenshots/3D/
• Unity: https://alastaira.wordpress.com/2013/11/12/importing-dem-terrain-heightmaps-for-unity-using-gdal/
• new to me (2019). AMES Stereo Pipeline (ASP) using their point2mesh routine. It is always good to dig through their routines for gems. https://ti.arc.nasa.gov/tech/asr/groups/intelligent-robotics/ngt/stereo/

Since I am a GDAL user, I usually first dig around for a method using it (since it has support for several planetary formats). When I ran across this script by Jake it was just the simple solution I was after. I wanted GDAL support in a command-line version with no fancy interface or setup. I would let Jake's original post stand on its own but there are a couple steps a I had to take to get it working for me.

1. It currently doesn't support NoData values. It supports them but it doesn't strip them. This seems like a tangible thing to handle in the code but I haven't tried so below is a work-around.
2. I assume you will want to be careful on the size of the input DEM. Thus here, for my Mars HiRISE example (DTEEC_025123_2045_026811_2045_A01.IMG in PDS v3 format), I sub-sample it to 5m/p. This DEM is located in Mawrth Valles. Image credit: NASA/JPL/University of Arizona/USGS.

NoData and Resampling Workflow

First get yourself a GDAL/Python environment running. I recommend Anaconda Python. Once Anaconda is installed type "conda install gdal". Done.

Second head over to gis.stackexchange and copy Jake's gdal_rastertotrn.py code. I find that filename hard to read "*totrn.py" (where trn must mean terrain). I would prefer simply the script name "gdal2PLY.py" (which is now also on github). Okay now to deal with NoData and the large HiRISE DTM listed above.

          1.) find minimum Z value

$ gdalinfo -mm DTEEC_025123_2045_026811_2045_A01.IMG
Computed Min/Max=-2790.594,-2063.563
NoData Value=-3.4028226550889045e+038  

Most 3D apps only want to ingest positive integer values. (e.g. 16 unsigned PNG). This file will kill most 3D apps: it is 32bit floating point (first strike), has negative elevation range (2nd strike), and look at that giant NoData value (3rd strike). Fortunately the only thing we need to account for is the NoData value. The negative value elevation range or floating point type will not be an issue.

2.) map the DEMs Nodata to a nice round value below the minimum (for our example say -2791). Note during this step we can also sub-sample the file to a lower resolution. We could create a giant PLY but I doubt many apps will be able to open it.

$ gdalwarp -dstnodata -2791 -tr 5 5 -r bilinear DTEEC_025123_2045_026811_2045_A01.IMG temp_DTEEC_5m_nodata.tif

where -dstnodata = destination Nodata value of -2791
where -tr 5 5 means target resolution of 5m in X, Y
where -r bilinear means the resampling method (bilinear should be smooth enough)

3. ) Now simply run Jake's script in your GDAL/Python environment.  

    $ python gdal2PLY.py temp_DTEEC_5m_nodata.tif out_DTEEC_025123_2045_026811_2045_A01.ply

         4.) Load in Meshlab to visualize and/or convert to other formats (like STL or OBJ).

Figure: resultant binary PLY formatted surface in meshlab. I added some "radiance scaling" under Render - Shaders to highlight the topography breaks. Note this will still have a NoData plane but since it is near the minimum value it doesn't mess-up the exaggeration.

Future: Still need a better method to handle the NoData...