ISIS3 on Windows 10 WSL Ubuntu (Bash)

-------->>>> Please see NEW 2024 version here. 

WARNING: text below is deprecated -------

Based on a tweet by MarsLakes, I also decided to give ISIS3 on Windows 10 a try. Note this is for Windows 10 Anniversary update or higher. Anyway, here is my current cheat sheet. It works better that I thought it would. Warning: this requires a lot of downloading.

Updated July 2019, to support ISIS3 (v3.6+) in an Anaconda environment.

(1) get Linux for Windows 10

I recommend Ubuntu 18 LTS installed via the Windows 10 app store, but other WSL flavors may work also.

And if you don't have access to the App store, try to download:

https://docs.microsoft.com/en-us/windows/wsl/install-manual

(2) For Windows, install Xming and launch (it runs in the background)

http://www.straightrunning.com/XmingNotes/

(3) update Ubuntu (this just gets your Ubuntu system up to date)

> sudo apt-get update
> sudo apt-get upgrade

(4) From Windows, start bash and let's get the GUI working ( can be one line also, listing them all).

> sudo apt-get install -y -q xclip gnome-themes-standard gtk2-engines-murrine
> sudo apt-get install -y -q dbus dbus-x11 x11-apps libgl1-mesa-glx
-- test X11
-- add this line to bottom of .bashrc (or type):
> echo export DISPLAY=:0 >> ~/.bashrc
> source ~/.bashrc
> xclock

(5) Get Anaconda

> wget https://repo.anaconda.com/archive/Anaconda3-2019.03-Linux-x86_64.sh
-- or get latest version from download page: https://www.anaconda.com/distribution/
> bash Anaconda3-2019.03-Linux-x86_64.sh
-- follow install instructions and allow installer to update your shell (~/.bashrc) when it asks.
-- once installed, restart the Ubuntu terminal from Windows start.

(6) Get ISIS3 applications

-- here you can follow the ISIS3 installation pages: https://github.com/USGS-Astrogeology/ISIS3
or these steps:
> conda create -n isis3 python=3.6
> conda activate isis3
> conda config --env --add channels conda-forge
> conda config --env --add channels usgs-astrogeology
> conda install -c usgs-astrogeology isis3
> conda activate isis3
-- execute the ISIS3 variable initialization script with default arguments.
-- this script prepares default values for: $ISISROOT, $ISIS3DATA, $ISIS3TESTDATA
> python $CONDA_PREFIX/scripts/isis3VarInit.py

(7) Get ISIS3 base data
> conda activate isis3
> cd $ISIS3DATA
-- note the "." at the end of the next line. This is the destination and it is required (or a dir. path).

> rsync -azv --delete --partial isisdist.astrogeology.usgs.gov::isis3data/data/base .
-- note the "." at the end. This is the destination location and it is required (or a dir. path).
-- try opening something
> qview $ISIS3DATA/base/dems/MSGR_DEM_USG_EQ_I_V02_prep.cub
-- done for the minimal install. You may need lots more SPICE data for a particular mission, read more about that towards the end for the install page: https://github.com/USGS-Astrogeology/ISIS3
-- So the next time you launch an Ubuntu terminal, just type "conda activate isis3" to run ISIS3 apps.

(8) to run a test CTX image in ISIS3

-- first get some additional MRO files

> cd $ISIS3DATA

-- note the space and period at the end for the next line:
> rsync -azv --delete --partial --exclude='kernels' isisdist.astrogeology.usgs.gov::isis3data/data/mro .
-- now get a scripted example
> cd ~

> sudo apt-get install csh

> mkdir isis3_test

> cd isis3_test

> wget https://raw.githubusercontent.com/USGS-Astrogeology/ISIS3_scripts/master/MRO_CTX_process_all/CTX_process_all.csh

> chmod +x CTX_process_all.csh

> wget -nd http://pdsimage.wr.usgs.gov/Missions/Mars_Reconnaissance_Orbiter/CTX/mrox_0674/data/P22_009816_1745_XI_05S073W.IMG

> printf "Group = Mapping\n  ProjectionName  = SimpleCylindrical\n  LongitudeDomain = 180\n  CenterLongitude = 0.0\n  CenterLatitude  = 0.0\nEnd_Group\nEnd" > simp0.map

> ./CTX_process_all.csh simp0.map  0
-- wait several minutes for ISIS3 to do its thing.

> qview P22_009816_1745_XI_05S073W.lev1eo.cub

figure 1. qview with CTX P22_009816_1745_XI_05S073W calibrated only

> sudo apt-get install gpicview

