Introduction

Using Plotly, which supports various languages including Python and allows you to draw interactive diagrams relatively easily, you can use elevation data to rotate around like Google Earth. Let's create a globe.

The result looks like this. https://rkiuchir.github.io/3DSphericalTopo

3 points

Acquisition and reading of elevation data
Convert latitude / longitude information expressed in Cartesian coordinate system to spherical coordinate system ([Plotly Chart Studio: Heatmap plot on a spherical map](https://chart-studio.plotly.com/~empet/14813/heatmap] -see plot-on-a-spherical-map / # /))
Draw with Plotly

Execution environment

macOS 10.14 Mojave
Python 3.6 (Anaconda)
Jupyter notebook

0. Install Plotly

It can be installed with pip install plotly.

1-1. Acquisition of altitude data

Elevation data (ETOPO1 Global Relief Model) with a resolution of 1 minute globally is available in NOAA (National Oceanic and Atmospheric Administration). This site](https://www.ngdc.noaa.gov/mgg/global/) is available.
This time, download the grid-registered netCDF format file (ETOPO1_Ice_g_gdal.grd) of ETOPO1 Ice Surface.

1-2. Reading elevation data

First, read the three data of latitude / longitude and altitude from the data file in netCDF format.

Input defined here is the latitude (lat_area), longitude (lon_area), and resolution (resolution) of the specified area.
Output is the latitude (lon), longitude (lat), and resolution (topo) of the mesh type in the specified area.
The resolution part is specified in the format of ° (degree), and the highest resolution specification is 0.0167 °.

Regarding the resolution, the data is skipped according to the specified resolution value.

import numpy as np
from netCDF4 import Dataset

def Etopo(lon_area, lat_area, resolution):
    ### Input
    # resolution: resolution of topography for both of longitude and latitude [deg]
    # (Original resolution is 0.0167 deg)
    # lon_area and lat_area: the region of the map which you want like [100, 130], [20, 25]
    ###
    

    ### Output
    # Mesh type longitude, latitude, and topography data
    ### 

    # Read NetCDF data
    data = Dataset("ETOPO1_Ice_g_gdal.grd", "r")
    
    # Get data
    lon_range = data.variables['x_range'][:]
    lat_range = data.variables['y_range'][:]
    topo_range = data.variables['z_range'][:]
    spacing = data.variables['spacing'][:]
    dimension = data.variables['dimension'][:]
    z = data.variables['z'][:]

    lon_num = dimension[0]
    lat_num = dimension[1]
    
    # Prepare array
    lon_input = np.zeros(lon_num); lat_input = np.zeros(lat_num)
    for i in range(lon_num):
        lon_input[i] = lon_range[0] + i * spacing[0]

    for i in range(lat_num):
        lat_input[i] = lat_range[0] + i * spacing[1]

    # Create 2D array
    lon, lat = np.meshgrid(lon_input, lat_input)

    # Convert 2D array from 1D array for z value
    topo = np.reshape(z, (lat_num, lon_num))

    # Skip the data for resolution
    if ((resolution < spacing[0]) | (resolution < spacing[1])):
        print('Set the highest resolution')
    else:
        skip = int(resolution/spacing[0])
        lon = lon[::skip,::skip]
        lat = lat[::skip,::skip]
        topo = topo[::skip,::skip]

    topo = topo[::-1]

    # Select the range of map
    range1 = np.where((lon>=lon_area[0]) & (lon<=lon_area[1]))
    lon = lon[range1]; lat = lat[range1]; topo = topo[range1]
    range2 = np.where((lat>=lat_area[0]) & (lat<=lat_area[1]))
    lon = lon[range2]; lat = lat[range2]; topo = topo[range2]
    
    # Convert 2D again
    lon_num = len(np.unique(lon))
    lat_num = len(np.unique(lat))
    lon = np.reshape(lon, (lat_num, lon_num))
    lat = np.reshape(lat, (lat_num, lon_num))
    topo = np.reshape(topo, (lat_num, lon_num))
   
    return lon, lat, topo

2. Convert latitude / longitude information to spherical coordinates

(See Plotly Chart Studio: Heatmap plot on a spherical map)

Here, the latitude / longitude information expressed in the orthogonal system coordinates prepared above is converted to the spherical coordinate system.

def degree2radians(degree):
    # convert degrees to radians
    return degree*np.pi/180

def mapping_map_to_sphere(lon, lat, radius=1):
    #this function maps the points of coords (lon, lat) to points onto the  
    sphere of radius radius
    
    lon=np.array(lon, dtype=np.float64)
    lat=np.array(lat, dtype=np.float64)
    lon=degree2radians(lon)
    lat=degree2radians(lat)
    xs=radius*np.cos(lon)*np.cos(lat)
    ys=radius*np.sin(lon)*np.cos(lat)
    zs=radius*np.sin(lat)
    return xs, ys, zs

3. Draw with Plotly

Now, let's actually draw the 3D data of latitude, longitude, and altitude represented by the spherical coordinate system in Plotly.

First, call the function prepared in 1-2. To read the global elevation data. If you read at too high a resolution, the amount of data will increase on the order of the cube, so this time the resolution is set to 0.8 °.

# Import topography data
# Select the area you want
resolution = 0.8
lon_area = [-180., 180.]
lat_area = [-90., 90.]

# Get mesh-shape topography data
lon_topo, lat_topo, topo = ReadGeo.Etopo(lon_area, lat_area, resolution)

Next, convert it to the spherical coordinate system with the function prepared in 2.

xs, ys, zs = mapping_map_to_sphere(lon_topo, lat_topo)

And from here, we will actually move on to drawing.

First, define the color scale used to draw the elevation data.

# Import color scale
import Plotly_code as Pcode
name = "topo"
Ctopo = Pcode.Colorscale_Plotly(name)
cmin = -8000
cmax = 8000

Then, draw using Plotly. Here you enter the input data and color scale.

topo_sphere=dict(type='surface',
            x=xs, 
            y=ys, 
            z=zs,
            colorscale=Ctopo,
            surfacecolor=topo,
            cmin=cmin, 
            cmax=cmax)
            )

Erase the shaft etc. so that it looks good.

noaxis=dict(showbackground=False,
            showgrid=False,
            showline=False,
            showticklabels=False,
            ticks='',
            title='',
            zeroline=False)

Finally, use layout to specify the title and background color. This time, the background color is black, with a little awareness of Google Earth.


import plotly.graph_objs as go

titlecolor = 'white'
bgcolor = 'black'

layout = go.Layout(
      autosize=False, width=1200, height=800,
      title = '3D spherical topography map',
      titlefont = dict(family='Courier New', color=titlecolor),
      showlegend = False,
      scene = dict(
        xaxis = noaxis,
        yaxis = noaxis,
        zaxis = noaxis,
        aspectmode='manual',
        aspectratio=go.layout.scene.Aspectratio(
            x=1, y=1, z=1)),
        paper_bgcolor = bgcolor,
        plot_bgcolor = bgcolor)

Then, draw using the prepared one (html output here).

from plotly.offline import plot

plot_data=[topo_sphere]

fig = go.Figure(data=plot_data, layout=layout)
plot(fig, validate = False, filename='3DSphericalTopography.html', 
     auto_open=True)

in conclusion

With this, I think that I could draw a plot of a globe that can be turned like the beginning. I use this by overlaying the distribution of earthquakes on top of this.

This example can be easily applied when plotting on a sphere, and I think it has a wide range of uses. In addition, it is useful not only for elevation data but also for creating a two-dimensional map, and it is also possible to draw three-dimensional images that express the height with elevation data.

[Python] Draw elevation data on a sphere with Plotly and draw a globe that can be rotated round and round