PC Align

Last month we had a new release of Ames Stereo Pipeline, version 2.3! We’ve performed a lot of bug fixing and implementing new features. But my two new prized features in ASP are pc_align and lronac2mosaic. Today I’d like to only introduce pc_align, a utility for registering DEMs, LIDAR points, and ASP point clouds to each other. All you have to do is specify an input, a reference, and an approximate estimate for how bad you think the misplacement is.

Does that sound like magic? Under the hood, pc_align is performing an implementation of the iterative closest point algorithm (ICP). Specifically we are using internally libpointmatcher library from ETH. What ICP does is iteratively attempt to match every point of the input with the nearest neighbor in the reference point cloud. ICP then solves for a transform that would globally reduce the distances between the current set of matches. Then it repeats itself and performs a new round of matching input to reference and then again solves for another global step. Repeat, repeat, repeat until we no longer see any improvements in the sum of distances between matches.

This means that failure cases for PC_align are when the reference set is too coarse to describe the features seen in the input. An example would be having a HiRISE DEM and then only having a single orbit of MOLA that intersects. A single line of MOLA does nothing to constrain the DEM about the axis of the shot line. The 300 meter post spacing might also not be detailed enough to constrain a DEM that is looking at small features such as dunes on a mostly flat plane. What is required is a large feature that both the DEM and the LIDAR source can resolve.

CTX to HRSC example

Enough about how it works. Examples! I’m going to be uncreative and just process CTX and Gale Crater because they’re fast to process and easy to find. I’ve processed the CTX images P21_009149_1752_XI_04S222W, P21_009294_1752_XI_04S222W, P22_009650_1772_XI_02S222W, and P22_009716_1773_XI_02S223W using stereo options “–alignment affineepipolar –subpixel-mode 1”. This means correlation happens in 15 minutes and then triangulation takes an hour because ISIS subroutines are not thread safe. I’ve plotted these two DEMs on top of a DLR HRSC product, H1927_0000_DT4.IMG. It is important to note that this version of the DLR DEM is referenced against the Mars ellipsoid and not the Aeroid. If your data is referenced against the Aeroid, you’ll need to use dem_geoid to temporarily remove it for processing with pc_align who only understand ellipsoidal datums.

In the above picture, the two ASP created DEMs stick out like sore thumbs and are misplaced by some 200 meters. This is due to pointing information for MRO and subsequently CTX being imperfect (how much can you ask for anyway?). You could run jigsaw and that is the gold standard solution, but that takes a lot of manual effort. So instead let’s use pc_align with the commands shown next.

> pc_align --max-displacement 200 H1927_0000_DT4.tif P22-PC.tif \
        --save-transformed-source-points  -o P22_align/P22_align \
> point2dem --t_srs "+proj=sinu +lon_0=138 +x_0=0 +y_0=0 \
        +a=3396000 +b=3396000 +units=m +no_defs" \
        --nodata -32767 P22_align/P22_align-trans_source.tif

There are 3 important observations to make from the command line above. (1) We are using the max displacement option to set the upper bound of how bad we think we are, 200 meters. (2) I’m feeding the ASP PC file as the input source instead of ASP’s DEM. This is because in (3) with the save transformed source points we’ll be writing out another PC file. PC align can only export PC files so we always have to perform another round of point2dem. It is possible to run stereo -> point2dem -> PC Align -> point2dem, but it means you are unneccesarily resampling your data once. Using the PC file directly from stereo saves us from potential aliasing and removes a point2dem call.

Here’s the final result where both CTX DEMs are plotted on top of the HRSC DEM. Everything looks really good except for that left edge. This might be because the DEM was rendered with incorrect geometry and the output DEM is subtly warped from a perfect solution.

Another cool feature that pc_align author Oleg Alexandrov added was recording the beginning and ending matching errors in CSV files. They’re found with the names <prefix>-{beg,end}_errors.csv. You can load those up in QGIS and plot theirs errors to visualize that pc_align uniformly reduced matching error across the map. (Thanks Ross for showing me how to do this!)


