Thursday, June 8, 2017

Using GDAL and/or NASA ASP to create a colorize hillshade

There are many good posts out there to help create colorshaded reliefs. Here I just want a quick method to create a "data" ready colorize hillshade. This means not a figure but a geospatial ready colorshade (e.g. GeoTiff) to use in a GIS or WMS server. And I need the capability to process very large GeoTiffs ( > 4GB).

Get the Tools

I recommend GDAL or NASA Ames Stereo Pipeline (ASP). Now ASP v2.6.0 actually contains the GDAL commands for this tutorial or you can use their own versions.

For GDAL binaries I recommend the Anaconda Python Environment (for all OSs). I am still digging out of Python 2.7 into Python 3.x but GDAL should work in either. First install the Anaconda environment and then run "conda install gdal" to get the gdal binaries. .

If you just want ASP there are binaries for Mac and Linux:

Get some Data

Raster digital elevation models (DEMs / DTMs) are getting pretty easy to find. And when needed you can even create your own in ASP for some instruments. For a simple example, lets grab a small 8MB HiRISE 32bit DTM from UofA:


Let's check out the image metadata
ast{91}> gdalinfo -mm DTEEC_001414_1780_001612_1780_U01.IMG
Driver: PDS/NASA Planetary Data System
Files: DTEEC_001414_1780_001612_1780_U01.IMG
Size is 1279, 1694
Coordinate System is:
Origin = (-401.706464281687772,-120698.120618077693507)
Pixel Size = (1.011855073757400,-1.011855073757400)
Corner Coordinates:
Upper Left  (    -401.706, -120698.121) (  5d30'24.40"W,  2d 2'10.50"S)
Lower Left  (    -401.706, -122412.203) (  5d30'24.40"W,  2d 3'54.60"S)
Upper Right (     892.456, -120698.121) (  5d29' 5.80"W,  2d 2'10.50"S)
Lower Right (     892.456, -122412.203) (  5d29' 5.80"W,  2d 3'54.60"S)
Center      (     245.375, -121555.162) (  5d29'45.10"W,  2d 3' 2.55"S)
Band 1 Block=1279x1 Type=Float32, ColorInterp=Undefined
    Computed Min/Max=-1452.895,-1373.310
  NoData Value=-3.40282265508890445e+38

Time to process using GDAL colorshade

1.) Hillshade
gdaldem hillshade DTEEC_001414_1780_001612_1780_U01.IMG Victoria_shade.tif -z 2
--note -z is the exaggeration amount

2.) Color ramp -- here we are using  a inverted diverging Blue to Brown from Color Brewer:
this look-up uses percent (as calculated min / max) and then Red Green Blue (8bit values)
create a new file and copy and paste this into it - filename: BlueBrown.lut
nv   0 0 0       //NoData
0%   1   102  94 //brown
20%  90  180 172 //
40%  199 234 229 //light brwon
60%  246 232 195 //light blue
80%  216 179 101 //
100% 140  81  10 //blue

here we are using  an inverted Yellow to Green to Blue from Color Brewer:
create a new file and copy and paste this into it - filename: YlGnBlue.lut
nv   0 0 0       //NoData
0%   37  52 148  //dark blue
20%  44  127 184 //med blue
40%  65  182 196 //light blue
60%  127 205 187 //light blue/green
80%  199 233 180 //light yellow/green
100% 255 255 204 //light yellow
Below you can also use discrete values in the lookup - not just percentages
nv 0 0 0               //NoData
-1452.8945  37  52 148 //dark blue
-1436.9777  44 127 184 //med blue
-1421.0608  65 182 196 //light blue
-1405.1440 127 205 187 //light blue/green
-1389.2272 199 233 180 //light yellow/green
-1373.3104 255 255 204 //light yellow
--note nv is for the NoData value - here set to black. Use 255 255 255 for white.

3.) Apply color ramp to DEM - this is a temporary file which can be deleted.
gdaldem color-relief DTEEC_001414_1780_001612_1780_U01.IMG BlueBrown.lut Victoria_ClrTemp.tif

4.1) merge the colorized file and hillshade (ASP routine)
hsv_merge Victoria_ClrTemp.tif Victoria_shade.tif -o Victoria_colorshade.tif

4.2) merge the colorized file and hillshade (GDAL Python routine in Anaconda - slower)
python Victoria_ClrTemp.tif Victoria_shade.tif Victoria_colorshade.tif

