What is the ISIS 3 cube format

ISIS 3 Cube Format - http://isis.astrogeology.usgs.gov/ or Wikipedia

Overview: ISIS (Integrated Software for Imaging Spectrometers) is a generalized software system that has been designed to optimize cartographic and scientific processing of images in planetary datasets. ISIS versions 2 and version 3 both have their own format called an ISIS cube which contains a label area and a data area. The label describes the data and uses a "keyword=value" text format. The label is normally embedded at the top of the file, but ISIS also supports "detached" labels where the label and image block are divided into separate files.

The data area is a 3-dimensional image with axis: samples, lines, and bands. Typically, the sample and line dimensions are used to represent spatial information while the band dimension represents spectral information. Pixels can be stored in 8-bit unsigned integer, 16-bit signed integer, or 32-bit floating-point format. When cubes are created in ISIS3, unless specifically specified, the file will be written using tiles (blocks of pixels). Tiles can be variable in size as defined in the label. In addition to valid numerical pixel values, five types of "special" pixel values are supported: null, high/low instrument saturation, and high/low representation saturation.

Type	Description	8bit Unsigned	16bit Signed	32bit Floating-Point
Null	Pixel has no data available	0	-32768	\xFF7FFFFB (-3.40282265508890445e+38)
Lis	Pixel lower bound was saturated on the instrument	0	-32767	\xFF7FFFFC
His	Pixel higher bound was saturated on the instrument	0	-32766	\xFF7FFFFD
Lrs	Pixel lower bound was saturated during processing	255	-32765	\xFF7FFFFE
Hrs	Pixel higher bound was saturated during processing	255	-32764	\xFF7FFFFF

Applications with read support for ISIS 3 cubes

1. IsisDlm (ISIS IDL Interface Documentation) - http://www.astro.cornell.edu/~carcich/LRO/doc/html/main.html 
2. Geospatial Data Abstraction Library (GDAL). http://www.gdal.org/, http://www.gdal.org/frmt_isis3.html 
• Applications which use GDAL and thus have read support for ISIS cubes 
• GDAL conversion utilities - http://www.gdal.org/gdal_utilities.html 
• QGIS - http://www.qgis.org/ 
• TuiView (simple viewer) - https://bitbucket.org/chchrsc/tuiview 
• ArcMap/ArcGIS Pro - http://www.esri.com/software/arcgis 
• SAGA GIS - http://www.saga-gis.org/ 
• UDig GIS - http://udig.refractions.net/ 
• GRASS GIS v7 - http://grasswiki.osgeo.org/wiki/Planetary_mapping 
• Generic Mapping Tool (GMT) - http://gmt.soest.hawaii.edu/ 
• Mirone - http://w3.ualg.pt/~jluis/mirone 
• more - http://trac.osgeo.org/gdal/wiki/SoftwareUsingGdal 
3.  Davinci (AZ State University) - http://davinci.asu.edu/index.php?title=ISIS3_I/O 
• JMars - http://jmars.asu.edu/ 

Simple command-line GDAL Conversion examples 

• Installers: http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries
• Binary tools in a Python environment (my favorite method to get binaries). Install the free Anaconda Python and then once install type "conda install gdal".

Examples:

$  gdal_translate –of GTiff input_32bit.cub output_32bit.tif

$  gdal_translate –of GTiff input_16bit.cub output_16bit.tif

$  gdal_translate –of GTiff input_8bit.cub output_8bit.tif

$  gdal_translate –of Jpeg –ot Byte –scale input.cub output_8bit.jpg

$  gdal_translate –of PNG –ot Byte –scale input.cub output_8bit.png

more from ISIS Wiki: https://isis.astrogeology.usgs.gov/fixit/projects/isis/wiki/Isis_Cube_Format
 related: http://planetarygis.blogspot.com/2014/06/using-isis-23-image-cubes-in-gis.html

more examples: https://isis.astrogeology.usgs.gov/IsisSupport/index.php/topic,2172.0.html

Figure 1. Illustration of an ISIS image cube, showing the three dimensions of an Isis image: width (samples), height (lines), and depth (bands).

Testing a Catalog Service (CSW) for a Planetary Portal

The topic for this post is to discuss a potential method for multiple data portals to share resources (and data) via OGC Catalog Services. Now I am using portal here loosely as many of the sites that could be a proper data portal are currently just a web site (even if build from a dynamic database).

