I ended up making sample data and worked on both the 2D and the 3D case. I'll start with the 2D case that you already have working because the extension to the 3D case is then very simple.
2D
First, let's create some random sample data. Note that I import everything I need for later here
import numpy as np
import matplotlib.pyplot as plt
import cartopy
from cartopy.crs import PlateCarree
from matplotlib.colors import Normalize
def create2Ddata():
'''Makes some random data'''
N = 2000
lat = 10 * np.random.rand(N) + 40
lon = 25 * np.random.rand(N) - 80
z = np.sin(4*np.pi*lat/180.0*np.pi) + np.cos(8*np.pi*lon/180.0*np.pi)
return lat, lon, z
# Create Data
lat, lon, z = create2Ddata()
This will serve as some random, scattered, geospatial, data that we want to plot using the histogram function. The next step is then to both create bins that make sense followed by actually binning.
def make2dhist(lon, lat, z, latbins, lonbins):
'''Takes the inputs and creates 2D histogram'''
zz, _, _ = np.histogram2d(lon, lat, bins=(
lonbins, latbins), weights=z, normed=False)
counts, _, _ = np.histogram2d(lon, lat, bins=(lonbins, latbins))\
# Workaround for zero count values to not divide by zero.
# Where counts == 0, zi = 0, else zi = zz/counts
zi = np.zeros_like(zz)
zi[counts.astype(bool)] = zz[counts.astype(bool)] / \
counts[counts.astype(bool)]
zi = np.ma.masked_equal(zi, 0)
return lonbins, latbins, zi
# Make bins
latbins = np.linspace(np.min(lat), np.max(lat), 75)
lonbins = np.linspace(np.min(lon), np.max(lon), 75)
# Bin the data
_, _, zi = make2dhist(lon, lat, z, latbins, lonbins)
Then, we plot both the scattered data and the binned data as follows.
def plotmap():
'''background map plotting'''
ax = plt.gca()
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='None',
linewidth=0.5, facecolor=(0.8, 0.8, 0.8))
ax.spines['geo'].set_linewidth(0.75)
fig = plt.figure()
# Just plot the scattered data
ax = plt.subplot(211, projection=PlateCarree())
plotmap()
plt.scatter(lon, lat, s=7, c=z, cmap='rainbow')
# Plot the binned 2D data
ax = plt.subplot(212, projection=PlateCarree())
plotmap()
plt.pcolormesh(
lonbins, latbins, zi.T, shading='auto', transform=PlateCarree(),
cmap='rainbow')
plt.show()
Figure 2D, not allowed to embed figures yet...
At the top, the scattered data, at the bottom the binned data.
3D
Let's continue with the 3D case. Again, let's create some random scattered data that varies in time:
def create3Ddata():
''' Make random 3D data '''
N = 8000
lat = 10 * np.random.rand(N) + 40
lon = 25 * np.random.rand(N) - 80
t = 10 * np.random.rand(N)
# Linearly changes sign of the cos+sin wavefield
z = (t/5 - 1) * (np.sin(2*2*np.pi*lat/180.0*np.pi)
+ np.cos(4*2*np.pi*lon/180.0*np.pi))
return lat, lon, t, z
# Create Data
lat, lon, t, z = create3Ddata()
Now, instead of using histogram2d here, we will use histogramdd, which is just the N-dimensional version of the same function.
def make3dhist(lon, lat, t, z, latbins, lonbins, tbins):
'''Takes the inputs and creates 3D histogram just as the 2D histogram
function'''
zz, _ = np.histogramdd(
np.vstack((lon, lat, t)).T,
bins=(lonbins, latbins, tbins),
weights=z, normed=False)
counts, _ = np.histogramdd(
np.vstack((lon, lat, t)).T,
bins=(lonbins, latbins, tbins))
# Workaround for zero count values tto not get an error.
# Where counts == 0, zi = 0, else zi = zz/counts
zi = np.zeros_like(zz)
zi[counts.astype(bool)] = zz[counts.astype(bool)] / \
counts[counts.astype(bool)]
zi = np.ma.masked_equal(zi, 0)
return lonbins, latbins, tbins, zi
# Create bins
latbins = np.linspace(np.min(lat), np.max(lat), 75)
lonbins = np.linspace(np.min(lon), np.max(lon), 75)
tbins = np.linspace(np.min(t), np.max(t), 5)
# Bin the data
_, _, _, zi = make3dhist(lon, lat, t, z, latbins, lonbins, tbins)
Finally, we plot both the scattered data and the binned data side by side in respective time bins. Note the normalization that is used to make sure variations in time are easily observed.
Note that there are three loops (I could have put them in a single one, but this is nicer for readability).
- The first loop bins the data in time and plots the binned data in one slice each.
- The second loop bins the data in time, then in space using the 2D histogram function from earlier, and plots a slice for each time bin.
- The third function uses the already 3D binned data from above and plots the slices as well by accessing slices in the 3D matrix.
# Normalize the colors so that variations in time are easily seen
norm = Normalize(vmin=-1.0, vmax=1.0)
fig = plt.figure(figsize=(12, 10))
# The scattered data in time bins
# Left column
for i in range(4):
ax = plt.subplot(4, 3, 3*i + 1, projection=PlateCarree())
plotmap()
# Find points in time bins
pos = np.where((tbins[i] < t) & (t < tbins[i+1]))
# Plot scatter points
plt.title(f'{tbins[i]:0.2f} < t < {tbins[i+1]:0.2f}')
plt.scatter(lon[pos], lat[pos], c=z[pos], s=7, cmap='rainbow', norm=norm)
plt.colorbar(orientation='horizontal', pad=0.0)
# Center column
for i in range(4):
ax = plt.subplot(4, 3, 3*i + 2, projection=PlateCarree())
plotmap()
plt.title(f'{tbins[i]:0.2f} < t < {tbins[i+1]:0.2f}')
# Find data points in time bins
pos = np.where((tbins[i] < t) & (t <= tbins[i+1]))
# Bin the data for each time bin separately
_, _, zt = make2dhist(lon[pos], lat[pos], z[pos], latbins, lonbins)
plt.pcolormesh(
lonbins, latbins, zt.T, shading='auto', transform=PlateCarree(),
cmap='rainbow', norm=norm)
plt.colorbar(orientation='horizontal', pad=0.0)
# Right column
for i in range(4):
ax = plt.subplot(4, 3, 3*i + 3, projection=PlateCarree())
plotmap()
plt.title(f'{tbins[i]:0.2f} < t < {tbins[i+1]:0.2f}')
plt.pcolormesh(
lonbins, latbins, zi[:, :, i].T, shading='auto', transform=PlateCarree(),
cmap='rainbow', norm=norm)
plt.colorbar(orientation='horizontal', pad=0.0)
plt.show()
Figure 3D, not allowed to embed figures yet...
In the left column, the scattered, random, geospatial data, where the titles indicate the bins. In the center column, the 2D histograms using "by hand" time-binned data. In the right column, the slices that were binned using a 3D histogram.
As expected center and right columns show the exact same thing.
Hope this solves your problem.