OR process using single ASP routine - color blend might be a little darker

colormap -s Victoria_shade.tif --legend -o Victoria_colorshade.tif --colormap-style BlueBrown.lut DTEEC_001414_1780_001612_1780_U01.IMG
--note there is a legend ramp but you will need to add min and max text values in a photo editor app. See the min and max values from the "gdalinfo" run above.


This is not the best DEM for these colors but feel free to check out color brewer to make your own!

Wednesday, March 1, 2017

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.
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 code. I find that filename hard to read "*" (where trn must mean terrain). I would prefer simply the script name "" (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 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...

Sunday, September 18, 2016

Irregular Bodies and GIS

For using irregular bodies in typical GIS mapping applications -- the current situation is not good. Without playing tricks, I would currently have to recommend using a specialize true 3D GIS-like application (see a couple listed below).

So deriving the body shape is currently fairly well understood (Nefian - presentation, paper, Machtoub - paper). This can be done with lidar and/or stereo, slope, or gravity (some older references).

But for GIS, the issue is that most mapping applications are not well suited for working with irregular shapes (what we call true 3D. We call raster DEMs/TINs only 2.5D). Now there are some GIS-like initiatives from our colleagues at APL (small body mapping tool), JAXA, JPL (vestatrek) and ASU's (JMars).

So what to do if you still wish to use a more typical 2.5D GIS application...

1.) For bodies that are not crazy differently from an ellipse/triaxial shape, the easiest method would be to use a best fit sphere. This is essentially what the IAU recommends for most odd bodies. Now for peanuts (or dumbbell) shaped bodies, this won't work too well (e.g. comet 67P/Churyumov-Gerasimenko or images from Rosetta).

2.) So for "peanuts", there has been research to take these irregular bodies and unwrap it for mapping applications. This generally means warping the irregular body into one Cartesian system (see: Stooke, Morphographic). As a graphic example, I like this video for unwrapping a face for a 3D texture in Blender - although sort-of creepy (see 10.5 minutes into the video).

3.) Others are testing methods to split up the body into multiple map projections. How to best split up the body would still be difficult and may require several different map projections to minimize large distortions for each part of the body. I think this might be explained for comet 67P by F. Preusker et al. but the paper is locked down. This seems the most reasonable solution since one can use well known map projection equations and applications will be able to still run stats (accurate measures) on the different parts. But displaying all the part together might be a challenge since they would essentially be defined as independent bodies. 

A color image of Comet 67P/Churyumov-Gerasimenko composed of three images taken by Rosetta's scientific imaging system OSIRIS. The images were taken on 6 August 2014 from a distance of 120 kilometres from the comet. Credits: ESA/Rosetta/MPS for OSIRIS Team MPS/UPD/LAM/IAA/SSO/INTA/UPM/DASP/IDA

I would love to hear other ideas for this topic.


Wednesday, August 3, 2016

Does the NASA PDS image format make for a good interoperable format?

It could be but that is not its purpose in life.  This post was inspired by question, "why not make PDS exporting (or conversion from Geotiff) available from the NASA AMES Stereo Pipeline." This was my answer. Please note this is IMHO - please feel free to disagree.