Goals:

1. Create “proper” geospatial metadata for science-ready (map projected) derived products including file-base (PDS, GeoTiffs, etc.) and live maps (WMS, WFS, etc).
2. Test a freely available, open source OGC CSW service to host this metadata.
3. Using this service, enable the ability for one site to catalog another site, called “harvesting”.
4. When other facilities harvest these records, make sure proper credit is given to the originator and allow downloads from the original sites. Here we are trying to minimize data redundancy (and issues like who’s is newer or different).
5. While websites can facilitate data discovery across facilities, still teach users how to query using these CSW via available command-line methods and from GUIs (like QGIS).

---

• Create proper metadata

To keep this discussion short I am going to skip over tons of metadata topics, but to be able to utilize existing CSW infrastructure, I am going to recommend using a FGDC metadata standard. Why – because we want to catalog geospatial data (and 99% of the standard works for planetary). But we could also more simply use the Dublin Core or NASA’s DIF which both seem well supported in the existing catalog software.

As part of a proof-of-concept, I have very rough program called gdal2metadata.py (code in github) to help generate (hopefully) valid FGDC metadata. I still recommend using the USGS tool “mp” to actually validate output. Only ~30% of metadata can automatically gathered from any one map projected file (in our case a GeoTiff or PDS image). Thus this code requires a “template” metadata record to be mostly written by the team responsible for the data set. If you keep the abstract, purpose, and other parts pretty generic this can be done for many related files. Once that is written, the script can run using this template on many files. It needs a lot of updating, for example to support more projections, but it is a start. Also in the github repository there are a couple template examples for Lunar DEMs. Anyway, once a Python environment is installed (w/ GDAL) you can simply run: (BTW, I recommend the awesome Anaconda environment and install GDAL using “conda install gdal”).

> ./gdal2metadata.py NAC_DTM_ATLAS5.TIF ASU_DTM_Template.xml

Warning: currently I don’t have the actual download location in the metadata but I will add it in soon. Also any data particulars for specific images, like accuracy or a listing of images used would have to be inserted via this script or another script to update those areas.

---

• OGC CSW service

Using the generated metadata record above, I was able to quickly get running with the free pycsw software. And I also cheated and used the OSGeo Live (v8) virtual machine for an immediate working environment. Fortunately, along-with the virtual machine, the pycsw team has published a walk-though workshop (which I followed specifically from this point).  Here are the few steps I had to do to load in this example FGDC record (It was fairly immediate once I figured out pycsw only wanted a lower case .xml extensions! )

> cd /var/www/pycsw
> export PYTHONPATH=`pwd`
> sudo vi default.cfg
       and change this line to include fgdc:
       profiles=fgdc,apiso
while editing, update all the custom name, place, email, etc. for the facility

Backup sqlite example database and create a new one. 

Note the location and name of the database is located in default.cfg
> sudo rm  /var/www/html/pycsw/tests/suites/cite/data/records.db
normally I would back this up - not simple delete 

> sudo python ./bin/pycsw-admin.py -c setup_db -f default.cfg
> sudo ./bin/pycsw-admin.py -c load_records -f default.cfg -p ~/myTestMetadata/LROC -r
where "-r" is recursive and will load and files it finds in the directory tree.

That is basically it…

Now start Firefox and point to http://localhost/pycsw/tests/index.html, change the pulldown request to GetCapabilities and server to “../csw.py?config=/var/www/html/pycsw/default.cfg” and hit Send. Here are other examples but also try to send a “GetRecord by bbox” using the same server and

<gml:lowerCorner>-90 -180</gml:lowerCorner>
<gml:upperCorner>90 180</gml:upperCorner>

