# Import Libraries
import pygmt
import pandas as pd
import numpy as np
import sys
import scipy
from scipy.ndimage import generic_filter
import xarray as xr

# Show versions
print("Python version:", sys.version)
print("PyGMT version:", pygmt.__version__)
print("Pandas version:", pd.__version__)
print("NumPy version:", np.__version__)
print("SciPy version:", np.__version__)
print("Xarray version:", xr.__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
Pandas version: 2.3.3
NumPy version: 2.3.3
SciPy version: 2.3.3
Xarray version: 2025.10.1

# Define region of map and create grid
region=[-65, -30, 45, 60] #Approximate coordinates of Sinus Roris region
moon_grid = pygmt.datasets.load_moon_relief(resolution='02m', region=region, registration=None) #Load moon elevation data

# FIRST FIGURE: Topographic map
fig = pygmt.Figure()

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

# Plot grid
fig.grdimage(
    grid=moon_grid,
    shading=shade,
    region=region,
    projection="L-47.5/52.5/45/60/12c",
    frame=["a", "+tSinus Roris Region (Moon)"],
    cmap="rainbow",
)

# Add star in proposed location of Port Roris
fig.plot(
    x=-48.604, 
    y=48.329, 
    style="a0.3c",  # star of radius 0.2 cm
    fill="red", 
    pen="black"     # black outline for visibility
)

# Add text label next to the dot
fig.text(
    x=-48.604 - 2.25,   # slightly to the right of Port Roris
    y=48.329 + 0.05,  # slightly above Port Roris
    text="Port Roris",
    font="6p,Helvetica-Bold,black"
)

#Add map scale
fig.basemap(
    map_scale="jBL+w500k+o0.5c/1c+f+l"
)

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

fig.shift_origin(xshift="16c")

# SECOND FIGURE: Local Flatness
moon_grid_array = moon_grid.values  # Convert xarray DataArray to NumPy array

# Define function to calculate standard deviation within a surrounding window of pixels
def local_std(window):
    return np.std(window)

# #Define window size and perform standard deviation calcuation
window_size = 10 
std_grid = generic_filter(moon_grid_array, local_std, size=window_size)

#Convery NumPy array back to an Xarray DataArray
std_da = xr.DataArray(
    std_grid,
    coords={"lat": moon_grid.lat, "lon": moon_grid.lon},
    dims=["lat", "lon"]
)

# Plot grid
fig.grdimage(
    grid=std_da,
    region=region,
    projection="L-47.5/52.5/45/60/12c",
    frame=["a", "+tLocal Flatness"],
    cmap="rainbow",
)

# Add colorbar
fig.colorbar(frame=[f"x+lLocal @~s@~ of elevation ({window_size} px window)"], position="JMR+o1c/0c+w8c")

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