Phocube missing essential output options (local phase angle, slope, slope azimuth)

[lavalaz]

Description

The name of the phocube application suggests it is for doing photometric analyses which is the comparison of phase angle versus reflectivity. The most important photometric parameter is the phase angle, corrected for topography. This would be called "LOCALPHASE" in the parlance of phocube. Phocube has LOCAL_EMISSION and LOCAL_INCIDENCE which means all the math to calculate local phase angle is being done, but there is no option to output the one parameter that really matters. It is essential to also output the local incidence and local emission angle values to allow proper error analysis of the local phase angles.

The request is to add three (3) new output options in phocube:

• LOCALPHASE (phase angle computed for each image pixel, using the DEM. Same math as calculating the phase angle without the DEM, but using the local incidence and local emission angle)

• SLOPE (slope from the DEM that is used to calculate LOCALEMISSION, LOCALINCIDENCE, and LOCALPHASE; this is already being computed for each pixel, just asking for an option to output it as a band)

• SLOPEAZIMUTH (azimuth of the downslope direction from the DEM that is used to calculate LOCALEMISSION, LOCALINCIDENCE, and LOCALPHASE; this is already being computed for each pixel, just asking for an option to output it as a band)

Example

phocube from=EW0131773041G_cal.cub to EW0131773041G_cal.pho.cub dn=yes localphase=yes slope=yes slopeazimuth=yes

[AaronBoyd]

This is incorrect. The phase angle calculation is independent of topography, and a "localphase" does not need to be calculated in isis, as it is the same calculation as the normal phase.

Phase angle is the angular seperation between the camera look vector and the solar vector, and does not change with local slopes.

NO CHANGE NEEDED FOR PHASE ANGLE.

While the idea of being able to calculate local slope with phocube is enticing, it is likely better for the user to know what the baseline is of the slope being calculated. If the user were to create a slopemap with known dimensions, the user could run map2cam to get the slopes in image space.

I would recommend not adding slope to phocube.

[lavalaz]

Well, duh, Aaron is completely correct. I was being a total doofus and thank you for commenting to catch that. So, LOCALPHASE is a stupid non-existent concept, please ignore.

However, I would really still like to see SLOPE and SLOPEAZIMUTH because I am having trouble sorting out the source of mismatches between where I see craters in the image (i.e., the DNs) and where the LOCALEMISSION and LOCALINCIDENCE seem to be responding to slopes. I can run slpmap, but I cannot confirm that everything is identical to how phocube does its slope calculations. Exposing the numbers phocube is using for topography would be very helpful in figuring out where the issue is.

Thank you!

[jessemapel]

@lavalaz I am very hesitant to expose information so that users can try and debug the software. Would you be okay with a software developer checking that phocube and slpmap are using equivalent methods?

For a little bit more information, here's how phocube computes the local photometric angles:

1. Compute the intersection with the shapemodel (usually a DEM) at the middle of the pixel
2. Compute the normal vector for the shapemodel at the intersection
3. Computes the vector from the intersection to the center of the sun at the time the time it was viewed
4. Computes the vector from the intersection to the sensor at the time the time it was viewed
5. Compute the angles between the three vectors

Here's the code for reference

The two main reason your local photometric angles may not line up with the topography you see in the imagery are:

1. Your shapemodel is lower resolution that your imagery
2. Your imagery is misaligned with your shapemodel

[lavalaz]

I need the slope information to compute correlation statistics between slope and DN values for my photometric analysis. While I can generate a slope map outside of phocube and mosaic it into the phocube output, this breaks an existing rigorous chain between the slopes and the DN values at each pixel. And is a lot of extra work to get numbers that are already being computed within phocube. This is not a matter of a user trying to do a programmer's job, it is a case of a scientist trying to do a scientist's job.

[jlaura]

This issue is still open and is a candidate for discussion at the product support sprint prioritization meeting. It looks to me like a solution was presented 2/28 and then a response came in 5/3.

Can we continue the discussion in here prior to the meeting?

[lavalaz]