Quite a few MOLA shots can be found inside the CTX footprints. Above is a plot of all MOLA PEDR data for Gale crater. I was able to download this information in CSV format from the MOLA PEDR Query tool from Washington University St. Louis. Conveniently PC_align can read CSV files, just not in the format provided by this tool. PC_align is expecting the data to be in format long, lat, elevation against ellipsoid. What is provide is lat, long, elevation against aeroid, and then radius. So I had to manually edit the CSV in Excel to be in the correct order and create my elevation values by subtracting 3396190 (this number was wrong in first draft) from the radius column. The other added bit of information needed with CSV files is that you’ll need to define the datum to use. If you don’t, pc_align will assume you’re using WGS84.

> pc_align --max-displacement 200 P22-PC.tif mola.csv\
        -o P22_mola/P22_mola --datum D_MARS --save-inv-trans \
> point2dem --t_srs "+proj=sinu +lon_0=138 +x_0=0 +y_0=0 \
        +a=3396000 +b=3396000 +units=m +no_defs" \
        --nodata -32767 P22_mola/P22_mola-trans_reference.tif

Two things to notice in these commands, the inputs are backwards from before and I’m saving the inverse transform.  You can keep things in the same order as when I was aligning to HRSC,  it is just that things will run very slowly. For performance reasons, the denser source should be considered the reference and then you must request the reference to be transformed to the source. You’ll likely routinely be using this inverse form with LIDAR sources.

In the end I was able to reduce alignment error for my CTX DEMs from being over 50 meters to being less than 15 meters against MOLA and from over 100 meter to 40 meters error against HRSC. A result I’m quite happy with for a single night processing at home. You can see my final composited MOLA registered CTX DEMs on the left. The ASP team will have more information about pc_align in LPSC abstract form next year. We also hope that you try out pc_align and find it worth regular use in your research.


I goofed in the MOLA example! Using D_MARS implies a datum that is a sphere with 3396190 meter radius. I subtracted the wrong number from MOLA’s radius measurement before (the value 3396000). That probably had some effect on the registration result shown in the pictures, but this mistake is smaller than the shot spacing of MOLA. Meaning the horizontal registration is fine, but my output DTMs are 190 meters higher than they should have been. FYI, D_MOON implies a datum that is a sphere with radius 1737400 meters.

Mass CTX Processing

A few weeks back, Ross Beyer presented my blog posts on autonomous HiRISE DEM processing at the HiRISE team meeting in Tucson, AZ. This brought about a question of could this be performed for CTX. Of course! ASP can be applied to do bulk processing for any of the missions it supports. Earth or any place ISIS and CSpice have defined coordinate system for. Just put in some safeties about run time into the processing scripts because ASP still occasionally goes mad and eats a whole bunch of processing time for no output. (Were working on it!)

Processing CTX stereo pairs however is in fact a little more difficult compared to the HiRISE processing I was doing before. HiRISE lists all their stereo pairs from their website or Dr. Shane Byrne’s website. There’s no equivalent for CTX. Luckily for me, some folks at UofA wrote PairendipityCTX (Chris Schaller?). They provided Ross and I a detailed report of overlapping files and other statistics. I cut out everything but the filenames for my own use and you can get a copy of the list with its 1,542 stereo pairs here.

Another difference was how these two missions stored their data on PDS. I can look at a HiRISE filename and work out its download path in PDS with no trouble. CTX on the other hand seems to have arbitrary volume breaks, thus the download URL is not predictable. My solution is a bad solution, but a quick solution. I wrote a python script that scraped PDS’s servers and identified all the CTX images it has and what their URLs are. I then just ‘grep’ for the URL in the processing scripts. Here’s the resulting text file that lists all of the 50,708 CTX images that existed at the time of my scraping. This is a mean trick because my script can make HTTP requests much faster than a human can. In a sense, I’m DOS’ing the PDS servers. So please copy my results rather than my methods.

Processing scripts

Previously with my autonomous HiRISE processing efforts, I just wrote a Bash script and then spawned it multiple times on my server using GNU parallel. That’s great, but looking around my office I saw a lot of unused computers that I’d like to do my bidding. What I wanted was a job management system like PBS, which is what they use on the super computer. But PDS is a little too heavy and sometimes cost money. So I instead found two alternative projects that seemed to fit the bill, Gearman and Celery. Gearman is the minimalist. Celery required a database backing and multiple ports open with each slave worker. I decided to use Gearman since it seemed simpler to learn.