So I assume there are many PDS3 (and now some PDS4) writers out there. I actually don't know of a robust GeoTiff to PDS image converter. In general, I view the PDS (3 or 4) format as really an "archival" format and not necessarily a good format to support interoperability between applications. Okay so what's the difference? Archive formats are meant to be "simple" to use and should have verbose labels to access the data for years into the future. Sounds great - what's the problem? In particular, PDS archives can have metadata spread across multiple files which describe the mission, image label particulars, the map projection equations (which may not be  typical equations, if even map projected at all), etc., which makes a single image file less self-contained. Thus not only should a PDS *.IMG (image with label) be created, but to be a valid PDS archive, there should also be other informational files created in a structured set of directories. When GeoTiffs are used for archives (more common for Earth data sets), they actually have the same metadata problem, which is usually solved by using an FGDC/ISO metadata file sitting directly next to the GeoTiff (extension *.tif.xml).

Okay - I don't care about metadata, I just want to write out a PDS *.IMG format. Above, I stated there are probably dozen of PDS writers out there. The main reason there are so many is because the folks writing these files must target their labels for their particular mission. It would be a challenge to create a generic PDS label creation/conversion solution. Many portions of the required metadata cannot be 100% automated (whether PDS or FGDC metadata). For example for PDS, there are many keywords that the user should supply that would need to be (manually) defined during conversion (e.g. PRODUCED_ID, PRODUCER_FULL_NAME, INSTRUMENT_HOST_NAME, etc.). Check out any *.lbl from here.

All that said, both ISIS3 and VICAR can convert to a "generic" PDS v3 label, but it will not be a fully archive-able PDS label. To help with the mission particulars, you will notice ISIS3 has a unique PDS export routine for each instrument (when that instrument is supported for PDS export, e.g. lrocnac2pds, lrowac2pds, mrf2pds - ).

Anyway, I am happy to share a very simple (hacky/brute force) gdal2PDS (v3) Python writer but it would need to be highly tweaked depending on the mission you are supporting. And I currently only support a couple projections as used by the LMMP project. Did I say it was very hacky -- use with caution:

Could a generic PDS3 writer be added to GDAL - yes. But PDS3 is slowly getting supplanted by PDS4 so that would need to be added too. PDS4 is XML based which can complicate simple writers. There are plans to add in a PDS4 reader to GDAL once the format stabilizes. A reader should be much more tangible than a generic writer (for map projected products at least). Heck, we could even support a "raw" Geotiff with a PDS 3 or 4 label pointing inside at the pixels, but raw Tiffs lose many benefits like tiles and compression capabilities.


lots more information about PDS4:

Sunday, July 17, 2016

more HiRISE conversion tips (until labels are fixed)

The HiRISE Team is still planning to fix a few label issues so these work arounds may go away.

background explanation for HiRISE projection issues:,2339.msg9393.html#msg9393

First you need to get some excellent DTMs (aka DEMs) and their orthorectified GeoJpeg2000 images:

HIRISE label issues
(1) MOLA-orthorectified RDRs
For images without partnered DTMs (typical MOLA-rectified images), they have an older map projection definition (where std par is defined as the center lat). There is a simple script called "fix_jp2" to switch the values in-place in the Jpeg2000 header.

download source or Windows binary:
alternative source code:

> gcc -o fix_jp2 fix_jp2.cpp
> fix_jp2 PSP_001918_1735_RED_C_01_ORTHO.JP2
Success, file updated.

(2) partnered DTM-orthorectified images
2.1. first run fix_jp2 as in step 1.

2.2. So a more subtle issue (only in the partnered-DTM ortho-images *.JP2s), is that the radius is not set correctly in the geoJpeg2000 header. This doesn't impact the other map projection parameters (like cellsize or offsets) -- it is really just an incorrect value. To fix this for ArcMap, you can see this link. To fix this for other GDAL environments you need to override the radius in the GeoJpeg2000 file (this generated vrt header also works for ArcMap).

using gdal tools you can create a corrected virtual header for the GeoJpeg2000 file using the projection defined in the DTM.
> gdalsrsinfo DTEEC_001918_1735_001984_1735_U01.IMG -o wkt > fixme.prj

