How to make an inset plot with mollweide projection?

Question

I want to make a skymap using the Mollweide projection for a main set of axes and for an inset axes. This is easy for the main axes but not for the inset. I've tried a few different things but it doesn't work for the inset. Please help!

Here you can find the latitude and longitude data, and here you can find the sky location probability density data.

First, I make the main plot:

xmin = min(l)
xmax = max(l)
ymin = min(b)
ymax = max(b)
X, Y = np.meshgrid(np.linspace(xmin, xmax, 100), np.linspace(ymin, ymax, 100))

mpl.rcParams["text.usetex"] = True
fig = plt.figure(1)
fig.set_figheight(8)
fig.set_figwidth(8)

ax = plt.axes(projection='mollweide')
ax.grid()
# skypost is the sky location probability-density data accessible above
plt.contour(X, Y, skypost, colors='blue', levels=[5, 50, 95])

which works fine. Next, I define the inset axes and plot the contours, however there seems to be no way that completely works for this. What I want is for the inset to zoom-in on the contours while keeping the mollweide projection. I've tried to do as the example on ligo.skymaps, i.e.,

axesinset = plt.axes(
    [0.0, 0.2, 0.25, 0.25], 
    projection='astro degrees zoom', 
    center='110d +20d', 
    radius='10 deg' )
plt.sca(axesinset)
axesinset.contour(X, Y, skypost, colors='blue', levels=[5, 50, 95])
axesinset.grid()

but this doesn't work since the contours don't even appear! I don't understand why they don't appear. I also do not understand why the x-axis of the inset is backwards?

Instead, I've tried just plotting a new mollweide projection in the inset, and restricting the xlim and ylim, but it says these options are not supported for the mollweide projection. Is there a way around this to restrict the axes limits?

Lastly, I've tried just doing a regular inset without the mollweide, which works, however the shape of the contours are distorted relative to the contours on the main mollweide plot which is physically relevant for my case. So this is very sub-optimal.

Any suggestions and advice are greatly appreciated.

Answer 1

I'm a bit late to the party, but I thought its worth mentioning that I've created a nice inset-map functionality for EOmaps...

It lets you create inset-maps in arbitrary projections and you can add whatever features you want!

from eomaps import Maps

m = Maps(Maps.CRS.Mollweide())
m.add_feature.preset.coastline()

# create a rectangular inset-map that shows a 5 degree rectangle
# centered around a given point
inset = m.new_inset_map(xy=(6, 43), xy_crs=4326,
                        radius=5, radius_crs=4326,
                        inset_crs=Maps.CRS.Mollweide(),
                        shape="rectangles")
inset.add_feature.preset.coastline()
inset.add_feature.preset.ocean()
inset.add_feature.cultural_10m.urban_areas(fc="r", ec="none")

m.apply_layout(
    {'0_map': [0.01, 0.17333, 0.98, 0.65333],
     '1_map': [0.05, 0.11667, 0.43341, 0.76667]})

Answer 2

To have the axis in the correct way, you can rotate the subplot by using rotate.

Concerning the fact that your contour are not shown, it is probably because you have to add the transform keyword. If you don't specify it, it is plotted in pixel coordinates by default (https://docs.astropy.org/en/stable/visualization/wcsaxes/overlays.html).

The example below shows that the desired point (in blue) is obtained by adding ax.get_transform("world"). The blue and green points are in the lower right corner because of the rotate.

I guess that it should be the same for contour.

ax = plt.subplot(111, projection='geo degrees zoom',
              center="0d - 0d", radius='10 deg', rotate='180 deg')
ax.grid()
ax.set_xlabel(r"$\phi \, [deg]$")
ax.set_ylabel(r"$\theta \, [deg]$")

ax.scatter(0,0, color = "blue")
ax.scatter(100,0, color = "green")
ax.scatter(0,0, color = "red", transform = ax.get_transform("world"))