My code for this little project is available on Github in my Mars3DGearman project. I hope to eventually add HiRISE support. Here’s how it plays out. All machines make a working directory that contains the CTX stereo list, the CTX url lookup list, and the folders DEM and DRG. One machine is designated the server, in my case, my workstation at home. It starts the only instance of ctx_processor.py with the help of a backing ‘gearmand’ executable (gearman daemon). All of the slaves then SSH back to my server and forward 4730, the port used by Gearman for communication. Then all the slave machines can start one or more instances of ctx_worker.py.

The part I haven’t worked out is how to send home the completed results to my main server. The ctx_worker script produces a DEM and orthophoto and then just dumps it locally into the DEM and DRG folder. Gearman does allow sending binary strings back the main server, but I’m betting a 100 MB long string would break something. I chose a different route. Since I already have all the slaves SSH’ing back to my main server, I decided to simply rsync all the slaves’ DEM and DRG folder back home. I don’t have to re-enter my password as I’ve enabled SSH ControlMaster which re-uses previous connections. For now, I just put that rsync in a watch command that fires every 2 hours. Alternatively it could be inside the ctx_worker script. A better bet would be to use SSH keys.

Another worthwhile improvement compared to my HiRISE processing scripts is the inclusion of a timeout method for each step. When it comes to CTX, if the correlation doesn’t finish in 2 hours, I’m not interested. This timeout is achieved through the run_cmd and process_timeout functions in the ctx_worker script. The Internet helped me a lot in figuring out how to make that a reality. Thanks Internet!


These last few days I’ve roped 4 machines into doing my bidding. So far they’ve produced ~260 DEMs. 5 new DEMs completed just while I was writing this article. There are still some hiccups in the process. But when the stars align, I seem to produce over 50 new DEMs every day. I’m not going to show you all of them as that felt like a lot of work to post on to the blog. Instead I’m just going to show off a couple screenshots of some interesting places in Google Mars. The color ramp is a little funky because someday I need to learn to reference everything against the Mars Aeroid and not just the sphere.

Not everything looks as great as those screenshots. Here are some places that failed to correlate. I’m not sure what went wrong and I unfortunately don’t have the time to investigate.

In conclusion, this is just another proof of concept of what is possible. I hope that someday someone else will attempt this and do a better job than me. ASP is not perfect, but it can achieve a lot of processing on its own that could be beneficial to the scientific community.


Shean, D. E., et al. “MRO CTX Stereo Image Processing and Preliminary DEM Quality Assessment.“ Lunar and Planetary Institute Science Conference Abstracts. Vol. 42. 2011.

David Shean, MSSS, and Larry Edwards actually already attempted this once before! In the above abstract you can see that they went above and beyond my little demo and processed 1180 stereo pairs. They also developed some initial steps for registering the output DEMs to MOLA and plotted the relationship convergence angle has on the outcome of a stereo pair.

Creating a registered LRO-NAC DEM without GCPs

I’ve previously written a tutorial on performing a bundle adjustment, jigsaw in ISIS lingo, using CTX. Let’s try something a little bit more complex with LRO-NAC. If you are not previously aware, LRO-NAC is actually two separate optical assemblies that are joined by the spacecraft bus. Ideally this spacecraft hardware would be solid as a rock and wouldn’t bend. Unfortunately that’s not the case, so what would have been a bundle adjustment between 2 cameras is now 4 in the case of LRO-NAC.

LRO-NAC Data Preparation

Picking a stereo pair has always been the most difficult part of the job for me. I honestly don’t know how to find stereo pairs. So instead I’m just going to recreate something ASU did in SOCET SET. If you are not aware, ASU has processed about 100 DEMs with LRO-NAC. I’m going to recreate their Moore F Crater DEM using M125713813 and M125720601. You can download the raw data via PDS.