> gpicview P22_009816_1745_XI_05S073W.png

figure 2. gpicview with CTX P22_009816_1745_XI_05S073W.png all map projected

(9) Now use GDAL to create a GIS-ready image (GeoTiff) not just a PNG. 

-- note ISIS3 and gdal cannot be in the same anaconda environment just yet. Let's make another independent environment for GDAL
> conda create -n gdal
> conda activate gdal
> conda install -c conda-forge gdal=3

-- first stretch to 8bit using ISIS3 (0.5% to 99.5% linear stretch -- optional):

> conda activate isis3

> stretch from=P22_009816_1745_XI_05S073W.lev2.cub to=P22_009816_1745_XI_05S073W_8bit.cub+8bit+1:254 USEPERCENTAGES=true pairs="0:1 0.5:1 99.5:254 100:254" null=0 lis=1 lrs=1 his=255 hrs=255

-- now get back to GDAL's environment. Jumping between conda environments can be scripted as shown here to bounce back and forth.

> conda activate gdal

> gdal_translate P22_009816_1745_XI_05S073W_8bit.cub P22_009816_1745_XI_05S073W_8bit.tif

> gpicview P22_009816_1745_XI_05S073W_8bit.tif

Updating to the next version of ISIS3 using the latest Anaconda methods will generally only take minutes now. You can follow ISIS3 Release Notes on Astrodiscuss to get alerted to updated versions.

But when you are ready, simply run:
> conda update isis3

And if you even need to downgrade (or specify a particular version), simply list the version number like:
> conda install isis3=3.7.1

tip: Now run ISIS commands from (really a shell script) from Powershell. I guess this means ArcMap or Socet GXP might be able to run ISIS3 commands.... :-)

From Windows, start powershell (not the (x86) version), from the horrible blue background (do yourself a favor and turn the background black), run:

> bash -i -c "conda activate isis3 ; /home/thare/isis3_test ; ./CTX_process_all.csh simp0.map  0"

-- so the trick here is the "-i" which setups the env. for ISIS3 and "-c" runs the command in quotes.

-- Also the Windows locations for /home/thare is here (you may need to show hidden files in Windows to see it): C:\Users\<WIN_USER>\AppData\Local\lxss\home\thare\

-- In Bash, your Windows files are: /mnt/c or /mnt/d, etc.
-- more tricks: https://www.howtogeek.com/262086/how-to-run-linux-commands-from-outside-the-bash-shell-on-windows-10/

reference: https://msdn.microsoft.com/en-us/commandline/wsl/faq

tip: more Python 3.x packages for Anaconda with GDAL

> conda activate gdal

> conda install -c conda-forge jupyter basemap matplotlib xerces-c kealib

-- note putting all the packages on the same line, can help with getting all dependencies correct.

When you run from your Ubuntu side:
> jupyter notebook
you can then run your browser on the Windows-side using the link provided in Ubuntu - pretty cool.

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. https://store.continuum.io/cshop/anaconda/ .

• For a more plain GDAL environment Windows binaries see: https://trac.osgeo.org/osgeo4w/
• For Mac: http://www.kyngchaos.com/software:frameworks
• For Linux: It's not too hard to download and build yourself.

If you just want ASP there are binaries for Mac and Linux: https://ti.arc.nasa.gov/tech/asr/intelligent-robotics/ngt/stereo/

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: http://www.uahirise.org/dtm/dtm.php?ID=PSP_001612_1780

wget http://www.uahirise.org/PDS/DTM/PSP/ORB_001400_001499/PSP_001414_1780_PSP_001612_1780/DTEEC_001414_1780_001612_1780_U01.IMG

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:
PROJCS["EQUIRECTANGULAR MARS",
    GEOGCS["GCS_MARS",
        DATUM["D_MARS",
            SPHEROID["MARS_localRadius",3396190,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",354.5],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-401.706464281687772,-120698.120618077693507)
Pixel Size = (1.011855073757400,-1.011855073757400)
Metadata:
  BANDWIDTH=
  CENTER_FILTER_WAVELENGTH=
  DATA_SET_ID="MRO-M-HIRISE-5-DTM-V1.0"
  FILTER_NAME=
  INSTRUMENT_ID="HIRISE"
...
  TARGET_NAME="MARS"
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:
http://colorbrewer2.org/#type=diverging&scheme=BrBG&n=6
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:
http://colorbrewer2.org/#type=sequential&scheme=YlGnBu&n=6
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)

wget https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/hsv_merge.py
python hsv_merge.py 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.

Done.

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

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...