Hello. Could someone please help me with the following issue? I am trying to visualize a very detailed DMR (Digital Terrain Model) for the Slovak Republic, which I have as a TIFF file on my hard drive. I would like to set the plane color to grey, but the color fill does not appear. I’m not sure what I’m doing wrong. When I used a similar script with pygmt.datasets.load_earth_relief as the data source, everything worked correctly. I suspect the problem lies in my raster file.
> Blockquote
import os
os.environ["PROJ_LIB"] = r"C:\Users\maria\anaconda3\pkgs\proj-9.6.2-h4f671f6_0\Library\share\proj"
import pygmt
import rasterio
import numpy as np
# Set output directory
output_dir = r"C:\Users\maria\Documents\data_analysis\Anaconda_scripts\pygmt_dmr_zobrazenie\results"
os.makedirs(output_dir, exist_ok=True)
# Path to your DEM (GeoTIFF)
dem_path = r"C:\Users\maria\Documents\data_analysis\Anaconda_scripts\pygmt_dmr_zobrazenie\dmr_raster_squares\pdbat_3701_1_A_square.tif"
# Get bounds and min/max values
with rasterio.open(dem_path) as src:
bounds = src.bounds
dem_data = src.read(1, masked=True)
min_elev = float(np.nanmin(dem_data))
max_elev = float(np.nanmax(dem_data))
print(f"DEM has a minimum elevation of {min_elev:.1f} and maximum of {max_elev:.1f} meters relative to sea level.")
# Set color scale before plotting
pygmt.makecpt(cmap="geo", series=[min_elev, max_elev, 1], continuous=True)
# Set projection and region (adjust region as needed, here it is based on bounds)
region = [bounds.left, bounds.right, bounds.bottom, bounds.top, min_elev, max_elev]
# Calculate aspect ratio for projection
width = bounds.right - bounds.left
height = bounds.top - bounds.bottom
aspect_ratio = width / height
proj_width = 15 # cm
proj_height = proj_width / aspect_ratio
projection = f"X{proj_width}c/{proj_height:.2f}c"
fig = pygmt.Figure()
pygmt.config(FORMAT_GEO_MAP="dddF", MAP_FRAME_TYPE="fancy")
frame = ["xaf", "yaf", "zaf", "+b"]
with pygmt.config(FONT="22p"):
fig.grdview(
grid=dem_path,
region=region,
projection=projection,
frame=frame,
perspective=[120, 25],
zsize="1.5c",
surftype="i",
cmap="geo", # <-- explicitly set cmap
plane=f"{60 + min_elev}+ggrey",
shading="+a20+nt1"
)
fig.basemap(
perspective=True,
rose="JCL+w5c+l+o1c/0c"
)
with pygmt.config(FONT="16p"):
fig.colorbar(frame=["a1000", "x+lElevation", "y+lm"], position="+w15c")
output_path = os.path.join(output_dir, "dem_3d_antarctica.jpg")
fig.savefig(output_path, dpi=300)
print(f"3D image has been saved to {output_path}")
# Automatically open the image in the default viewer
os.startfile(output_path)
> Blockquote