Definitely a start. But now try QGIS with the CSW plug-in also on the OSGeo Live VM (using http://localhost/pycsw/csw.py for the server). Here is a snapshot after running a spatial search for records:

Now need to test the ability for pycsw to automatically harvest layers from a WMS getcapabilities. Also on this page is how to access these records using Python on the command-line.

---

• Summary (thus far)

There are still plenty items to test including making sure the metadata is pushed through properly and WMSs are cataloged properly. But this seems to an method using mostly existing tools to share data across facilities. There are also tools which use pycsw to nicely display the catalog like GeoNode, Open Data Catalog, CKAN (used by data.gov)

LPSC abstract on the topic.

to be continued…

HiRISE or LROC NAC DEM to printable 3D model

I have not actually had a 3D model produced on a 3D printer but this method should create a file ready to print. The issue I encountered was not the ability to create a 3D scene and/or converted to a proper format, but a 3D model as a solid required by 3D printers. When making a model using a digital elevation model (DEM) also called digital terrain model (DTM) it is easy to make an undulating surface but "flat". 3D printers need to be a solid which is a 3D surface with a "skirt" under it (see image).

For this tutorial I used this LROC NAC DTM called unnamed fresh crater east of Lents: http://wms.lroc.asu.edu/lroc/view_rdr_product/NAC_DTM_NEARLENTS1

So the steps are:
1. Create a 3D model in ArcScene, Global Mapper, Terragen, Blender, etc. Note in the image the surface is just rendered – meaning it is not yet a solid model. Tip: square up 3D scene to clip out any NoData area. In ArcScene you use the AOI button to clip the scene to a drag-able 3D box.


2. Export Scene to a 3D format – e.g.; VRML (*.wrl), Standard Tessellation Language (*.stl), or Polygon File Format(*.ply, also know as the Stanford Triangle Format).

So now for the trick - make the undulating (but flat surface) into a solid using the free software ZPrint.

3. Open ZPrint and import the 3D model. This is where you have to know your printer's capable dimensions in width and height. For this model I set millimeters and scale by 2% which ends up being about 300mm (~1 foot in Y).

4. In ZPrint, select Edit – "Make Solid". Here the defaults generally work. I now have a solid (see image) which can be printed on a compatible ZPrinter or exported to a STL or PLY for printing elsewhere.


If you need another 3D format for your type of 3D printer, import the PLY or STL 3D solid into the excellent and free Meshlab and export to another format like Alias Wavefront (*.obj). PLY seems to work best for Meshlab. Lastly, it seems like it might be possible to make a solid in Meshlab but I haven't been able to get it work as easy as ZPrint.

Have fun 3D printing.

---

Update - I later created a model of the full crater and it was printed in 3D at peak-solutions-llc.com. It is nice to see the image printed on the model - not all 3D printers can do that yet.

Update 2 - Now the online Moontrek viewer (JPL) can output STL or OBJ for you!

http://sservi.nasa.gov/articles/moon-trek-beta-released/ 

ArcMap volume methods for small features

Over the years I have calculated volumes for various features. This tutorial will deal with tips for trying to calculate the volumes of dunes. This is more a manual method than a automatic batch method for many dunes and is more applicable to features found in high-resolution data sets like HiRISE DEMs and LROC DEMs.

There are several methods in ArcMap to do this from just using image statistics to tins, terrains, profiles, etc.  For example, in batch, you can do "zonal stats to table" to get the mean elevation (or other parameter) for each digitized dune polygon and then times the returned mean by area of dune (using a table join). However, this would be very much an approximate for the dune's volume.

Again I will explain manual methods for ArcMap users who have the 3D Analyst extension. (1) calculate volume for a dune on a flat surrounding surface (2) calculate volume from a sloped surface.

Tip: For HiRISE, if the orthophoto doesn't register to DEM see: https://isis.astrogeology.usgs.gov/IsisSupport/index.php/topic,3729.msg14926.html#msg14926 or you can try fix_jp2.
Tip: Example for testing download the DEM and full-res Ortho: http://www.uahirise.org/dtm/dtm.php?ID=ESP_021893_1315

1) Calculate volume for a feature on a flat surface.

1.1) Create an empty polygon file for dunes (ideally in Feature Class or maybe just a Shapefile)
http://resources.arcgis.com/en/help/main/10.1/index.html#/in_a_feature_dataset/002200000004000000/

1.2) Once all dunes are bounded by polygons, select single polygon, export to graphic:
http://resources.arcgis.com/en/help/main/10.1/index.html#//00s90000001v000000
--Turn off polygon layer and select graphic using standard black arrow tool
--With single dune graphic selected, Right click on DEM - Export to clip DEM
(method 2 easy): https://www.esri.com/arcgis-blog/products/product/analytics/clipping-an-image-or-raster-in-arcgis/

