# Import libraries
import pygmt
import sys

# 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

# Define region of map and create grid
region = [14.70, 14.99, -23.93, -23.68]  # Dunes in Namibia's coastal region
grid = pygmt.datasets.load_earth_relief(resolution='01s', region=region) #Load earth elevation Xarray DataArray for the region

#---MAIN MAP---

# Initialize figure
fig = pygmt.Figure()

# Set coordinates for elevation profile line
lonA, lonB = 14.70, 14.99
latA = latB = -23.78

# Create subtle hillshade effect
hillshade = pygmt.grdgradient(grid=grid, azimuth="0/0", normalize="t0.2")

#Plot elevation grid
fig.grdimage(
    grid=grid,
    shading=hillshade,
    region=region,
    projection="M12c",
    cmap="copper",
    frame="af"
)

# Plot the survey line over the grid
fig.plot(
    x=[lonA, lonB], 
    y=[latA, latB], 
    pen="1p,red,solid"
)

# Draw circles at the ends of the line
fig.plot(
    x=[lonA, lonB],
    y=[latA, latB],
    style="c0.3c",   # circle of 0.3 cm
    fill="red",
    pen="black"      # optional black outline
)

# Label "A" for western end
fig.text(
    x=lonA,
    y=latA,
    text="A",
    offset="0.2c/0.4c",
    font="15p,red"
)

# Label "B" for eastern end
fig.text(
    x=lonB,
    y=latB,
    text="B",
    offset="-0.2c/0.4c",
    font="15p,red"
)

# Add a map scale
fig.basemap(
    map_scale="jBR+w5k+o0.5c/0.75c+f+l"
)

# Add a colorbar
fig.colorbar(frame=["x+lElevation (m)"], position="JMR+o1c/0c+w8c")

# ---INSET MAP---

# Define Namibia region for inset map
namibia_region = [11.5, 25.5, -29.0, -16.5]  # lon_min, lon_max, lat_min, lat_max

# Add inset map
with fig.inset(position="jBL+o0.5c+w4c/3.8c", box="+p4p,gold"):
    
    # Plot Namibia coastline
    fig.coast(
        region=namibia_region,
        projection="M4c", 
        land="lightgray",
        water="lightblue",
        borders="1",
        shorelines="0.5p,blue",
        frame=False
    )
    
    # Define and plot map region rectangle within inset map 
    rectangle = [[region[0], region[2], region[1], region[3]]]
    fig.plot(
        data=rectangle, 
        style="r+s", 
        pen="2p,red"
    )

    # Add "Namibia" text
    fig.text(
        text="Namibia",
        region=namibia_region,
        projection="M4c",
        x=17.5,  # approximate center longitude
        y=-20.5,  # approximate center latitude
        font="8p,Helvetica-Bold,black",
        justify="CM"  # Center-Middle alignment
    )
    
#---ELEVATION PROFILE---
fig.shift_origin(yshift="12c")

# Define region
fig.basemap(
    region=[lonA, lonB, 100, 600],  # x_min, x_max, y_min, y_max 
    projection="X12c/3c", # Cartesian projection with a width of 12 centimeters and a height of 3 centimeters
    frame=0
)

# Draw line on map associated with the elevation data underlying it
profile_df = pygmt.project(
    center=[lonA, latA],  # Start point of survey line (longitude, latitude)
    endpoint=[lonB, latB],  # End point of survey line (longitude, latitude)
    generate=0.0001,  # Output data in steps of 0.1 degrees
)

# Create a new data column 'elevation' to go with the line
profile_df = pygmt.grdtrack(grid=grid, points=profile_df, newcolname="elevation")

# Plot elevation along the survey line
fig.plot(
    x=profile_df.r, #x value in degrees
    y=profile_df.elevation,
    fill="tan",
    pen="1p,darkbrown,solid",
    close="+y100",  # Force closed polygon
)

# Add 'A' and 'B' to ends of the line
fig.text(
    x=[lonA, lonB],
    y=[650, 650],
    text=["A", "B"],
    no_clip=True,  # Do not clip text that fall outside the plot bounds
    font="10p,red",
)

# Draw black frame and title around the plot
fig.basemap(frame=["WSrt+tElevation Profile of Namibia Dunes", "ya100+lElevation (m)"])

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