Plotting 4D data as contours

Question

I want to create layered contour plots using (x,y,z) coordinates and a fourth dimension, represented by color and intensity. Each layer-related data is in the CSV file (data) with lat, lon, levels, and dust_mixing_ratio columns. How to create and superimpose all the heatmaps of each level along the Z-axis? Or is there any other way to plot my data in 4D like the sample image? So far I have tried to plot using the following code as mentioned in Plot 4D data as layered heatmaps in Python but I am unable to understand how to proceed further?

def grids_maker(filepath):
# Get the data
    df = pd.read_csv("D:\DATA\2015\MIXING_RATIOe.csv", sep=' ')
    # Make things more legible
    xy = df[['lat', 'lon']]
    x  = xy.lat
    y  = xy.lon
    z  = df.levels
    g  = df.dust_mixing_ratio
    reso_x = reso_y = 50
    interp = 'cubic' # or 'nearest' or 'linear'
    # Convert the 4d-space's dimensions into grids
    grid_x, grid_y = np.mgrid[
        x.min():x.max():1j*reso_x,
        y.min():y.max():1j*reso_y
    ]

    grid_z = si.griddata(
        xy, z.values,
        (grid_x, grid_y),
        method=interp
    )

    grid_g = si.griddata(
        xy, g.values,
        (grid_x, grid_y),
        method=interp
    )

    return {
        'x' : grid_x,
        'y' : grid_y,
        'z' : grid_z,
        'g' : grid_g,
    }

Answer 1

Based on the answer that you linked to, you could do something like:

import matplotlib.pyplot as plt
import scipy.interpolate as si
from matplotlib import cm
import pandas as pd
import numpy as np

# read in data
data = pd.read_csv("MIXING_RATIOe.csv")

# create colormap
scam = plt.cm.ScalarMappable(
    norm=cm.colors.Normalize(data["dust_mixing_ratio"].min(), data["dust_mixing_ratio"].max()),
    cmap="jet",
)

# make the figure
fig = plt.figure()
ax = fig.add_subplot(projection="3d")

# get longitude/latitude grids
lonmin = data["lon"].min()
lonmax = data["lon"].max()
latmin = data["lat"].min()
latmax = data["lat"].max()

# use 100 points across latitude and longitude (chanhe this as you want)
lats = np.linspace(latmin, latmax, 100)
lons = np.linspace(lonmin, lonmax, 100)

gridlat, gridlon = np.meshgrid(lats, lons)

# transparancy
alpha = 0.25

for level, values in data.groupby(["levels"]):
    lon = values.lon
    lat = values.lat
    
    latlon = values[["lat", "lon"]]
    dust = values["dust_mixing_ratio"]

    # interpolate dust mixing ratios onto the grid
    griddust = si.griddata(latlon, dust.values, (gridlat, gridlon), method="cubic")

    scam.set_array([])
    fc = scam.to_rgba(griddust)

    # set transparancy
    fc[:, :, -1] = alpha

    ax.plot_surface(
        gridlat,
        gridlon,
        np.full_like(gridlat, level[0]),
        facecolors=fc,
        antialiased=True,
    )

ax.set_xlabel("Latitude")
ax.set_ylabel("Longitude")
ax.set_zlabel("Level")

plt.show()

which produces:

Obviously, the levels here are very close together, so it's hard to see individual layers and what's going on inside, but you could expand the z-axis a bit (although with 72 different levels it's always going to be tricky to view them all). You could make an animation that adds a level at a time.

Update

If you want to limit the ranges (set_x/y/zlim does not produce attractive results), you can just limit the grid ranges and only plot for a certain range of levels. Applying the ranges from the comment can be done with:

...

# explicitly set the grid ranges
lats = np.linspace(20, 30, 100)
lons = np.linspace(70, 80, 100)

...

for level, values in data.groupby(["levels"]):
    # only plot levels between 42 and 72
    if level[0] < 42 or level[0] > 72:
        continue

    lon = values.lon
    lat = values.lat
    ...