1.3) With clipped DEM, now run "Surface Volume" from Toolbox (requires 3D Analyst):
http://resources.arcgis.com/en/help/main/10.1/index.html#//00q900000027000000
--Done, volume reported in output table in m³


2) calculate volume from tilted surface

Summary: The idea here is to actually remove (clip out) the dune from DEM, interpolate across dune hole in DEM (using one of several methods available) , then use "cut fill" tool to subtract dune from tilted interpolated-over dune hole.

optional tip: Using a graphic box with bounds dune with a little buffer, clip DEM to smaller DEM with Dune and some extra. Or at any processing stage you can set the processing extent to the "Current Display" if zoomed in or another layer with the bounding extent.

Do the same steps 1.1 and 1.2 above.

2.3) From 1.2, again select same single dune graphic, Right click on DEM – Export to clip DEM (but use "outside" to create dune hole).
more clipping help: http://blogs.esri.com/esri/arcgis/2011/10/05/clipping-imagery-and-raster-data-in-arcgis/

2.4) convert created DEM (dune hole) to multi-points:
http://resources.arcgis.com/en/help/main/10.1/index.html#/Raster_To_Multipoint/00q900000067000000/

2.5) run Create Tin with points. This will fill the hole with Tin triangles - spanning the gap – making a generally very flat surface where the dune use to be.
http://resources.arcgis.com/en/help/main/10.1/index.html#/Create_TIN/00q90000001v000000/

2.6) Convert Tin back to raster which should return a new DEM with again the flatten and tilted surface where to dune was.
http://resources.arcgis.com/en/help/main/10.1/index.html#/TIN_To_Raster/00q900000077000000/

2.2.7) Run "Cut Fill" on "1.2" (polygon clipped dune DEM) and "2.6" (hole-filled DEM).
http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000vt000000
--Open resulting table for new raster layer and select the largest volume in the table. The volume will again be m³.

No example pictures added just yet but I will try.

-Trent

Using ISIS 2/3 image cubes in GIS software

Below contains recommended methods to directly use or convert ISIS 3 (and ISIS v2) into a well-formed GIS compatible format.
Recommended Rules of the Road:

• Before ingesting any file (PDS, ISIS, other) into a GIS it should be map projected.
• If file is to meant for visual purposes only (e.g. mapping on), it is highly recommended to stretch the images to 8bit (see below)
• Don't use funky projections (e.g.; ISIS's oblique cylindrical)
• Best to use positive East longitudes and a –180 to 180 longitude system.
• For Longitude systems 180 or 360 it is highly recommended to:
• if the set longitude system = 180 (-180 to 180), then set center longitude = 0 (or center of the image).
• if the set longitude system = 360 (0 to 360), then set center longitude = 180 (or center of the image).
• Don't use ISIS's method, isis2std, to convert files to Tiffs, Jpeg2000, Jpgs or PNGs. Why? It does not create a fully-formed GIS format with defined projection. It supports a GIS worldfile (simply the image registration in the current projected Cartesian plane) but not the projection definition.
• Do use GDAL's tools, gdal_translate, to convert files from ISIS3 to one of many output formats (GeoTiff, GeoJpeg2000, Jpeg, PNG, ENVI, ect.).

Direct support: There is fairly good support for reading ISIS3 cubes (and version 2) files in the GDAL library. This means packages which use GDAL can usually read most map projected ISIS 2,3 and PDS files directly. This includes packages like QGIS, ArcMap GIS, Mirone, Saga GIS, Opticks, GRASS, and soon GMT, more). Since the GDAL reader only support one NoDATA value for 32bit and 16bit files there might be issues for the application to understand the valid data range. In this case, it is best to use stretch in ISIS to set the 5 supported special pixel values to valid or NoDATA (see below).
Conversion support: As stated above, I recommend using the GDAL utility application gdal_translate to convert from ISIS2, 3 or PDS format to another format for your application. Forth coming thread on GDAL tips for planetary data…
a few conversion examples:

• gdal_translate –of GTiff input_32bit.cub output_32bit.tiff
• gdal_translate –of GTiff input_16bit.cub output_16bit.tiff
• gdal_translate –of GTiff input_8bit.cub output_8bit.tiff
• gdal_translate –of Jpeg –ot Byte –scale input.cub output_8bit.jpg
• gdal_translate –of Jp2kak –co quality=100 input_16bit.cub output_16bit.jp2

