3

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?

Trenton McKinney
  • 56,955
  • 33
  • 144
  • 158
  • It's usual to ask new questions as new StackOverflow questions otherwise things get confusing, but in brief: (1) My guess: OpenGL; try pyvista, mayavi or another 3D viz package; (2) maybe switch the plotting order? (3) make sure you save before plt.show(); (4) I think it is. – Matt Hall Dec 01 '21 at 18:01
  • Thanks for your answer – and for the intro to Stackoverflow :o) I'm much more used to the MathWorks Support and Answers. There you can add follow-up questions that makes a lot more sense to me than opening another case. – Martin H. Trauth Dec 02 '21 at 07:57

2 Answers2

2

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

Matt Hall
  • 7,614
  • 1
  • 23
  • 36
1

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

amoodie
  • 323
  • 2
  • 10
  • Beautiful, Andrew, I'll buy that! Best, Martin – Martin H. Trauth Nov 29 '21 at 17:13
  • Hmm, but I think now you are switching LON and LAT... I think the difference you're seeing is that the OP did a `np.flipud()` on the data, and I didn't. It's confusing. Anyway, I updated my answer to reverse the `lat` linspace, and you inspired me to add the shading I was too lazy to do before. – Matt Hall Nov 29 '21 at 21:14
  • No worries, these digital terrain data are defined differently and we have to flipud, fliplr and tranpose them all the time, in addition to the confusing difference in indexing Python vs. MATLAB. – Martin H. Trauth Nov 30 '21 at 12:50