Search code examples
pythonscatter-plotmatplotlib-basemapz-order

Basemap scatterplot with zorder determined by size

I'm looking to make a basemap with a scatterplot on top with the zorder of the points determined by the size of each individual point, so that no point ends up totally covered by another. (The ideal end result would look like a bullseye.)

Let's say I have the following code:

Ca_data = array([0.088, 0.094, 0.097, 0.126, 0.112, 0.092, 0.076, 0.105])
SO4_data = array([0.109, 0.001, 0.001, 0.007, 0.214, 0.005, 0.008, 0.559])
longitude = linspace(-101, -100, 8)
latitude = linspace(34.5, 35, 8)

m=Basemap(llcrnrlon=-101,llcrnrlat=34.5,urcrnrlon=-100,urcrnrlat=35,resolution='c', epsg=4326)
m.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D', xpixels=1500, ypixels=1500, zorder=1)
m.scatter(longitude, latitude, latlon=True, s=6000*Ca_data,c='r',marker="o",label='Ca')
m.scatter(longitude, latitude, latlon=True, s=6000*SO4_data,c='b',marker="o",label='SO4')
plt.show()

As it is, anywhere SO4 is larger than Ca, I would only see SO4. I've considered going in and adding a zorder to each line, but I don't think that would work well since I have several more elements to add in with the same issues. Any ideas?

Solution

  • The only way I can see how to do this is to plot both your datasets at once. Otherwise there's no way to not plot the second dataset on top of the first.

    Most of this can be done easily, the only issue is the label kwarg of scatter. We can give elementwise (array_like) sizes and colours for scatter, but we can't do that with labels. Still, as the visualization of the data is much more important, I'd go down this route, hacking up label-related issues (a legend, mostly) later:

    import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

Ca_data = np.array([0.088, 0.094, 0.097, 0.126, 0.112, 0.092, 0.076, 0.105])
SO4_data = np.array([0.109, 0.001, 0.001, 0.007, 0.214, 0.005, 0.008, 0.559])
longitude = np.linspace(-101, -100, 8)
latitude = np.linspace(34.5, 35, 8)

# this is the new part: concatenate all the data    
plotlon = np.tile(longitude,2)
plotlat = np.tile(latitude,2)
plotdat = np.concatenate((6000*Ca_data,6000*SO4_data)) # sizes
cdat = np.repeat(('r','b'),longitude.size)             # colors

# determine reverse sorting order
inds = np.argsort(plotdat)[::-1]

plt.figure()
m = Basemap(llcrnrlon=-101,llcrnrlat=34.5,urcrnrlon=-100,urcrnrlat=35,resolution='c', epsg=4326)
m.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D', xpixels=1500, ypixels=1500, zorder=1)
# use a single scatter() call, with ordered-concatenated data
m.scatter(plotlon[inds], plotlat[inds], latlon=True, s=plotdat[inds],c=cdat[inds],marker="o") # label has been removed!

plt.show()

    result

    Generalization to n datasets is straightforward: you need to tile your latitudes and longitudes n times instead of 2, and you need a colour for each dataset that you then repeat.