Solar module inverts shading when used with -Rxmin/xmax/ymin/ymax

[RunningWild89]

Running the solar module while specifying -R with coordinates rather than -Rg results in the shading of the wrong area. In other words, Instead of the areas away from the sun being dark, the areas NEAR the sun are dark.

gmt begin solar_test_2
	#gmt coast -JH0/30c -Rg -Dh -W -Bafg -A500 -Gtan
	gmt coast -JM0/30c -R-15/15/40/60 -Bafg -Dh -W -A500 -Gtan
	gmt solar -Td+d2020-11-04T16:30:00+z7 -Gnavy@80
	gmt solar -Tc+d2020-11-04T16:30:00+z7 -Gnavy@70
gmt end show

No error messages are produced.

Compare Outcomes

solar_test_1.pdf
solar_test_2.pdf

System information

• Operating system:
• GMT version 6.1 on windows. Script run on Git with bash.

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

[liamtoney]

I just discovered this bug as well. Here's another example (PyGMT). This code:

import pygmt
from obspy import UTCDateTime

now = UTCDateTime()  # Get current UTC time

fig = pygmt.Figure()

for region, projection, xshift in zip([(-180, 180, 60, 90), 'g'], ['S-150/90/5i', 'R-150/8i'], ['0', '7i']):

    fig.coast(shorelines=True, land='lightgray', area_thresh='0/0/1', region=region, projection=projection, xshift=xshift)

    # Use a context manager to call solar
    with pygmt.clib.Session() as session:
        session.call_module('solar', f'-Td+d{now.isoformat()} -Gnavyblue@60 -Ba60f15g30')
            
fig.show()

produces this:

The left-hand figure's shading is inverted.

Output of pygmt.show_versions()

PyGMT information:
  version: v0.3.0
System information:
  python: 3.8.8 | packaged by conda-forge | (default, Feb 20 2021, 16:12:38)  [Clang 11.0.1 ]
  executable: /opt/miniconda3/envs/gi-pygmt-2021/bin/python
  machine: macOS-10.16-x86_64-i386-64bit
Dependency information:
  numpy: 1.20.1
  pandas: 1.2.2
  xarray: 0.16.2
  netCDF4: 1.5.6
  packaging: 20.9
  ghostscript: 9.53.3
  gmt: 6.1.1
GMT library information:
  binary dir: /opt/miniconda3/envs/gi-pygmt-2021/bin
  cores: 8
  grid layout: rows
  library path: /opt/miniconda3/envs/gi-pygmt-2021/lib/libgmt.dylib
  padding: 2
  plugin dir: /opt/miniconda3/envs/gi-pygmt-2021/lib/gmt/plugins
  share dir: /opt/miniconda3/envs/gi-pygmt-2021/share/gmt
  version: 6.1.1

[PaulWessel]

@joa-quim, we need to fix this one I think. I am guessing the polygon is completely covering this box and we need to detect that.

[PaulWessel]

Looked at this some more and fixed a minor memory leak (#5101) unrelated to the failure. Here is a wider view to see what the global polygon looks like and what a small region looks like. Trying to run pssolar with the small outlined region and Mercator gives the wrong sense of fill. Here is the view centered on the antipole of the point that pssolar draws a small circle around:

In contrast, selecting the small red box with Mercator gives the inverse result:

The small circle routine takes the location of the sun at the given time, which is ~33E, 15S. So that should return a global polygon that is an oblique polar cap around that point. But that is the sun-lit side; we want the other side. I do not see anything in pssolar that specifies this, yet the painting for global maps is correct (but fails for smaller). Perhaps the bug is not that it fails for smaller regions, but that is is implemeneted backwards and that happens to work for global views due to some other bug? What do you think, @joa-quim ?

[PaulWessel]

Checking and finding that the clipping in the Mercator example above works great except it picks the wrong side. Trying to think of simple ways to defeat that. Given the nature of the polygons, it may work to split it into two closed polygon on either hemisphere (either in latitude or longitude), then fill/clip using the two (and draw optional outline with the original polygon).

[joa-quim]

So, it's not because we want the dark side of the polygon?

[PaulWessel]

We want the dark polygon since that is what we must fill. The code creates a global oblique small-circle polygon. This clearly gets filled correctly for global projections but flips (at least in the above case - not sure how general it is) for smaller regions. I was hoping that splitting along Greenwhich would work but not (yet at least). I also tried testing various spherical triangles by connecting sections of the line to the pole point and fill those, but gets in to wrapping longitudes at some point. So still a fail and may have to remain a fail for rc1.