Python: Why is plot_surface transparent?

Viewed 383

I tried to create a surface plot with Python 3 and now I am wondering why it is transparent? Any ideas? I expected it to look like the plot I created with MATLAB using the same data set ETOPO1. A second question, changing the aspect ratio isn't possible with plot_surface,right?

Best, Martin

import numpy as np
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

ETOPO1 = np.flipud(ETOPO1)

lon = np.arange(30,60+1/60,1/60)
lat = np.arange(-20,20+1/60,1/60)
LON,LAT = np.meshgrid(lon,lat)

fig, ax2 = plt.subplots(subplot_kw={"projection": "3d"})
ls = LightSource(270,45)
rgb = ls.shade(ETOPO1,
        cmap=cm.gist_earth,
        vert_exag=0.1,
        blend_mode='hsv')
ax2.plot_surface(LON,LAT,ETOPO1,
        rstride=1, cstride=1,
        linewidth=0,
        facecolors=rgb,
        antialiased=True,
        shade=True)
ax2.view_init(60, 20-90)
ax2.tick_params(axis='both', labelsize=6)
ax2.grid(False)
plt.show()

3D surface from python

3D surface from matlab

Edit 1 Dec 2021

Here are a few related questions:

  1. The code above needs 15 sec if I use rstride=5, cstride=5 but 318 sec = 5.3 min if I use the full resolution rstride=1, cstride=1. This is surprising as MATLAB needs 0.14 sec for the full resolution – why?

  2. I tried to add contours on top of the surface using

    v = np.array([500,1000,2000,3000])
    ax2.contour(LON,LAT,ETOPO1+1,
         levels=v,linewidths=0.3,
         colors='r',linestyles='solid')
    

    but, despite adding 1, 10, or 100 to ETOPO1, they are always hidden by the surface.

  3. I tried to save the figure using ETOPO1, I

    plt.savefig('etopo1_python.png',dpi=300)

    but got an empty image in the PNG file. Any ideas?

  4. Since

    antialiased=True

    causes transparency, the question arises, is it a bug?

2 Answers

According to this answer it's a known issue with the plot_surface() function.

This other answer suggests setting antialiased=False in plot_surface() and this worked for me on your data. I did not look closely to see if it caused any issues, however.

As shown below, you can adjust the aspect ratio by adding something like ax.set_box_aspect((1, 1, 0.1)) (below I compute the Z parameter to achieve some desired approximate vertical exaggeration, but you get the idea).

import numpy as np
from matplotlib import cm
from matplotlib.colors import LightSource
import matplotlib.pyplot as plt

url = 'http://141.89.112.21/wp-content/uploads/2021/11/etopo1_data_python.txt'
ETOPO1 = np.loadtxt(np.DataSource().open(url), skiprows=5)

h, w = ETOPO1.shape
lon = np.linspace( 30,  60, w)  # linspace recommended for non-integer intervals.
lat = np.linspace( 20, -20, h)  # Reverse this instead of flipping the array.
LON, LAT = np.meshgrid(lon, lat)

fig, ax = plt.subplots(figsize=(20, 15),
                       subplot_kw={"projection": "3d"})

ve = 200  # Approx. vertical exaggeration.
ax.set_box_aspect((1, 1, ve/1850))

rgb = LightSource(270, 45).shade(ETOPO1,
                                 cmap=cm.Blues_r,
                                 blend_mode='soft',
                                 vert_exag=ve/1850)

ax.plot_surface(LON, LAT, ETOPO1,
                rstride=2, cstride=2,
                facecolors=rgb,
                antialiased=False)

ax.view_init(20, -70)
ax.tick_params(axis='both', labelsize=6)
ax.grid(False)
plt.show()

This produces:

The result of the code

Building on @kwinkunks and your attempt with shade gets pretty close to what you're after. Note: I think LAT, LON were switched out of meshgrid in kwinkunks answer.

h, w = ETOPO1.shape
lon = np.linspace( 30, 60, w)
lat = np.linspace(-20, 20, h)
LAT, LON = np.meshgrid(lon, lat)

fig, ax = plt.subplots(
    figsize=(10, 10),
    subplot_kw={"projection": "3d"})

ve = 100  # Approx. vertical exaggeration.
ax.set_box_aspect((1, 1, ve/1850))

ls = LightSource(270,45)
cs = ls.shade(
    ETOPO1, cmap=cm.gist_earth,
    vert_exag=ve/1850)
ax.plot_surface(
    LON, LAT, ETOPO1,
    rstride=5, cstride=5,
    facecolors=cs)

ax.view_init(20, 15)

example

Related