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

for MOLA shots: http://ode.rsl.wustl.edu/mars/indexTools.aspx?Displaypage=molapedr

for LOLA shots: http://ode.rsl.wustl.edu/lunar/indexTools.aspx?Displaypage=lolapedr

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.

-Trent

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