# Import libraries
import sys
import pygmt

# Show versions
print("Python version:", sys.version)
print("PyGMT version:", pygmt.__version__)

Python version: 3.11.13 | packaged by conda-forge | (main, Jun  4 2025, 14:39:58) [MSC v.1943 64 bit (AMD64)]
PyGMT version: v0.17.0

# Download elevation grid
grid = pygmt.datasets.load_earth_relief(resolution="15m", region="g")

# Create custom colormap (saved as a .cpt file in the directory the script is run in)
cpt_text = """# desert.cpt
-8000     210/180/140    8000    210/180/140
B         210/180/140
F         210/180/140
N         210/180/140
"""

# Write the above text to a .cpt file
# Using this "with" format suppresses a print of the number of characters to the terminal
with open("desert.cpt", "w") as f:
    f.write(cpt_text)

# Do the same for a new colormap
cpt_text = """# desert_with_ocean.cpt
-8000     0/0/139    -4000    173/216/230
B         210/180/140
F         210/180/140
N         210/180/140
"""
with open("desert_with_ocean.cpt", "w") as f:
    f.write(cpt_text)

# FIRST FIGURE

# Initialize figure
fig = pygmt.Figure()

# Create hillshade effect
shade = pygmt.grdgradient(grid=grid, azimuth="0/0", normalize="t1")

# Create basemap
fig.basemap(
    region="d", # "d" for whole globe
    projection="N12c", # Robinson projection
    frame=["WSne+tDry Oceans", "xa40", "ya40"], # Captial W and S mean that number lablels are added to those axes
)

# Add elevation grid
fig.grdimage(
    grid=grid, # Use earth relief grid
    shading=shade, 
    cmap="desert.cpt",
)

fig.shift_origin(xshift="14c")

# SECOND FIGURE

#Same basemap as firt figure
fig.basemap(
    region="d",
    projection="N12c",
    frame=["WSne+tOcean Surface at -4000 m", "xa40", "ya40"],
)

# Add elevation grid
fig.grdimage(
    grid=grid,
    shading=shade,
    cmap="desert_with_ocean.cpt",
)

#Add a colorbar to this figure for the ocean depth
fig.colorbar(frame=["x+lElevation (m)"])

#Show and save figure
#fig.show()
#fig.savefig("city_and_the_stars.jpg")