Scatter plot on Robinson projection chops off part of map

[AgilentGCMS]

Description

I am trying to make a scatter plot on a map using cartopy. I calculate the figure size based on the map aspect ratio in order to completely fill up the figure. If I omit the scatter call, this works fine. However, with the scatter call, a portion of the globe to the right gets chopped off.

Steps to reproduce

1. Download the code and a sample data file and run the code
2. The map looks like this, with a little bit chopped off to the right
3. Comment the line sc = plot_ax.scatter in the code and run again. This time, the map looks like this. It shows the whole globe, no chopping off.

Environment summary

• Python 3.9.10 (installed via macports)
• Matplotlib 3.5.1 with MacOSX backend
• Numpy 1.22.0
• Cartopy 0.19.0.post1
• macOS Big Sur 11.6.2

Other comment

I do not think this has anything to do with me trying to fill up the figure area. If I replace the line

plot_ax = plt.axes([0., 0., 1., 1.], projection=map_proj)

with

plot_ax = plt.axes([0., 0., 0.9, 0.9], projection=map_proj)

the problem still occurs, as evidenced here. So I think the problem is with the scatter call.

[AgilentGCMS]

Update

This problem also occurs with the EckertIV projection but not with the EckertIII projection. Maybe that's a clue?

[greglucas]

Could you please copy/paste the code here using triple backticks in a code block, and you can copy-paste images into the issue as well. It should be a minimal self-contained example. If you are only doing a scatter plot, try to make it with fake data and not require an external data file if possible.

[mathause]

You may have to call ax.set_global() after scatter as it restricts the axes to the extent of the data (and I think matplotlib does the same).

[AgilentGCMS]

@greglucas Here's a MWE:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import pyplot as plt
import numpy as np

# get the data
lats, lons, nobs = np.loadtxt('ch4_data.txt', unpack=True)

lon_center = 0.0
map_proj = ccrs.Robinson(central_longitude=lon_center)
data_proj = ccrs.PlateCarree(central_longitude=lon_center)

fig = plt.figure()
plot_ax = plt.axes(projection=map_proj)
plot_ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.7, edgecolor='k')
plot_ax.add_feature(cfeature.LAKES.with_scale('50m'), linewidth=0.5, edgecolor='k', facecolor='none')
plot_ax.gridlines(xlocs=np.arange(-180.,181.,30.), ylocs=np.arange(-90.,91.,30.), linewidth=0.5, linestyle=(0,(5,2)), color='0.25')

sc = plot_ax.scatter(lons, lats, s=0.01*nobs, c='0.75', edgecolors='0.25', marker='o', alpha=0.5, linewidths=1.0, transform=data_proj)

plt.show()

Unfortunately, I can't provide a MWE without the specific data. The data file is a three column text file with 247 lines, so I could have just made the numpy array inline (within the code), but I don't see what purpose it would serve except make the code long. Perhaps there are some data points in there that cause this problem, but then again, that doesn't help me, because that's the data I need to plot. As far as I can tell, all those are bona fide latitudes and longitudes.

[greglucas]

I believe @mathause has the solution. I assume that using these lats and lons shows the same issue you're seeing?

lats = np.linspace(-90, 90)
lons = np.linspace(-180, 180)

You probably don't have data in the region of ((-180, 0), (180, 0)). Try adding a scatter with those points and the extent should grow to include those bounds then.

[AgilentGCMS]

@mathause Indeed, ax.set_global() did the trick. However, this is not default matplotlib behavior for scatter. E.g., the following code produces this plot:

fig = plt.figure()
ax = plt.axes()
ax.set_xlim(-3.,3.)
ax.set_ylim(-4.,4.)
Xarr = 1.6*(np.random.random(10)-0.5)
Yarr = 2.0*(np.random.random(10)-0.5)
ax.scatter(Xarr, Yarr, s=100, c='0.75', edgecolors='0.25', marker='o', alpha=0.5, linewidth=1.0)
plt.show()

That is, a call to scatter does not restrict existing axes to the extent of the data. Why would cartopy do this? At least for global projections set_global() offers a way around, but I can see this being a problem for regional projections.

I can see how this extent recalculation would be a desired feature for simple X-Y plots such as

ax = plt.axes()
ax.scatter(X, Y, ...)

because the user has not indicated any preferences for the required extents. However, when one plots data on a map, one usually spends more time setting up the map projection and map features before plotting the data. In that case, shouldn't the assumption be that the user has already decided on the axes extents they want?

[AgilentGCMS]

@greglucas Indeed, @mathause has the right solution. However, I'm wondering if automatically resizing the axes is the correct thing to do when plotting data on a map.

[greglucas]

It should only expand the limits I believe. If you set the limits before the scatter, it does respect those limits as far as I can tell. Additionally, you can disable autoscaling of axes.

plot_ax.set_xlim(-1e6, 1e6)
plot_ax.set_ylim(-1e6, 1e6)
# plot_ax.autoscale(False)

[AgilentGCMS]

@greglucas OK, understood. However, in the case of a scatter plot on a map, it actually contracted the limits. Also, turning off autoscaling does not work for cartopy axes. If I turn of autoscaling just before calling scatter, like so:

plot_ax.autoscale(False)
sc = plot_ax.scatter(lons, lats, s=0.01*nobs, c='0.75', edgecolors='0.25', marker='o', alpha=0.5, linewidths=1.0, transform=data_proj)

then it produces this empty plot. The point I'm trying to make is that the solution suggested by @mathause works, but cartopy shouldn't be resizing the map axes like it currently does. Unless there's a good reason why it should, which I don't understand at the moment.

[CommonClimate]

Weighing in to support everything @AgilentGCMS said in this thread. I am glad ax.set_global() could solve the issue, but it really shouldn't be an issue in the first place! As a fellow open-source developer of scientific software, my golden rule is that if I can't rationalize defaults to my constituents, then they are the wrong defaults. I hope this post goes a small way to making better defaults for Cartopy. In any case, I am grateful to all the time and effort put into this project, which is tremendously useful to the Python geoscience community and beyond.

[rcomer]

I think what Cartopy is doing here is completely consistent with Matplotlib. If you create a new axes in Matplotlib, you get default x- and y-limits of [0, 1]:

In [1]: import matplotlib.pyplot as plt

In [2]: fig, ax = plt.subplots()

If the limits of the first artist you add are smaller than this default, the limits will shrink:

In [3]: ax.scatter([0.5, 0.6], [0.5, 0.6])
Out[3]: <matplotlib.collections.PathCollection at 0x7f4d2df97df0>

The OP wants a global map because they have points all over the world but if all the points were over, say, Australasia, they would likely want a different extent. So I think "have the extent fit the data" is a reasonable default.