{"id":2523,"date":"2025-12-04T07:00:54","date_gmt":"2025-12-04T13:00:54","guid":{"rendered":"https:\/\/wollen.org\/blog\/?p=2523"},"modified":"2026-02-01T08:14:38","modified_gmt":"2026-02-01T14:14:38","slug":"matterhorn-3d","status":"publish","type":"post","link":"https:\/\/wollen.org\/blog\/2025\/12\/matterhorn-3d\/","title":{"rendered":"The Matterhorn in 3D"},"content":{"rendered":"

I’ve wanted to post a Matplotlib 3D visualization for a while, but sharing it on the web is difficult. It’s easy to save 2D plots and embed them on this page. For a 3D plot, it works better when you can do plt.show()<\/code> and explore it through the figure window.<\/p>\n

I think I found a good subject and a satisfying way to present it. The plan is to do plt.savefig()<\/code> like usual, and then create an animation that rotates 360\u00b0 around the plot.<\/p>\n

The Matterhorn<\/a> is a famous mountain peak in the Alps, just north of the Swiss-Italian border. Its summit rests at 14,692 feet above sea level. That makes it only the 12th-tallest mountain in the Alps, but its four steep faces that tower over the surrounding area make it one of the most recognizable in the world.<\/p>\n

The Swiss Federal Office of Topography provides an extremely precise digital elevation model (DEM) of their whole country. I think visualizing the Matterhorn is a perfect way to jump into 3D data with Matplotlib.<\/p>\n

\"\"<\/a>
Matterhorn. Image credit: Smithsonian Magazine.<\/figcaption><\/figure>\n
\n

1. Download the data.<\/h4>\n

You can access swissALTI3D DEM data here<\/a>.<\/p>\n

They provide the DEM in a couple different formats. I like the idea of working with bare-bones XYZ<\/em> files. With these files, you get a bunch of three-dimensional (x, y, z) coordinates, hopefully in a nice rectangular pattern, written in plain text like a CSV<\/em> file. Here’s a sample of the data:<\/p>\n

X Y Z\r\n2615001 1090999 3495.52\r\n2615003 1090999 3495.18\r\n2615005 1090999 3494.36<\/pre>\n
Z represents elevation above sea level (meters). X and Y are reported in Switzerland’s LV95 coordinate system (also in meters). You could convert the data to a more recognizable latitude-longitude system but it doesn’t seem necessary.<\/p>\n
The DEM divides Switzerland into a grid of 1 km\u00b2 sections. I’ve located the Matterhorn and selected a 3×3 square around the peak. swissALTI3D provides a text file containing URLs for each section’s data. I’ll link that file at the bottom of this post.<\/p>\n
\n
Downloading nine files manually isn’t that difficult but I would still rather use Python. It will make it easy to grab hundreds of files for a future map.<\/p>\n
The project’s folder structure is shown below. We have two subfolders: one for downloaded ZIP<\/em> files and another for XYZ<\/em> files after they’ve been decompressed.<\/p>\n
project_folder\/\r\n\u251c\u2500\u2500 main.py\r\n\u251c\u2500\u2500 zip_files\/\r\n\u2514\u2500\u2500 xyz_files\/<\/pre>\n
Use requests<\/a> and zipfile<\/a> to grab the data. Append file paths to a list so we can reference them later.<\/p>\n
import requests\r\nfrom zipfile import ZipFile\r\nfrom os import listdir\r\nimport pandas as pd\r\n\r\nfile_path_list = []\r\n\r\ndf = pd.read_csv(\"ch.swisstopo.swissalti3d-7Nq9hXI9.csv\", names=[\"url\"])\r\n\r\nfor r, row in df.iterrows():\r\n    response = requests.get(row['url'], stream=True)\r\n    with open(f\"zip_files\/{r}.zip\", \"wb\") as f:\r\n        for chunk in response.iter_content(chunk_size=8192):\r\n            f.write(chunk)\r\n\r\n    with ZipFile(f\"zip_files\/{r}.zip\", 'r') as z:\r\n        z.extractall(f\"xyz_files\/{r}\")\r\n\r\n    file_path_list.append(f\"xyz_files\/{r}\/{listdir(f'xyz_files\/{r}')[0]}\")<\/pre>\n
\n
2. Prepare the data.<\/h4>\n
I’m most comfortable with pandas but you could just as easily read the data with NumPy. It will be converted into NumPy arrays anyway so that approach would make sense.<\/p>\n
Load the nine XYZ<\/em> files into individual pandas DataFrames and append them to a list. Then call pd.concat<\/code> to stack them into one large DataFrame.<\/p>\n
df_list = []\r\n\r\nfor file in file_path_list:\r\n    df_list.append(pd.read_csv(file, delimiter=\" \"))\r\n\r\ndf = pd.concat(df_list)<\/pre>\n
This is optional but I want to convert the Z column (elevation) from meters to feet. Multiply the column by a conversion factor (\u22483.28).<\/p>\n
df['Z'] = df['Z'] * 3.28084<\/pre>\n
df.head()<\/code> looks like this:<\/p>\n
         X        Y             Z\r\n0  2615001  1090999  11468.241837\r\n1  2615003  1090999  11467.126351\r\n2  2615005  1090999  11464.436062\r\n3  2615007  1090999  11454.527926\r\n4  2615009  1090999  11446.817952<\/pre>\n
We’re looking at a small plot of terrain about 3,000 feet below Matterhorn’s peak.<\/p>\n
Notice that X increments by 2. Each 1 km square is further divided into a 500×500 grid of 2-meter squares. Every row of the DataFrame represents one of those squares. Matplotlib will plot the locations in 3D space to visualize the mountain.<\/p>\n
If you do the math, each 1 km square has 500*500=250,000 points. We stacked nine of them together so df<\/code> has 2.25 million rows. That’s easily the largest dataset I’ve worked with on the blog. It’s only ~60 MB so memory shouldn’t be an issue, but it might take your PC a couple minutes to run the script.<\/p>\n
A lot of people say, “Why bother, I can just use Excel.” This is a good example of when you definitely can’t<\/em> do the same thing in Excel.<\/p>\n
\n
Now we have to convert the data into something Matplotlib can understand.<\/p>\n
It wouldn’t make sense to plot all 2.25 million points on a small 1000-pixel image. It would take significantly longer to process and we simply don’t have enough screen resolution to see it.<\/p>\n
Instead we’ll decide on a lower number of divisions\u2014maybe 200\u2014and divide the XY area into a grid of larger squares. Then we’ll use scipy.interpolate.griddata<\/code> to calculate Z values for the new grid map.<\/p>\n
np.linspace<\/code> creates an evenly spaced array ranging from minimum to maximum values. np.meshgrid<\/code> is something I’ve written about in greater detail<\/a> before. It essentially reshapes a pair of 1D arrays into a pair of 2D arrays with corresponding points. You can think of them like a map of all the grid points.<\/p>\n
To visualize terrain, we combine them with another 2D array that holds elevation data, which is where griddata<\/code> comes in. This method accepts the original (x, y, z) data and the new grid arrays.<\/p>\n
import numpy as np\r\nfrom scipy.interpolate import griddata\r\n\r\nx = df['X']\r\ny = df['Y']\r\nz = df['Z']\r\n\r\nx_grid = np.linspace(x.min(), x.max(), 200)\r\ny_grid = np.linspace(y.min(), y.max(), 200)\r\nx_grid, y_grid = np.meshgrid(x_grid, y_grid)\r\n\r\nz_elevation = griddata((x, y), z, (x_grid, y_grid), method=\"cubic\")<\/pre>\n
You’ll notice diminishing returns after 200 divisions. Even this plot won’t look significantly better than 100 but we can err on the high side.<\/p>\n
Here’s a diagram that shows a 1 km\u00b2 chunk of Switzerland at various resolutions. You can spot finer features with higher resolution, but obviously there’s a tradeoff with CPU time.<\/p>\n
\"\"<\/a><\/p>\n
You might think, “If I set the argument in linspace() to 500, that would be like using all the data in the XYZ<\/em> file. Why wouldn’t I just generate the Z array myself?”<\/p>\n
You could do that. But you’d probably<\/em> end up writing a nested for<\/code> loop that’s significantly slower than scipy.interpolate.griddata<\/code>. If you’re a wizard then feel free. I know when I’m outmatched so I’ll happily rely on SciPy.<\/p>\n
\n
3. Plot the data.<\/h4>\n
I’ll stick with default Matplotlib and skip creating a custom mplstyle. I want to write the simplest possible example (within reason) to help someone get started with 3D data.<\/p>\n
Instantiating a 3D Axes is a little more awkward than 2D. It’s best to first create a Figure and then<\/em> create an Axes. “111” is a quick way of saying, “1 row, 1 column, 1st position.” Use subplots_adjust<\/code> to specify the plot’s outer margins. Units are relative to the entire figure (0.0 to 1.0).<\/p>\n
import matplotlib.pyplot as plt\r\n\r\nfig = plt.figure(figsize=(10, 8.4))\r\nax = fig.add_subplot(111, projection=\"3d\")\r\nfig.subplots_adjust(left=0.04, right=0.99, bottom=0.01, top=0.99)<\/pre>\n
The important business happens in plot_surface<\/code>. Here we pass the three 2D NumPy arrays we just created.<\/p>\n
cmap<\/code> refers to color map, which will automatically color the surface according to Z value. In this application Z represents elevation, so we’ll see a gradient from the base of Matterhorn to its peak. linewidth<\/code> and edgecolor<\/code> give the surface a sort of wire frame that makes it easier to see when rendered on screen. You could omit them but I think they improve the visualization.<\/p>\n
surface = ax.plot_surface(x_grid, y_grid, z_elevation, cmap=\"terrain\", linewidth=0.5, edgecolor=\"#444\")<\/pre>\n
Let’s also create a colorbar<\/code>, which is like a legend for Z values. It tells us exactly what elevation is represented by a given color. By default colorbar<\/code> steals part of its parent’s Axes, so it will slide in to the surface’s right.<\/p>\n
shrink<\/code> is the color bar’s size relative to its Axes. 0.5 means it will be half the overall Axes height. aspect<\/code> defines height-width ratio, so this color bar will be tall and skinny. There’s no way to set a color bar title directly but we can call a method on the instance after it’s created.<\/p>\n
color_bar = fig.colorbar(surface, shrink=0.5, aspect=10)\r\ncolor_bar.set_label(\"Elevation (ft)\")<\/pre>\n
Set a title for the plot as well.<\/p>\n
ax.set_title(\"Matterhorn Terrain Map\", size=11)<\/pre>\n
view_init<\/code> positions the camera. It wouldn’t be necessary if we were calling plt.show()<\/code>. Once the figure window appeared you could click and drag the camera wherever you wanted. But it’s very important to locate the camera when producing a static PNG image. We only get one look at the surface.<\/p>\n
elev<\/code> (elevation) moves the camera up and down. azim<\/code> (azimuth) rotates it like the moon rotates around the Earth. The chosen values get close to matching the Matterhorn picture on its Wikipedia page (shown below).<\/p>\n
ax.view_init(elev=10, azim=45)<\/pre>\n
Finally, save the figure. It might take a minute or two for the script to run, depending on your internet and CPU speeds.<\/p>\n
plt.savefig(\"matterhorn.png\", dpi=200)<\/pre>\n
\n
4. The output.<\/h4>\n
This shows a 3 km x 3 km area surrounding the Matterhorn. It isn’t the highest peak in the world but it’s arguably the most striking. You can see why it’s so famous.<\/p>\n
\"\"<\/a><\/p>\n
Real life for comparison:<\/p>\n
\"\"<\/a>
Matterhorn. Image credit: Wikipedia.<\/figcaption><\/figure>\n
There’s a lot of empty white space in the output image. It’s there because when you do plt.show<\/code> and move the camera around, axis labels need a lot of space. There’s only so much you can do about it natively in Matplotlib. If it bothers you (like it does me), one solution is to use OpenCV<\/a> to crop out white space.<\/p>\n
That’s what I did to create this GIF:<\/p>\n
\"\"<\/a><\/p>\n
I love how the animation turned out. Mountains are so big that it’s tough to get a good look around them. The code is essentially the same as above except it generates 180 images while incrementing the azimuth two degrees each frame.<\/p>\n
I like this script because it takes the most basic possible three-dimensional spatial data\u2014(x, y, z) coordinates in a CSV<\/em> file\u2014and outputs a nice visualization. You don’t need a proprietary database format and a whole software package to work with it. Even in 3D, you can get the job done with open-source software.\n<\/span>\n\n<\/p>\n
\n
Download the data.<\/strong><\/a><\/p>\n
Full code:<\/strong><\/p>\n
import requests\r\nfrom zipfile import ZipFile\r\nfrom os import listdir\r\nimport pandas as pd\r\nimport numpy as np\r\nfrom scipy.interpolate import griddata\r\nimport matplotlib.pyplot as plt\r\n\r\n\r\nfile_path_list = []\r\n\r\ndf = pd.read_csv(\"ch.swisstopo.swissalti3d-7Nq9hXI9.csv\", names=[\"url\"])\r\n\r\nfor r, row in df.iterrows():\r\n    response = requests.get(row['url'], stream=True)\r\n    with open(f\"zip_files\/{r}.zip\", \"wb\") as f:\r\n        for chunk in response.iter_content(chunk_size=8192):\r\n            f.write(chunk)\r\n    \r\n    with ZipFile(f\"zip_files\/{r}.zip\", 'r') as z:\r\n        z.extractall(f\"xyz_files\/{r}\")\r\n\r\n    file_path_list.append(f\"xyz_files\/{r}\/{listdir(f'xyz_files\/{r}')[0]}\")\r\n\r\ndf_list = []\r\n\r\nfor file in file_path_list:\r\n    df_list.append(pd.read_csv(file, delimiter=\" \"))\r\n\r\ndf = pd.concat(df_list)\r\n\r\ndf['Z'] = df['Z'] * 3.28084\r\n\r\nx = df['X']\r\ny = df['Y']\r\nz = df['Z']\r\n\r\nx_grid = np.linspace(x.min(), x.max(), 200)\r\ny_grid = np.linspace(y.min(), y.max(), 200)\r\nx_grid, y_grid = np.meshgrid(x_grid, y_grid)\r\n\r\nz_elevation = griddata((x, y), z, (x_grid, y_grid), method=\"cubic\")\r\n\r\nfig = plt.figure(figsize=(10, 8.4))\r\nax = fig.add_subplot(111, projection=\"3d\")\r\nfig.subplots_adjust(left=0.04, right=0.99, bottom=0.01, top=0.99)\r\n\r\nsurface = ax.plot_surface(x_grid, y_grid, z_elevation, cmap=\"terrain\", linewidth=0.5, edgecolor=\"#444\")\r\n\r\ncolor_bar = fig.colorbar(surface, shrink=0.5, aspect=10)\r\ncolor_bar.set_label(\"Elevation (ft)\")\r\n\r\nax.set_title(\"Matterhorn Terrain Map\", size=11)\r\n\r\nax.view_init(elev=10, azim=45)\r\n\r\nplt.savefig(\"matterhorn.png\", dpi=200)<\/pre>\n
 <\/p>\n","protected":false},"excerpt":{"rendered":"
I’ve wanted to post a Matplotlib 3D visualization for a while, but sharing it on the web is difficult. It’s easy to save 2D plots and embed them on this page. For a 3D plot, it works better when you<\/p>\n","protected":false},"author":1,"featured_media":2552,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[19,561,279],"tags":[507,506,387,517,520,39,515,519,526,278,22,122,518,530,529,438,522,524,72,521,73,271,24,509,388,510,75,514,30,512,505,25,523,513,527,508,511,528,516,141,525],"class_list":["post-2523","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-maps","category-animated","category-nature","tag-2d","tag-3d","tag-animated","tag-azimuth","tag-cmap","tag-code","tag-colorbar","tag-colormap","tag-concat","tag-coordinates","tag-data","tag-dataset","tag-elevation","tag-geography","tag-geology","tag-gif","tag-griddata","tag-interpolate","tag-latitude","tag-linspace","tag-longitude","tag-map","tag-matplotlib","tag-matterhorn","tag-meshgrid","tag-mountain","tag-numpy","tag-opencv","tag-pandas","tag-plot_surface","tag-programming","tag-python","tag-scipy-interpolate-griddata","tag-surface","tag-swissalti3d","tag-switzerland","tag-terrain-map","tag-topography","tag-view_init","tag-visualization","tag-xyz"],"_links":{"self":[{"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/posts\/2523","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/comments?post=2523"}],"version-history":[{"count":34,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/posts\/2523\/revisions"}],"predecessor-version":[{"id":3361,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/posts\/2523\/revisions\/3361"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/media\/2552"}],"wp:attachment":[{"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/media?parent=2523"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/categories?post=2523"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/wollen.org\/blog\/wp-json\/wp\/v2\/tags?post=2523"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}