-
Notifications
You must be signed in to change notification settings - Fork 214
/
grdview_surface.py
54 lines (45 loc) · 1.59 KB
/
grdview_surface.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
"""
Plotting a surface
------------------
The :meth:`pygmt.Figure.grdview()` method can plot 3-D surfaces with
``surftype="s"``. Here, we supply the data as an :class:`xarray.DataArray` with
the coordinate vectors ``x`` and ``y`` defined. Note that the ``perspective``
parameter here controls the azimuth and elevation angle of the view. We provide
a list of two arguments to ``frame`` - the first argument specifies the
:math:`x`- and :math:`y`-axes frame attributes and the second argument,
prepended with ``"z"``, specifies the :math:`z`-axis frame attributes.
Specifying the same scale for the ``projection`` and ``zscale`` parameters
ensures equal axis scaling. The ``shading`` parameter specifies illumination;
here we choose an azimuth of 45° with ``shading="+a45"``.
"""
import numpy as np
import pygmt
import xarray as xr
# Define an interesting function of two variables, see:
# https://en.wikipedia.org/wiki/Ackley_function
def ackley(x, y):
return (
-20 * np.exp(-0.2 * np.sqrt(0.5 * (x**2 + y**2)))
- np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y)))
+ np.exp(1)
+ 20
)
# Create gridded data
INC = 0.05
x = np.arange(-5, 5 + INC, INC)
y = np.arange(-5, 5 + INC, INC)
data = xr.DataArray(ackley(*np.meshgrid(x, y)), coords=(x, y))
fig = pygmt.Figure()
# Plot grid as a 3-D surface
SCALE = 0.5 # in centimeters
fig.grdview(
data,
frame=["a5f1", "za5f1"],
projection=f"x{SCALE}c",
zscale=f"{SCALE}c",
surftype="s",
cmap="roma",
perspective=[135, 30], # Azimuth southeast (135°), at elevation 30°
shading="+a45",
)
fig.show()