crater shape/orientation with spice

[steo85it]

While making experiments with your libraries at the Moon's south pole, I noticed a (possible) issue when using lsp_spice. This should be the interior of a crater (see elevation map of the region below), with a smaller crater within it even projecting a shadow of its own (which is clearly not the expected behavior...).

I think that this could be due to the fact of being in the southern hemisphere and the transformation to cartesian, which the code ending up "confusing" a crater for a hill.

(above is the elevation from the Moon mean radius, while below is Z ~ signed distance from the Moon center)

You can check it by setting

i0, i1 = 2425, 2663  # 600, 700
j0, j1 = 1912, 2172  # 650, 725

as indexes to select a region from the whole south_pole grd "LDEM_80S_150M_adjusted.grd" in lsp_make_shape_model and running the lsp_shape.py example for the 6-Sep-2021. Since the region is a bit large, you can maybe subsample the map a bit by adding

IJ = IJ[::3]
V = np.array([x.ravel(), y.ravel(), z.ravel()]).T
V = V[::3]

to lsp_make_shape_model.py.

I tested several solutions/options:

• checked what happens if lat0 is set to 90 (North Pole) -> the illumination changes a bit (seasons?), but the craters persists behaving as hills
• using a (stereo-) projected map (instead of converting it to a cartesian point cloud) works in preserving the "crater behavior", but I don't know then how to get the actual illumination at the South Pole at the Moon (instead of at its center). At least in some tests, seasons seem to be opposite than expected (Sep 21 should be winter and Mar 21 Summer at the Moon South Pole, but I get a Sun with a higher elevation at these craters).
• following this test, I thought that in the end what counts is the sun direction at a given location, so that sun_dirs in lsp_spice.py could be tweaked to include the cartesian position of the crater at the South Pole. This doesn't seem to impact the result (the angular difference given by the radius of the Moon is too small? Or is there something I'm missing in how that sun_dir is then used by get_direct_irradiance?
• I also thought that the N(ormal) vector should be responsible of telling the code where the "outside" of the planet is, but that looks fine (~ 0,0,-1) in the non-working tests. Still, in some tests I noticed that introducing/removing the test

  N[(N*P).sum(1) < 0] *= -1

  can reproduce the same behavior of "hill vs crater".

Any hint here would be appreciated! :)

[sampotter]

OK, had a chance to start looking at this. We may need to go back and forth a bit to figure out what's going on.

From your post, it looks like there's some confusion (probably mine :^) ) about how to set sun_dirs.

My assumption is the following:

• possun is the position of the sun relative to the "position of the moon", where I'm assuming "position of the moon" = its centroid
• so, to compute sun_dirs, just normalize possun to get the correct direction

In the MATLAB SPICE script Erwan sent me, he did this using trigonometry, which I just translated directly. But the result is the same as normalizing. I verified that this is the case, and made the simplification (see my latest commits, in particular 84712bf7e31bd8ff4d210ced09b71e972fae4d0f).

Is my understanding of this correct? Let's get this squared away first...

[steo85it]

Great, thanks! This looks correct and, at least in the Moon case, it should be ok also if we are/our target is on the surface of the Moon instead of at its center of mass (distance_Sun_Moon >> R_Moon). And your simplification looks good, too (in principle, the second element of spkpos output should be the source to target norm and could be directly used, even though sometimes small additional/unwanted corrections could be included).

[sampotter]

Ah, I think I figured it out. I was passing the arguments to TrimeshShadeModel in lsp_spice.py. See the change in 83ce63f7bc6ece3ded71effaa0efbdf5c937307a. I'm testing this now to see if it fixes things, but I'm pretty sure it should.

In TrimeshShadeModel's constructor, if you just pass in V and F (which is what I was effectively doing), you compute N automatically by normalizing the cross product of two edges in each triangle indexed by F. This ends up with upside down normals at the LSP. Since I had the argument order wrong, the face normals were being computed automatically, resulting in the behavior that you observed.

[steo85it]

Ah, that would make sense! I was indeed playing with N arrays too, but never checked the order... let me check that out! Great job if it works!

[steo85it]

It works for me now, thanks!!