I picked this stereo pair blindly and it ended up being extremely difficult to register correctly. In order to help out, I added a fifth image to the bundle adjustment to give another orbit’s worth of spacecraft ephemeris into the measurements. That image was M110383422LE and I picked it because it overlapped 2 of my previous images and had the same lighting conditions as my previous images.

Once you have downloaded all the images, prep the files in ISIS with the following:

parallel "lronac2isis from={} to={.}.cub; lronaccal from={.}.cub to={.}.cal.cub; spiceinit from={.}.cal.cub; rm {.}.cub" ::: *IMG

* The ‘parallel’ command is GNU Parallel, a non-standard system utility that I love dearly.

Attaching a custom Elevation Source

From a previous article I wrote about, we saw that map projecting against the WAC Global DTM gave much better results than working against the sparse LOLA that is the default for map projecting by in ISIS. So I’m going to do that now. Once we have a successful jigsaw result, then we’ll map-project. If you forgot how to attach the WAC Global DTM, here are the commands I used. My image pair was at 37.3N and 185E. So the tile WAC_GLD100_E300N2250_100M.IMG from this link was what I needed.

pds2isis from=WAC_GLD100_E300N2250_100M.IMG to=WAC_GLD100_E300N2250_100M.cub
algebra from=WAC_GLD100_E300N2250_100M.cub to=WAC_GLD100_E300N2250_100M.lvl1.cub operator=unary A=1 C=1737400
demprep from=WAC_GLD100_E300N2250_100M.lvl1.cub to=WAC_GLD100_E300N2250_100M.lvl2.cub
rm WAC_GLD100_E300N2250_100M.cub WAC_GLD100_E300N2250_100M.lvl1.cub

parallel spiceinit from={} shape=USER model=$PATH_TO/WAC_GLD100_E300N2250_100M.lvl2.cub ::: *cal.cub

Making the control network

Making the control network is pretty straightforward and is just time consuming. I performed the commands I previously outlined here on this blog.


Group = PolygonSeederAlgorithm
      Name = Grid
      MinimumThickness = 0.01
      MinimumArea = 1
      XSpacing = 600
      YSpacing = 1000


Object = AutoRegistration
   Group = Algorithm
     Name         = MaximumCorrelation
     Tolerance    = 0.7

   Group = PatternChip
     Samples = 19
     Lines   = 19
     MinimumZScore = 1.5
     ValidPercent = 80

   Group = SearchChip
     Samples = 120
     Lines   = 100


parallel footprintinit from={} ::: *cal.cub
echo *.cub | xargs -n1 echo > cube.lis
findimageoverlaps from=cube.lis overlaplist=overlap.lis
autoseed fromlist=cube.lis overlaplist=overlap.lis onet=control.net deffile=autoseed.def networkid=lronac pointid=???? description=moon
pointreg fromlist=cube.lis cnet=control.net onet=control_pointreg.net deffile=autoRegTemplate.def
qnet (editting and fixing control_pointreg.net)

I had to spend a lot of time in Qnet for this stereo pair. About half of my automatic matches didn’t work. (Possibly my search range was too small?) In order to correct this, I had to then manually edit the output control network in Qnet. I then filter down all the control points to the ones marked ‘Ignored’. Then I just manually align the measurements. You just need to get the measures close to the right spot by clicking on the same feature. Afterwards, you can then click the ‘Register’ button to get subpixel accurate alignment.

Autoseed and AutoReg also have a bad habit of only matching LE images to LE, and RE images to RE. It is a good idea to run along the edge of an LE image and try to find common points that are seen in both RE images. These are going to be the few control points that have more than 2 control measures. I also had to do this between my 5th image and the other 2 images that it had overlap with. This problem was likely because the overlap between these 3 images didn’t meet the requirement of MinimumThickness in Autoseed.def.

Normally at this point I would then recommend to you to start picking out ground control points. I tried that. Then I went a little insane. I’m better now. Instead we’re going to do crazy stuff in the next section.