Conversion from 32bit to 8bit in ISIS before conversion:
My currently favorite method to convert to 8bit in ISIS (before using gdal to convert to geoTiff or GeoJpeg2000 or other format). I can't recommend the ISIS program bit2bit since it pushes valid data which is clipped into NoDATA (not good for viewing the image).

ISIS Stretch method 1:>stretch from=input.cub to=output_8bit.cub+8bit+1:254 USEPERCENTAGES=true pairs="0:1 100:254" null=0 lis=1 lrs=1 his=255 hrs=255
    This allows you to specify input percentages for the mapping pairs. Thus when
    USEPERCENTAGES=true is set pairs="0:1 100:254" means:
    map 0% to 1 (or the file's min value to 1) and 100% to 254 (file's max value).

ISIS Stretch method 2:This also means you can apply a recommended 0.5% clip to remove the potential extraneous lows and highs like:
>stretch from=input.cub to=output_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 once in an 8bit format convert to GeoTiff or other for easy and fast viewing. For huge files take advantage of BigTiff or Jpeg2000 lossy support

• gdal_translate –ot GTiff –co Bigtiff=if_safer –co tiled=yes input_8bit.cub output_8bit.tif
• gdal_translate –ot Jp2kak input_8bit.cub output_8bit.jp2

----simple example Csh Script to force output scale to 1 to 255.
mag{90}> more /usgs/cdev/contrib/bin/to8bit_gdal_tif.csh

Code: 

#!/bin/csh
#usage: to8bit_gdal_tif.csh input.ext output.tif
#Note the use of $3 in gdal_translate which just allows you to add in other options. By default it will just be blank.
set in=$1
set out=$2

set stats = `gdalinfo -stats $in > /tmp/xxxtemp_stats.junk`
set min = `cat /tmp/xxxtemp_stats.junk | grep MINIMUM | awk -F= '{print $2}' | sed 's/ //g'`
set max = `cat /tmp/xxxtemp_stats.junk | grep MAXIMUM | awk -F= '{print $2}' | sed 's/ //g'`

gdal_translate -ot byte -of GTIFF -co compress=lzw -ot Byte $3 -a_nodata 0 -scale $min $max 1 255 $in $out
/bin/rm -f /tmp/xxxtemp_stats.junk 
-Trent

Planetary GIS – why start a blog now?

For more than 10 years, apparently since 2002, I have used various discussion sites to list common GIS issues for the community. To get the ball rolling, many posts were just questions sent to me in emails that I would try to answer. I actually started posting the questions from a generic user called “planetary researcher” and then logged in as myself, “thare”, and answered the question. Using the wayback machine, it appears in the second discussion site incarnation, the top read threads included answers to:

• How do I read the MOLA MEGDR 128 [pixels/degree] images?
• How to get a level 2 ISIS image into a GIS?
• Polar glitches in the MDIM2.1 files?

Obviously the Mars MOLA topography was still the new kid on the block and very popular (frankly still is). And while some GISs now have direct support for ISIS files, loading map projected ISIS and PDS files is still a common issue for GIS users. And the “polar” glitch thread was more about the new MDIM2.1 projection style changing to Equirectangular from Sinusoidal and the fact the Jpeg images had too many pixels for Photoshop CS ( version 8 ) . BTW, Looking at those “small” file sizes now, newer images released continue to explode in pixel count and size. We now have several 1TB mosaics in the community! The challenge is not really how to open them, but how to deliver them.

Anyway, that second discussion site incarnation lasted only year before the ISIS started a discussion forum and all the GIS topics were co-located with the ISIS topics. But in that one year more than 100 thousand threads were viewed. I assume many hits were from general GIS users who quickly bounced once they figured out we were discussing planets (or false hits from indexers).

Anyway, to get to the point, our latest discussion software platform has been down for about a week from an apparent PHP breach. While many topics even on the latest site become stale, there are still plenty of useful threads for folks to stumble upon. I will start transferring some of the those older topics (but still relevant) and post new topics into this blog. We still get plenty of email questions weekly which posts here can be based on. I don’t see using a blogging tool really any better than a discussion site for keeping threads from becoming stale, but perhaps I can remove or edit sections more easily.

For my first real post I will try to summarize how to load/convert ISIS files – again.

-Trent