# %%
import pygmt
from pygmt.params import Position

# Define a study area with huge elevation changes
region = [141, 147, 36, 43]

# Download a grid for the Earth relief with a resolution of 10 arc-minutes
grid = pygmt.datasets.load_earth_relief(resolution="10m", region=region)

# Convert the grid into a pandas DataFrame
grd_df = pygmt.grd2xyz(grid=grid)
zmin = grd_df["z"].min() - 50
zmax = grd_df["z"].max() + 50

# Create a 3-D bar plot with color-coding
fig = pygmt.Figure()

fig.basemap(
    region=[*region, zmin, zmax],
    projection="M10c",
    zsize="8c",
    frame=["WSneZ", "xaf", "yag", "za1000f500+lElevation / m"],
    perspective=(195, 30),
)

pygmt.makecpt(cmap="SCM/oleron", series=(zmin, zmax))
fig.plot3d(
    data=grd_df,
    # Use "o" to plot bars and give the desired size
    # The base of the bars is set via "+b"
    style=f"o0.34c+b{zmin}",
    cmap=True,
    pen="0.01p,gray30",
    perspective=True,
    # Use the elevation for the color-coding of the bars
    # Zero-based indexing for columns
    incols=[0, 1, 2, 2],  # lon, lat, z for elevation, z for color-coding
)
fig.colorbar(
    frame=["xa1000f500+lElevation", "y+lm"],
    position=Position("TR", cstype="inside", offset=1.4),
    orientation="vertical",
    length=7,
    move_text="label",
    label_as_column=True,
)

fig.show()