This script takes a .tif layer and produces an elevation profile chart that always includes the highest point of the layer:

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

def create_elevation_profile(tif_file):
    dataset = gdal.Open(tif_file)
    band = dataset.GetRasterBand(1)
    elevation_data = band.ReadAsArray()

    # Find the coordinates of the highest point
    highest_point_indices = np.unravel_index(np.argmax(elevation_data), elevation_data.shape)
    highest_point_elevation = elevation_data[highest_point_indices]

    # Calculate the distance values for the elevation profile
    distance = np.arange(elevation_data.shape[1])

    # Extract the elevation profile along the highest point
    elevation_profile = elevation_data[highest_point_indices[0], :]

    # Plot the elevation profile with the highest point included
    plt.plot(distance, elevation_profile, color="black")
    plt.xlabel("Distance (pixels)")
    plt.ylabel("Elevation (m)")
    plt.title("Elevation Profile")
    plt.grid(True)

    # Add a marker for the highest point
    plt.scatter(highest_point_indices[1],
                highest_point_elevation,
                color="red",
                label=f"Highest Point: {round(highest_point_elevation, 1)}m")
    plt.legend()
    plt.show()

# Provide the path to your TIF file
tif_file_path = "where/is/your/.tif"

create_elevation_profile(tif_file_path)

Enjoy!

If you have any suggestions or find any issues with the code then get in contact at thespatialspace@gmail.com

Leave a comment

Latest Stories