The reason we’re not using GCPs, is that the only source to align against is the WAC Image mosaic. I tried this and I could easily get a bundle adjustment result that aligned well with the image mosaic. However the output DEM would always show a massive error between my height results and the results of LOLA and the WAC Global DTM. After adding the 5th image and clamping the trust of my GCP to 50 meters in radius, I got a result that matched LOLA and WAC DTM in an okay manner but didn’t really align at all to the WAC image mosaic. Things made sense once I compared the image mosaic to a hillshade of WAC DTM.

WAC’s Image mosaic doesn’t agree with WAC’s Global DTM! Ghaaaa! Madness! I was so incredibly disappointed. I didn’t know what to do at this point, so I just dropped the ground control points.

Bundle Adjusting

Alrighty then, time to start bundle adjusting with jigsaw. Likely on the first run things are going to blow up because bad measurements still exist in the control network. My solution to this was to first have jigsaw solve for only a few camera parameters and then check out the output control network for features that had the highest error. In Qnet you can filter control measures that have errors higher than a certain amount of pixels. After I fixed those measurements manually, I resave the control network as my input network and then redo the jigsaw.

Each time I redo jigsaw, I let jigsaw solve for more parameters. My pattern was:

jigsaw fromlist=cube.lis radius=no twist=no cnet= control_pointreg_gcp.net onet=control_jigsaw.net update=no
jigsaw fromlist=cube.lis radius=no twist=yes cnet= control_pointreg.net onet=control_jigsaw.net observations=no update=no
jigsaw fromlist=cube.lis radius=yes twist=yes cnet= control_pointreg.net onet=control_jigsaw.net observations=no update=no
jigsaw fromlist=cube.lis radius=yes twist=yes cnet= control_pointreg.net onet=control_jigsaw.net observations=no update=no camsolve=velocities

After iterating many times, I eventually produced a control network and jigsaw result that had a sigma0 of .34 pixels. However, if we’re to write this camera solution to file and then make a DEM, we would see that I have worse alignment with LOLA than if we hadn’t bundle adjusted at all.

So how can we get a good bundle adjustment result that aligns well with LOLA? One idea would be to pick GCPs directly against LOLA using their gridded data product. I know folks in the USGS who have done exactly this. However I don’t even see how they can determine correct alignment, so I’m going to ‘nix that idea. Another idea would be to have Jigsaw not solve for radius (radius=no), but instead use only the radius values provided by LOLA. This will keep problem in the ball park, but the resulting jigsaw solution will have a sigma0 around 15 pixels! Not solving for radius essentially ties our hands behind our back in the case of LRO-NAC when our radius source is so incredibly low res.

What we need instead is to simply constrain / rubber-band our Control Points to have radius values that are similar to LOLA. Unfortunately, how do we determine what correct Lat Long to sample the LOLA radius for each control measure? My solution was to not define that. I did two things: (1) I changed all my control measures from ‘free’ to ‘constrained’ with a sigma of 200,200,50. This can be performed with cneteditor. (2) I then wrote and ran the following bash script.

#/bin/env bash


cp $input_control control_loop.net

for i in {0..100}; do
    echo Iteration $i
    cnetnewradii cnet= control_loop.net onet= output.net model= $radius_source getlatlon= adjusted
    mv output.net control_loop.net
    jigsaw fromlist=cube.lis radius=yes twist=yes cnet= control_loop.net  onet= output.net update=yes spsolve= position camsolve= velocities
    mv output.net control_loop.net


Your mind should be blown at this point. I’ve essentially written an iterative closest point algorithm using the ISIS utilities jigsaw and cnetnewradii. This process is only able to work because we have a perfect control network that we vetted with the previous jigsaw runs. This also works because our LOLA data here is not flat, and has plenty of texture. Finally, LRO-NAC was close to the correct solution. Imagine images with horrible attitude data, like Apollo, would just blow up with this technique.

On each iteration, the sigma0 shouldn’t change. Or if it changes, it should be extremely small. That’s because a translation or scale change of the control points and camera location has no effect on the reprojection error. What should be reducing on each iteration is the standard deviation of the radius correction for each control point as listed in bundleout.txt. The sparse 3D model of the problem in Jigsaw is slowly fitting itself to the LOLA DEM but constrained only by the camera geometry. The standard deviation of the radius correction will never go to zero and will eventually plateau. This is because there will always be some error against LOLA since large portions of the measurements are interpolations between orbits.