The "local-phase-angle" request was stupid, please ignore that part.
I have done a work-around for the slope and azimuth, but it is a major pain and is really slow, because it is a lot of data to move around (for larger images). Since phocube has to access all that data already, it would be nice to avoid having to deal with it twice and output it from phocube. The other benefit is that there is no chance of user error pointing to a different DEM for the phocube and slope map.
Just for background, the issue is that the topographic and image data are almost never identical in scale. The resolution of the emission and incidence angles are tied to the DEM, not the image data. So you get a whole host of artifacts when you try to do photometric analysis at the pixel level in the image. In order to understand the issues, the analyst needs to be able to see if strange photometry corresponds to issues with the mismatch between image and topo data or if there is something strange about the surface. In principle, any representation of topography would help - shaded relief, radius, etc. But slope and slope azimuth most directly ties to emission and incidence angle so I suggest that it makes the most sense to use those angles.
I have been discussing visible images, but slope and slope azimuth is also very useful in understanding thermal data.
At the same time, there are many cases when one does not care about the slope and/or slope azimuth so this definitely should be an optional ouput.
Hope this added discussion is helpful!

[jessemapel]

There was some concern during the prioritization about adding more information to phocube but looking over what's already there, I think these would fit.

@lavalaz Can you provide a definition, including the math, for the slope azimuth? The sensor model has the following:

• The position of the sensor (and thus the look vector)
• The position of the sun
• The ground point
• The surface normal vector at the ground point (both from the DEM and the ellipsoid)

[lavalaz]

Short version: it is the compass direction in degrees of the surface normal vector (N,S,E,W in numerical form - most people have north = 0, east = 90, etc., but I believe ISIS has a different convention that I have to look up every time). Let me know if you need a more precise definition...

…

------------- Laszlo Kestay Research Geologist USGS Astrogeology Science Center