>cat fixme.prj
PROJCS["EQUIRECTANGULAR MARS",GEOGCS["GCS_MARS",DATUM["D_MARS",SPHEROID["MARS_localRadius",3396040,0]],PRIMEM["Reference_Meridian",0],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",283.07],PARAMETER["standard_parallel_1",-5],PARAMETER["false_easting",0],PARAMETER["false_northing",0]]

>gdal_translate -of vrt -a_srs fixme.prj PSP_001918_1735_RED_C_01_ORTHO.JP2 PSP_001918_1735_RED_C_01_ORTHO.vrt
Input file size is 6368, 10767

This process will work for the full-res "A" ortho files too, e.g. PSP_001918_1735_RED_A_01_ORTHO

BTW, if you don't like Jpeg2000 images you can use the HiRISE team JP2_PDS to convert to raw PDS format or the free Kakadu tools to convert to GeoTiff/GeoBigTiff. e.g.

> kdu_expand -i PSP_001918_1735_RED_C_01_ORTHO.JP2 -o PSP_001918_1735_RED_C_01_ORTHO.tif


To map project these all to the same coordinate system (use the corrected vrt) and get the latest GDAL binaries and for each file run the lines below (for GDAL binaries, I like Anaconda Python environment [once installed type: 'conda install gdal'] or Kyng Chaos for Mac, or Osgeo4W - Windows). These steps should work for DTM *.IMG files also.

-only if you don't have GDAL with jpeg2000 support:
> kdu_expand -i PSP_001918_1735_RED_C_01_ORTHO.JP2 -o PSP_001918_1735_RED_C_01_ORTHO.tif
-- again fix radius for converted Tiff
> gdalsrsinfo DTEEC_001918_1735_001984_1735_U01.IMG -o wkt > fixme.prj
> gdal_translate -of vrt -a_srs fixme.prj PSP_001918_1735_RED_C_01_ORTHO.tif PSP_001918_1735_RED_C_01_ORTHO.vrt

in Degrees (equivalent ~1m/p):--download "Mars 2000 Sphere.prj"
> wget
> gdalwarp -t_srs "Mars 2000 Sphere.prj" -tr 1.70670654e-5 1.70670654e-5 -r bilinear PSP_001918_1735_RED_C_01_ORTHO.vrt PSP_001918_1735_RED_C_01_ORTHO_degrees.tif

Or in Meters (e.g. Simple Cylindrical or Equirectangular w/ center at 0, out cellsize set to 1m):
> echo 'PROJCS["Mars_Equidistant_Cylindrical",GEOGCS["Mars 2000",DATUM["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,0]],PRIMEM["Reference_Meridian",0],UNIT["Decimal_Degree",0.0174532925199433]],PROJECTION["Equidistant_Cylindrical"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",0],UNIT["Meter",1]]' > Mars2000_Equi0.prj

> gdalwarp -t_srs Mars2000_Equi0.prj -tr 1 1 -r bilinear PSP_001918_1735_RED_C_01_ORTHO.vrt PSP_001918_1735_RED_C_01_ORTHO_1meter.tif


Tip: To check out your handy work for the  DTM and fixed Ortho use TuiView or 3D Analyst. For 3D Analyst, simply load the images in and under the DTM and Ortho, set the Base Height Tab to use the DTM - almost instant 3D:

Tuesday, July 7, 2015

Profiling MOLA / LOLA Shots (part 2)

Based on the previous post 1, I was pointed to a couple abstracts where the profile was plotted on the feature. See a figure from either abstract below.

abstract 1: or my ancient abstract.

I had created an old ArcView 3.x avenue script to convert a point to side profile line on the map but I didn't want to necessarily rewrite that script in ArcMap, QGIS or OGR. So here is a simple trick which should also work (perhaps with a little parameter tweaking).

1.) Get your MOLA / LOLA shapefile (see part 1).

2.) Like part 1, select a single orbit and now export to a new shapefile (right click on layer, data export...).

3.) Open the table for new shapefile and calculate statistics for "PtTopo", write that number down. For my orbit, it is -1380.45

4.) Start editing the exported track (so you can undo it if it doesn't work).

5.) In the table, right click on "Shape" field and select "Field Calculator...", Change the language to Python at the top and select "show code block".

