4D contour plot using .nc file

Question

I am trying to plot a 4D surface plot from netcdf data which has 4 dimensions: time, lat, long, and lev for 5 variables (DU001, DU002...005) (sample data). I have to plot the first variable DU001 vs lat, long, and levels (72 levels) such that the x-axis is lat, the y-axis is long, the z-axis is levels, and the DU001 will be represented with color. So far I have tried the below code but I am getting only one surface in my plot .

I think it is only taking one level. How to correct it?

import xarray as xr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import cm
path= 'D:\\DATA\\2015\\test_data'# Open the NetCDF file
data = xr.open_dataset('D:\\DATA\\2015\\test_data\\MERRA2_400.inst3_3d_aer_Nv.20150515.SUB.nc')
# Select the DU01 variable and the lat, long, and lev dimensions
lat = data['lat']
lon = data['lon']
lev = data['lev']
DMR = data['DU001'] 
# Reshape the data
du001_2d = DMR[:, :, :].squeeze()
dmr_values = du001_2d.values.squeeze()
# Create meshgrid for coordinates
lon_2d, lat_2d = np.meshgrid(lon, lat)
# Create the 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the surface for each level
for i, level in enumerate(lev):
    ax.plot_surface(lon_2d, lat_2d, dmr_values[i], cmap='viridis')
    # Set labels and title
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    ax.set_zlabel('Level')
    ax.set_title('3D Plot of Dust Mixing Ratio')

    # Set the z-limits based on the valid range of the lev array
    ax.set_zlim(lev[0], lev[71])  # Assuming lev is a 1D array

# Display the plot
plt.show()

I don't know where I am going wrong. I am very new to Python. Any help would be appreciated

sample of what I am trying to plot is

Have you tried `ax.plot_surface(lon_2d, lat_2d, level, facecolors=dmr_values[i], cmap='viridis')`? You may need to remap `dmr_values[i]` to be in a format for the color methods, see: https://stackoverflow.com/questions/32413371/color-of-a-3d-surface-plot-in-python. — jared, Jun 09 '23 at 05:53
Yes, @jared I have tried. But I am getting an error ValueError: Argument Z must be 2-dimensional. — shravani banerjee, Jun 09 '23 at 06:22

Stef · Accepted Answer · 2023-06-12T16:35:11.033

1

You plot all dmr_values one on the other, i.e. what you see is just dmr_values[-1]. The maximum height of a layer is about 6e-8, that's why it appears as just one flat surface when using a z-axis range of 0 - 72.

If you want to plot 72 colored layers, you need to provide 72 flat surfaces and color them yourself:

ax.plot_surface(lon_2d,
                lat_2d,
                np.full_like(lat_2d, level),
                facecolors=cm.ScalarMappable(cmap='viridis').to_rgba(dmr_values[i]),
                shade=False)

Update for comment below: If you want to interpolate you could use 3 filled contour plots for the 3 visible surfaces:

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

path = 'MERRA2_400.inst3_3d_aer_Nv.20150515.SUB.nc'
data = xr.open_dataset(path)

lat = data['lat']
lon = data['lon']
lev = data['lev']
DMR = data['DU001']

fig, ax = plt.subplots(subplot_kw=dict(projection="3d"))

x, y, z = np.meshgrid(lon, lat, lev, indexing='ij')

du = np.swapaxes(DMR[:, :, :].squeeze().values, 0, -1)
kw = {
    'vmin': du.min(),
    'vmax': du.max(),
    'levels': np.linspace(du.min(), du.max(), 20),
}

_ = ax.contourf(
    x[:, :, 0], y[:, :, 0], du[:, :, -1],
    zdir='z', offset=z.max(), **kw
)

_ = ax.contourf(
    x[:, 0, :], du[:, 0, :], z[:, 0, :],
    zdir='y', offset=y.min(), **kw
)

c = ax.contourf(
    du[-1, :, :], y[0, :, :], z[0, :, :],
    zdir='x', offset=x.max(), **kw
)

xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
zmin, zmax = z.min(), z.max()
ax.set(xlim=[xmin, xmax], ylim=[ymin, ymax], zlim=[zmin, zmax],
       xlabel='Longitude', ylabel='Latitude', zlabel='Level')
fig.colorbar(c, ax=ax, shrink=0.7)

edited Jun 12 '23 at 16:35

answered Jun 09 '23 at 12:18

Stef

28,728
2
24
52

Dear @Stef, I got this plot. Thanks. But, I wanted a smooth plot just like an interpolated one. Is there any way to make my plot not look like a plot with stacked images? I want to show the variation of DMR with levels. – shravani banerjee Jun 10 '23 at 10:48
I think the interpolated one is not interpolating the lower layer or levels. As seen from the stacked plot earlier and the interpolated ones, the difference is in the lower layers. – shravani banerjee Jun 15 '23 at 06:39
The interpolated version is correct, it normalizes colors over the entire volume. The stacked layers version normalizes colors per layer, if you normalize over all layers (`cm.ScalarMappable(mpl.colors.Normalize(vmin=dmr_values.min(), vmax=dmr_values.max()), cmap='viridis')`) you'll get the same result. – Stef Jun 15 '23 at 06:45
What if, I slice my data between a short range of levels as well as a shorter range of lat and long, then a good amount of information can be seen in the plot? I think then the interpolation will be done within a smaller volume instead of a larger volume. – shravani banerjee Jun 15 '23 at 07:08
the interpolation is always down on all data provided, what I'm talking of is the color coding. You can of course set `vmin` and `vmax` to arbitrary values to highlight special ranges. If your data vary to a great extent you may want to use a logarithmic norm (or symlog if data includes zeros). – Stef Jun 15 '23 at 07:19
But, the contourf only plots on specific planes, it does not plot every combination of lon, lat, and lev values. Is there any way to interpolate my data in the overall volume like cubic or nearest anything? – shravani banerjee Jun 15 '23 at 08:46
In one plot you can only show one set of planes, not every combination simultanuously. But of course you can use any (existing) lon, lat and level values for the planes, just change the indices in x, y, z and du and set the offset accordingly. For values between the existing values you'll need to interpolte (numpy, scipy). – Stef Jun 15 '23 at 08:58
I am trying to plot a smooth 3D plot. A sample of a plot that I actually need to plot is updated in my question now. Is there any way to plot it? – shravani banerjee Jun 15 '23 at 10:41
how to slice the lat and long and then plot with the earlier code, which is not the interpolated one? (lat can be between 20-30 and long between 70-80). – shravani banerjee Jun 19 '23 at 08:47
just select the required data `data = data.sel(lat=slice(20, 30), lon=slice(70, 80))` – Stef Jun 19 '23 at 09:06

4D contour plot using .nc file

1 Answers1