________________________________ From: Jesse Mapel <[email protected]> Sent: Friday, September 4, 2020 11:01 AM To: USGS-Astrogeology/ISIS3 <[email protected]> Cc: Kestay, Laszlo P <[email protected]>; Mention <[email protected]> Subject: [EXTERNAL] Re: [USGS-Astrogeology/ISIS3] Phocube missing essential output options (local phase angle, slope, slope azimuth) (#3645) This email has been received from outside of DOI - Use caution before clicking on links, opening attachments, or responding. There was some concern during the prioritization about adding more information to phocube but looking over what's already there, I think these would fit. @lavalaz<https://github.com/lavalaz> Can you provide a definition, including the math, for the slope azimuth? The sensor model has the following: * The position of the sensor (and thus the look vector) * The position of the sun * The ground point * The surface normal vector at the ground point (both from the DEM and the ellipsoid) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#3645 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AKMB27SSAMKMWHJI4BZNMJDSEETPTANCNFSM4KGE53GQ>.

[jessemapel]

Let me know if you need a more precise definition...

A more precise definition would be good. What two vectors is it the compass angle between; the surface normal and what?

[lavalaz]

Its not an angle between two vectors. Its a projection of the vector onto the plane of the map. And then the angle between that vector and a reference direction (north, if you use the normal convention, which isis does not). It is exactly the same as the sub-spacecraft and sub-solar azimuth... I can't give you the math because I don't know the way in which isis describes vectors, planes, map directions, etc. But I'll send some figures in a few minutes... ------------- Laszlo Kestay Research Geologist USGS Astrogeology Science Center

…

________________________________ From: Jesse Mapel <[email protected]> Sent: Friday, September 4, 2020 11:12 AM To: USGS-Astrogeology/ISIS3 <[email protected]> Cc: Kestay, Laszlo P <[email protected]>; Mention <[email protected]> Subject: [EXTERNAL] Re: [USGS-Astrogeology/ISIS3] Phocube missing essential output options (local phase angle, slope, slope azimuth) (#3645) This email has been received from outside of DOI - Use caution before clicking on links, opening attachments, or responding. Let me know if you need a more precise definition... A more precise definition would be good. What two vectors is it the compass angle between; the surface normal and what? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#3645 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AKMB27T4PJZ4CJMWYFNW7LLSEEUX7ANCNFSM4KGE53GQ>.

[lavalaz]

Perhaps the confusion is that slope "azimuth" and slope "aspect" are used interchangeably in the literature. Looks like isis follows the ArcGIS terminology and uses "aspect". Please look at the slpmap application, output = aspect, for the math and description. https://isis.astrogeology.usgs.gov/Application/presentation/Tabbed/slpmap/slpmap.html

[rbeyer]

The term "aspect" looked funny to me (because I usually think of this as azimuth, maybe my astronomy roots), but apparently this aliasing between aspect and azimuth is pervasive. There is a gdaldem aspect command, whose documentation uses the term "azimuth" to describe what the aspect command outputs, so there you go.

[jessemapel]

@lavalaz That's what I was looking for, thank you. We should be able to implement this based on that.

[jessemapel]

I've got an idea of what code needs to be added for this all worked out. @lavalaz @rbeyer Any thoughts on whether to call this slope azimuth or slope aspect? Aspect seems to be more common, but phocube in general uses azimuth for calculations like this (the similar options are subspacecraftgroundazimuth & subsolargroundazimuth).

[rbeyer]

I'd make sure to have both in the documentation, and be clear about it (unlike the gdal docs). Clearly, you have to pick one word for a key or an argument, but in the documentation be clear that "aspect" and "azimuth" are synonyms here. This way people familiar with just one term or the other can both "get it."

As far as which one to use, there are arguments both ways. If phocube already uses "azimuth" in this exact same way for other things, then maybe staying with that pattern makes sense?

[jessemapel]

I've got this all coded and I'm looking at the results and seeing some interesting things. This picture compares the slope computations from phocube with the slope computation from slpmap. This is the CTX image P07_003621_1980_XI_18N133W which has a nice view of some craters in it. Here's the steps I took to generate the new phocube image:

1. ingest with mroctx2isis
2. attach spice with spiceinit and make sure it picks up a DEM (MOLA in this case)
3. run phocube to compute slopes
4. map project the slope band to an equirectangular projection

Then to produce the slpmap comparison I did the following

1. use map2map to get a subsection of MOLA in the same projection as the projected slope image
2. use slpmap to compute slops on the new MOLA sub-section

In this image you can see that the slopes from phocube line up nicely with slpmap. The top left is the projected phocube output, the top right image is the slpmap output, the bottom left is the MOLA sub-section, and the bottom right is a 2d plot comparing the phocube and slpmap output. The slpmap output looks very strange with huge variations and a strange gridding happening. I'm not sure why that's happening, but the phocube output looks like a less noisy version of it.

In this image you can see the slope aspect/azimuth comparison. This one's a bit stranger. It is roughly the same region and you can see that on the walls of the crater rims the two line up again. The slpmap output also has an odd grid effect for the slope aspect and the phocube output is like a less noisy version of it. On the "flat" ground, though, things go a bit wild. The slope in these regions is very nearly 0, so it's not strange that the azimuth angles vary wildly, but it varies in a different way than the slpmap output. Is that okay?

[lavalaz]

Very interesting! My guess is that the grid on the MOLA is from running slpmap at a spatial resolution higher than the grid spacing (you seem to be resolving the edges of the grids in the gridded data). It does look like the way slope and azimuth/aspect are being computed is consistent. The slope values look reasonable. I'm not sure I can validate the aspect values from just the screen shot, but the range is a bit odd - I would expect 0 to 360 (or -180 to +180, but 0-360 is more typical).

[jessemapel]

The resolution explanation makes sense, the map2map output is MUCH MUCH higher res than MOLA. Here's a different screenshot that shows the aspect value range. This is the small crater at the bottom of the image and you can see the jump from 0 to 360 at the bottom of the crater.

[rbeyer]

Yeah, Laz is right, the MOLA data is definitely oversampled, and the gridwork is because of that. However, phocube must be using the same DEM, just sampling it differently. If you map-project the MOLA, but don't alter its resolution, the output of slpmap on it should then be commensurate with what phocube is now doing. The "aspect" image might even look more similar.

[lavalaz]

This is all looking good. My only concern is that the plots seem to show some pixels have negative values for azimuth/aspect. The way it flips from 360 to 0 looks correct. And FYI, I can confirm that this ability to visualize the slope and aspect behind the local emission and incidence angles is really going to help me understand what is real versus noise in the data. Thank you!

[jessemapel]

I regenerated the slpmap output with the correct resolution and it looks much better. Unfortunately, there appears to be an issue with the slope azimuth. The slpmap output looks quite nice, but the phocube output is a mess. I am going to look back into what could be going on and try another path to computing the azimuth that's more direct.

[jessemapel]

I tinkered with this more over the weekend but couldn't get something that computed the slope azimuth correctly on the flat portions of my test image. I've updated the status on the PR and am going to leave this issue open until it can be worked on more.

[jessemapel]

The incomplete code that produced the above results can be found here: https://github.com/jessemapel/ISIS3/tree/phoslope