6.) In the code block copy this code:
def shiftXCoordinateWithZ(shape, ptTopo, xScale, meanTopo):
    point = shape.getPart(0)
    centeredZ = ptTopo - meanTopo
    point.X += centeredZ * xScale
    return point

In the next area, which should be below "shape =" paste this and edit the values for Xscale and meanTopo. XScale worked for LOLA but depending on your topographic range you may beed to adjust this value.
shiftXCoordinateWithZ(!SHAPE!, !PtTopo!, 0.0002, -1380.5 )

7.) If the updated points makes for a good profile, save edits. If not, then stop editing and don't save edits. Start editing again and try to adjust the "xScale" (the second number in the function call). This value applies an exaggeration across the X (or longitude). A higher number will exaggerate more. The "meanTopo" values tries to center the points along the original line.

8.) If you like the profile size but want to show these points as a line, use the tool Points to Line to convert the points to a profile like. You can also smooth the line for a cleaner looking profile.


Profiling MOLA / LOLA Shots (part 1)

When MOLA was first released as PDERs (shot data), creating profiles in a GIS over surficial features was a little bit of a challenge. You had to download 20GB of the original shot "PEDR" data and then run the pedr2tab tool to generate an exported table. Edit the table to make it a recognized table format, load into ArcView 3.x, assign a Mars radius, and finally convert to a GIS shapefile. Many scripts were written to convert these tables to shapefiles for programs like ArcView 3.x, ENVI, and eventually GDAL's ogr2ogr. It is now as easy as visiting a PDS website to get the shot data. The GeoScience Node has a service to skip all these steps and go directly to a shapefile (or a shapefileZ - where each point has a geometry with X, Y, and Z).

With more than 600 million shots for Mars you will need to only request a region. I think they have a soft limit of about 6 million shots but it is better to obviously keep it smaller than that. For LOLA, the Team will probably reach a 8 billion shots by the end of the mission!

Huh? Why not just profile over the available MOLA and LOLA digital elevation models (DEMS)? Of course that is a very suitable and recommended method and can be done in many applications like ArcMap, JMars, ENVI, QGIS, Saga GIS, Mirone, GRASS, Matlab. Soon we are releasing a new geodesic DEM profiling tool for ArcMap which can also help with geologic cross-section generation. But this means you are profiling over filled or interpolated data where shots were not originally available. Sometimes profiling over the shots still makes sense.

Anyway, once you get the shot data from the PDS as GIS shapefiles or simply a table, profiling the data can be as easy as selecting one orbit in a GIS or even Excel and plotting the "PtLat" Latitude field versus the "PtTopo" elevation height field (for LOLA). Note that mixing latitude (in degrees) against elevation in meters sounds odd but it is a good proxy since the orbit of the instruments is/was polar. Ideally both axes being in meters would be best.

In ArcMap or other GIS, simply load your shapefile shots (via "add Data"). Query for a single orbit you want to profile. Now use the select by attribute to select that orbit and in the case of LOLA you should choose a single Shot Number (left). You really only want to profile over one orbit per graph. Now create a graph using the type "vertical line" or scatter. Set the Y axis as "PtTopo" and X axis as "PtLat" (right). In the graph wizard, hit "Next" and select the option "Show only selected features on the graph." (below). 

You can now finish the plot or continue to customize the chart's title or options. Note my example plot is garishly colorized since I symbolized the original shots in ArcMap. This is handy but adding the legend isn't very helpful on the chart so I turned it off. You can steal the legend from the table of contents for figures if you need to. To export this chart, you can right click on it and export to various graphic formats or PDF. 

Tip: Select a different orbit and the graph should update too. Select features in the graph and they should highlight on the map.

Tips: one of the benefits of a shapefileZ is the ability to also profile the whole lot in ArcScene or other capable 3D viewer. You can plot the non-Z shapefiles too by just selecting the "PtTopo" field as the height source.

Continue to part 2.