Disclaimer: I’ve only tried this technique for this LRO-NAC stereo pair. I don’t know if it will work for others. However this was the only technique I could come up with that would produce results that best aligned with LOLA. I couldn’t get anything good to come out of using WAC Image Mosaic GCPs.

Map projecting for Stereo

Now that we have an updated ephemeris, we’re ready to perform a map projection to speed up stereo. Be sure to not accidentally spiceinit your files, which would wipe your updated ephemeris. Notice that I’m map projecting everything at a low resolution. This is because I’m in a hurry for this demo. In a production environment for LRO-NAC, we would likely map project at 1 MPP or 0.5 MPP. This step and ASP’s triangulation are usually the longest part of this entire process.

parallel cam2map pixres=mpp resolution=5 from={} to={/.}.map.cub ::: ../*cal.cub

Running ASP

Nothing special here, I’m just processing the LE-LE, RE-RE, and LE-RE combination. RE-LE stereo pair wouldn’t resolve anything at 5 m/px. Afterwards, I gdalbuildvrt the outputs together. The massive lines separating the individual DEMs that make up this stereo pair may concern you. However that effect will diminish if I had map projected my inputs at a higher resolution. The reason being that, correlation is much like convolution, in that we erode the input image by a half kernel size along the image border. That border gets perceivably smaller if the image is of higher resolution.

parallel -n 2 stereo M125713813{1}.cal.map.cub M125720601{2}.cal.map.cub {1}{2}/{1}{2} ::: LE LE RE RE LE RE
parallel -n2 point2dem -r moon --t_srs \"+proj=eqc +lat_ts=37 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +a=1737400 +b=1737400 +units=m +no_defs\" --tr 5 --orthoimage {1}{2}/{1}{2}-L.tif {1}{2}/{1}{2}-PC.tif --nodata -32767 ::: LE LE RE RE LE RE
gdalbuildvrt aspba_drg.vrt */*-DRG.tif
gdalbuildvrt aspba_dem.vrt */*-DEM.tif
parallel gdal_translate -co TILED=yes -co COMPRESS=LZW -of GTiff {} {.}.tif ::: *.vrt

Conclusion and Comparison

Now the moment you’ve been waiting for. The point where I have to prove I’m not crazy. I have three versions of this stereo pair. There is the result produced by ASU using SOCET SET and then there’s ASP with and without Bundle Adjustment. The look of ASP’s result will be a little shoddy because I map projected my input at 5 m/px which prevent me from get the most detail out. It did however save me a lot processing time while I experimented this method out.

Here are the hillshades for those 3 datasets overlaid a hillshade of the LRO-WAC Global DTM which is 100 m/px. There’s not much to tell, but you can at least see that everyone seems to have decent horizontal alignment.

Here are the difference maps between those 3 datasets and the LRO-WAC Global DTM.

Here are the difference maps between those 3 datasets and the LOLA Gridded Data Product, which I resampled to 100 m/px. I think this shows straight forward the improvement bundle adjustment made to the output DEM.

Just for giggles, I also did a difference between the ASU DEM and the bundle adjusted ASP DEM. They’re pretty close to each other. There seems to just be a simple tilt between the two or that the ASP DEM is just shifted to a slightly lower latitude. Who’s more correct is anyone’s guess for now.


Ideally at this point I would then compare the ASU DEM and the BA’d ASP DEM against the raw LOLA shot points. This is where curiosity got me and where I move on to the next personal project. Maybe someone else would be interested in a more rigorous review.

Also, I’m not sure if the 5th image is required.

MER Data Introduction

Getting data from the NASA Planetary Data Services (PDS) can be a little intimidating. The most import thing a person needs to understand about the MER data is that it is sectioned into individual sites that the rovers visited during their trip. These sites are localized on interesting features like an individual bay into a large crater or a particularly large bolder. These sites are further subdivided into sols (Martian days) at the site. When we search for images we limit them by both sites and sols.

Gathering Data

It is now time for us to select a site for which we want to render in 3D. I find MER’s Analyst’s Notebook as a good tool for this task. Here’s what I did:

  1. Go to: http://an.rsl.wustl.edu/mer/
  2. Click on “Opportunity”.
  3. To find a site: click on the “Map” icon on the toolbar, and then click on “Traverse Map”.

Each one of the black dots represents a site location. The number after the slash is a sol number. I compared this map to Google Earth (GE) and decided that I wanted to render “Cape Agulhas”. You can do that too in GE if you click on the flag icon for a MER and select “load rover way points”. Just note that the map on GE doesn’t show as much information as the traverse map on Analyst’s Notebook.

  1. On the “Traverse Map”, click “91/1673”.
  2. Select only Navcam data products on the left.
  3. On the Data Products drop down box, select “Show for sol 1674”
  4. Click “Redraw Map”.

From here we can see that the Navcam did a full panorama. We are now ready to start downloading data from PDS. You can download from Analyst’s Notebook, however I was only able to download a single image at a time. On NASA PDS we can download a ‘wget’ script that will download all Navcam images for us in one go.

  1. Go to NASA PDS’s image search for MER, http://pds-imaging.jpl.nasa.gov/search/search.html#QuickSearch
  2. On the left, under “Select Instrument(s):”, select “NAVCAM.
  3. Next to “Instrument Host ID”, select “MERB/MER1/Opportunity”.
  4. Next to “Product Type”, select “EDR”.
  5. Next to “Image Type”, select “Regular”.
  6. Next to “Eye”, select “LEFT”.
  7. Next to “Planet Day Number”, type 1674 for both the Min and Max text boxes.
  8. On the left, click “Get Results”.

You now have a listing of all the Left EDR NAVCAM. We want all of these; so on the left select ‘WGET’ under download products. Then click download. Move your downloaded ‘atlas_wget_script’ to your work directory. Here is how you start the download in the terminal:

source atlas_wget_script

Later on I can show how to produce 3D maps from this imagery. However this seems like a good start.

Fly by night Cassini ISS processing

Anytime I want to process Cassini ISS data, I always run into problems. This is largely because the documentation is lacking. Today is the day that is going to change. Here are my notes on how to extract ISS data of Iapetus, all for the low cost of free. I assume you have USGS’s ISIS3 installed.

Iapetus, file N1483152391_1

Downloading from NASA PDS

As usual, all NASA data can be pulled free from their PDS (Planetary Data System). Since I’m looking for Iapetus, I click ‘Saturn’ and then ‘Cassini Image Search’. This will bring you to the Cassini Image Quick Search page. On the left side of the page, select the ISS instrument. Then select EDR under product type. EDR products are raw images while RDR are processed products like maps. Finally under target name, select the option Iapetus. Click ‘Get Results’ on the far right. At the time of writing, there was 4025 records, or images available.

Here’s where things get a little messy. A lot of these images are of Iapetus from very far away. So far away that they wouldn’t be good for my projects. I’d like to sort these images by their distance. On the far right, under SELECT PARAMETERS FOR REPORT OR TABLE COLUMNS; I can add the column ‘Target Distance’. On the new column, I can click the up arrow and hopefully it will sort the images by who is closest. However it doesn’t work. In the middle of the page; PDS says “508 of 4025″ products. The page is only sorting 1/8th of the total images of Iapetus. To correct this, keep clicking the ‘Get More’ button until it says “4025 of 4025″ products.

At this point we are ready to start downloading products. On the far right of the page, download the CSV file of the records they found along with the WGET product. When you open up the CSV file you’ll notice that they didn’t keep the images ordered in distance from target. Lame, I know. However, we use Linux and we can pull ourselves out of this. Here’s the bash commands I used to sort the file. Gurus can probably do a better job of this than me.

cat atlas_report.csv | awk -F "," '{print $3 $4}' | sort -n -k 2 > sorted_files

If you open up sorted_files in your favorite text editor, you can clean up the stray lines that were mis-sorted. I chose to keep the files with distances marked ‘-1.0e3′. I also hand deleted the lines where the distances were marked greater than 300 km out from the moon. This left me with 758 names of images that are pretty close to Iapetus. Now we need to pull these 758 names from the atlas_wget_script so that I only download the images I want. However when you look in that script, you’ll notice that they name the images slightly different from what we saw in the atlas_report.csv. Here’s an example:

1_N1466586085.118 versus N1466586085_1.IMG

Here’s how I filtered sorted_files down to only that mantissa like bit I wanted and then pulled out the respected lines I wanted out of atlas_wget_script.

cat sorted_files | awk '{print substr($1,3,11)}' > wanted_keywords
cat atlas_wget_script | grep -f wanted_keywords > wanted_get_script

Hooray! Finally we are ready to really download imagery files! If you look inside wanted_get_script, you’ll see the download commands for all the images of Iapetus that I want. To start downloading:

source wanted_get_script

Processing with ISIS 3

You probably think it’s all kisses and hugs from here. You have a direct full of IMG and LBL_label files do ya? Well the command we want to use is ciss2isis. Well then, don’t let me hold you back, read their instructions and give it shot on one of the images you downloaded.

ciss2isis from=N1483150900_1.LBL_label to=test.cub

What it broke?

**I/O ERROR** Unable to read PVL file [/Volumes/Data/Saturn/CASSINI_data/Iapetus/temp/N1483150900_1.LBL_label]
**PVL ERROR** Error in pvl on line [166]
**PVL ERROR** Unable to read keyword [DESCRIPTION    = "For the packet from which this telemetry header was Extracted:  the first line of this packet is a continuation of a line begun in the previous packet. OBJECT = COLUMN NAME = NULL_PADDING DATA_TYPE = MSB_UNSIGNED_INTEGER START_BYTE = 61 BYTES = 475 END_OBJECT = COLUMN END_OBJECT = TELEMETRY_TABLE OBJECT = LINE_PREFIX_TABLE INTERCHANGE_FORMAT = BINARY ROWS = 256 COLUMNS = 7 ROW_BYTES = 24 ROW_SUFFIX_BYTES = 512 OBJECT             = COLUMN NAME               = LINE_NUMBER DATA_TYPE          = LSB_UNSIGNED_INTEGER START_BYTE         = 1 BYTES              = 2 DESCRIPTION        = "The image line number of this record.  Maintained at]
**PVL ERROR** Keyword has extraneous data [The image line number of this record.  Maintained at] at the end

This looks darn scary. However if you download the original LBL files from the PDS manually without using the atlas_wget_scripts ‘s version, the problem becomes clear. The download script was accessing a new label file that had a poorly inserted ‘EXTENDED_ISS_SCIENCE_HEADER’. This was probably done by some batch processing job at PDS or Cassini imaging lab. This mistake on PDS managed to delete some of the lines defining the end of objects. This caused the ISIS3 command ciss2isis to freak out an die.

Yet don’t worry .. this can be corrected. I have written a patch for these LBL_label files that will make then usable. In a new file, write the following contents.

--- N1483150900_1.LBL_label     2010-10-19 19:46:20.000000000 -0700
+++ N1483150900_1.LBL_corr      2010-10-19 19:45:01.000000000 -0700
@@ -99,6 +99,7 @@
START_BYTE     = 1
BYTES          = 60

NAME           = CAMERA
@@ -144,7 +145,8 @@
BITS           = 1
DESCRIPTION    = "For the packet from which this telemetry header was
Extracted:  the first line of this packet is a
-                    continuation of a line begun in the previous packet.
+                    continuation of a line begun in the previous packet."
@@ -209,6 +211,8 @@
START_BYTE         = 9
BYTES              = 2
      LINES = 256

I saved the above as correction.patch. It can then be applied to all LBL_label files to make then work.

echo *label | xargs -n1 echo | xargs -n1 -I{} patch {} correction.patch

You will now find that ciss2isis works. Here’s how I do that with multiprocessors.

echo *label | xargs -n1 echo | awk -F "." '{print $1}' | xargs -n1 -I{} -P8 ciss2isis from={}.LBL_label to={}.cub

You now have a directory with cub files. You can view those with the qview command or export them with isis2std. There’s still more that needs to be done inorder clean up this imagery. I won’t tell you today, but I can hint what the next few steps will be. Have fun!

Iapetus plus stars