Terrain Analysis and Geolocation

Terrain shapes everything in military and intelligence operations: what you can see, where you can go, where you can hide. Digital Elevation Models (DEMs) let you analyze terrain computationally — computing viewsheds, slope, aspect, and line-of-sight from satellite data. Combined with shadow analysis, terrain features also let you geolocate imagery and measure object heights.

Connection to OSINT: Geolocation from imagery is a core OSINT skill. This note covers the computational/satellite side — matching terrain features to DEMs, using shadows for height measurement. Cross-reference your Security course’s Geolocation Techniques.

---

Digital Elevation Models (DEM)

A DEM is a raster where each pixel stores elevation above sea level (or above a reference ellipsoid).

Free Global DEMs

DEM	Resolution	Coverage	Source	Notes
SRTM	30m	56S-60N	NASA (2000 Shuttle mission)	Widely used, some voids in steep terrain
Copernicus DEM	30m (GLO-30)	Global	ESA/Airbus	Better quality than SRTM, newer
ASTER GDEM	30m	83S-83N	NASA/Japan	Noisier than SRTM
ALOS World 3D	30m	Global	JAXA	Good in areas where SRTM has voids

Note: Estonia at 59N is near the northern limit of SRTM. Copernicus DEM has better coverage at high latitudes.

DSM vs DTM vs DEM

• DEM (Digital Elevation Model): generic term for elevation data
• DSM (Digital Surface Model): includes buildings, trees — top of everything
• DTM (Digital Terrain Model): bare earth only (buildings/trees removed)
• SRTM and Copernicus GLO-30 are approximately DSMs (they see the top of canopy/buildings)

Loading and Displaying a DEM

import rasterio
import numpy as np
import matplotlib.pyplot as plt
 
def load_dem(dem_path):
    """Load a DEM and return elevation array + metadata."""
    with rasterio.open(dem_path) as src:
        elevation = src.read(1).astype(np.float32)
        transform = src.transform
        crs = src.crs
        nodata = src.nodata
 
    if nodata is not None:
        elevation[elevation == nodata] = np.nan
 
    return elevation, transform, crs
 
 
def display_dem(elevation, title="Digital Elevation Model", cmap="terrain"):
    """Display a DEM with contour lines."""
    fig, ax = plt.subplots(figsize=(10, 8))
 
    im = ax.imshow(elevation, cmap=cmap)
    plt.colorbar(im, ax=ax, label="Elevation (m)")
 
    # Add contour lines
    valid = np.nan_to_num(elevation, nan=0)
    contours = ax.contour(valid, levels=15, colors="black", linewidths=0.3, alpha=0.5)
    ax.clabel(contours, inline=True, fontsize=6, fmt="%.0f")
 
    ax.set_title(title)
    ax.axis("off")
    plt.tight_layout()
    return fig
 
# Usage:
# elevation, transform, crs = load_dem("srtm_n59_e024.tif")
# display_dem(elevation, "SRTM — Tallinn Area")
# plt.show()

---

Derived Products

Slope

Steepness of terrain at each pixel. Critical for mobility analysis (vehicles can’t climb steep slopes) and erosion assessment.

import numpy as np
 
def compute_slope(elevation, pixel_size_m=30):
    """
    Compute slope in degrees from a DEM.
    Uses numpy gradient for partial derivatives.
    """
    # Partial derivatives (dz/dx, dz/dy) in meters
    dy, dx = np.gradient(elevation, pixel_size_m)
 
    # Slope in radians, then degrees
    slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
    slope_deg = np.degrees(slope_rad)
 
    return slope_deg
 
# Mobility classification:
# 0-5 deg: easy for all vehicles
# 5-15 deg: moderate, wheeled vehicles struggle in wet conditions
# 15-30 deg: difficult, tracked vehicles only
# 30-45 deg: very difficult, dismounted infantry only
# >45 deg: impassable

Aspect

Direction the terrain faces (compass direction of the steepest downhill slope).

def compute_aspect(elevation, pixel_size_m=30):
    """
    Compute aspect (direction terrain faces) in degrees.
    0/360 = North, 90 = East, 180 = South, 270 = West.
    """
    dy, dx = np.gradient(elevation, pixel_size_m)
    aspect = np.degrees(np.arctan2(-dx, dy))  # negative dx because x increases east
    aspect = (aspect + 360) % 360  # convert to 0-360
    return aspect
 
# Intelligence note: south-facing slopes get more sun (Northern Hemisphere)
# = less snow, drier conditions. North-facing = more concealment from shadows.

Hillshade

3D-like visualization using simulated illumination. Makes terrain features visually intuitive.

def compute_hillshade(elevation, pixel_size_m=30, azimuth=315, altitude=45):
    """
    Compute hillshade for visualization.
    azimuth: sun direction (315 = NW, standard for maps)
    altitude: sun elevation above horizon (45 deg standard)
    """
    dy, dx = np.gradient(elevation, pixel_size_m)
    slope = np.arctan(np.sqrt(dx**2 + dy**2))
    aspect = np.arctan2(-dx, dy)
 
    az_rad = np.radians(azimuth)
    alt_rad = np.radians(altitude)
 
    hillshade = (
        np.sin(alt_rad) * np.cos(slope) +
        np.cos(alt_rad) * np.sin(slope) * np.cos(az_rad - aspect)
    )
 
    return np.clip(hillshade, 0, 1) * 255
 
 
def display_terrain_analysis(elevation, pixel_size_m=30, save_path=None):
    """Complete terrain analysis visualization."""
    slope = compute_slope(elevation, pixel_size_m)
    aspect = compute_aspect(elevation, pixel_size_m)
    hillshade = compute_hillshade(elevation, pixel_size_m)
 
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
 
    # Elevation
    im0 = axes[0, 0].imshow(elevation, cmap="terrain")
    plt.colorbar(im0, ax=axes[0, 0], label="Elevation (m)")
    axes[0, 0].set_title("Elevation")
 
    # Slope
    im1 = axes[0, 1].imshow(slope, cmap="YlOrRd", vmin=0, vmax=45)
    plt.colorbar(im1, ax=axes[0, 1], label="Slope (degrees)")
    axes[0, 1].set_title("Slope")
 
    # Aspect
    im2 = axes[1, 0].imshow(aspect, cmap="hsv", vmin=0, vmax=360)
    plt.colorbar(im2, ax=axes[1, 0], label="Aspect (degrees)")
    axes[1, 0].set_title("Aspect (N=0/360, E=90, S=180, W=270)")
 
    # Hillshade
    axes[1, 1].imshow(hillshade, cmap="gray")
    axes[1, 1].set_title("Hillshade")
 
    for ax in axes.flat:
        ax.axis("off")
 
    plt.tight_layout()
    if save_path:
        plt.savefig(save_path, dpi=150)
    plt.show()
 
# Demo with synthetic terrain
np.random.seed(42)
size = 300
x, y = np.meshgrid(np.linspace(0, 5, size), np.linspace(0, 5, size))
# Ridge + valley terrain
elevation = (
    200 + 80 * np.sin(x * 1.5) * np.cos(y * 0.8) +
    40 * np.exp(-((x-2)**2 + (y-3)**2) / 0.5) +  # hill
    np.random.normal(0, 2, (size, size))  # noise
)
display_terrain_analysis(elevation, pixel_size_m=30, save_path="terrain_analysis.png")

---

Viewshed Analysis

“What can you see from this point?” — critical for observation post placement, radar coverage, communication line-of-sight, and threat assessment.

import numpy as np
 
def compute_viewshed(elevation, observer_row, observer_col, observer_height_m=2.0,
                     pixel_size_m=30, max_range_m=None):
    """
    Compute viewshed: which pixels are visible from observer position.
    Simple ray-tracing approach along radial lines.
 
    observer_height_m: height above ground (2m for standing person, 10m for tower)
    Returns: boolean array (True = visible from observer)
    """
    rows, cols = elevation.shape
    observer_elev = elevation[observer_row, observer_col] + observer_height_m
 
    visible = np.zeros((rows, cols), dtype=bool)
    visible[observer_row, observer_col] = True
 
    if max_range_m is None:
        max_range_m = max(rows, cols) * pixel_size_m
 
    max_range_pixels = int(max_range_m / pixel_size_m)
 
    # Cast rays in all directions
    n_rays = 720  # one ray every 0.5 degrees
    for angle_idx in range(n_rays):
        angle = 2 * np.pi * angle_idx / n_rays
 
        dx = np.cos(angle)
        dy = np.sin(angle)
 
        max_angle = -np.inf  # maximum elevation angle seen so far along this ray
 
        for step in range(1, max_range_pixels):
            col = int(observer_col + dx * step)
            row = int(observer_row + dy * step)
 
            if row < 0 or row >= rows or col < 0 or col >= cols:
                break
 
            # Distance from observer
            dist_m = step * pixel_size_m
 
            # Elevation angle to this pixel
            target_elev = elevation[row, col]
            elev_angle = np.arctan2(target_elev - observer_elev, dist_m)
 
            # Visible if elevation angle exceeds all previous angles along ray
            if elev_angle > max_angle:
                visible[row, col] = True
                max_angle = elev_angle
            # If blocked, pixel is not visible (stays False)
 
    return visible
 
 
def plot_viewshed(elevation, visible, observer_row, observer_col, save_path=None):
    """Visualize viewshed result."""
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))
 
    # Hillshade + viewshed overlay
    hillshade = compute_hillshade(elevation)
 
    axes[0].imshow(hillshade, cmap="gray")
    axes[0].imshow(visible, cmap="Greens", alpha=0.3)
    axes[0].plot(observer_col, observer_row, "r*", markersize=20, label="Observer")
    axes[0].set_title("Viewshed (green = visible)")
    axes[0].legend()
 
    # Elevation profile along one direction
    profile_col = np.arange(elevation.shape[1])
    profile_elev = elevation[observer_row, :]
    profile_vis = visible[observer_row, :]
 
    axes[1].fill_between(profile_col, 0, profile_elev, alpha=0.3, color="brown")
    axes[1].plot(profile_col, profile_elev, "k-", linewidth=0.5)
    vis_elev = np.where(profile_vis, profile_elev, np.nan)
    axes[1].plot(profile_col, vis_elev, "g.", markersize=2, label="Visible")
    axes[1].axvline(observer_col, color="red", linestyle="--", label="Observer")
    axes[1].set_xlabel("Column (pixels)")
    axes[1].set_ylabel("Elevation (m)")
    axes[1].set_title("E-W Profile at Observer Row")
    axes[1].legend()
 
    for ax in axes:
        ax.grid(True, alpha=0.2)
 
    plt.tight_layout()
    if save_path:
        plt.savefig(save_path, dpi=150)
    plt.show()
 
    visible_pct = np.sum(visible) / visible.size * 100
    print(f"Visible area: {visible_pct:.1f}% of total area")
 
 
# Demo
np.random.seed(42)
size = 200
x, y = np.meshgrid(np.linspace(0, 4, size), np.linspace(0, 4, size))
elevation = (
    100 + 50 * np.sin(x * 2) * np.cos(y * 1.5) +
    30 * np.exp(-((x-2)**2 + (y-2)**2) / 0.3)
)
 
observer_row, observer_col = 100, 100  # center of map
visible = compute_viewshed(elevation, observer_row, observer_col,
                            observer_height_m=5, pixel_size_m=30)
plot_viewshed(elevation, visible, observer_row, observer_col,
              save_path="viewshed_demo.png")

---

Military Terrain Analysis: OCOKA

OCOKA is the NATO framework for terrain assessment. Every factor can be analyzed with a DEM.

Observation and Fields of Fire

• Analysis: Viewshed from defensive positions and likely approach routes
• Data: DEM + viewshed algorithm
• Question: “From position X, what can be observed? Where are blind spots?”

Cover and Concealment

• Cover: protection from fire (terrain, buildings, dense vegetation)
• Concealment: protection from observation (forests, urban areas, terrain folds)
• Data: DEM (terrain folds = dead ground), NDVI (vegetation density), building footprints
• Question: “Where can forces move without being observed from position X?”

Obstacles

• Analysis: Slope > 30 degrees, water crossings, dense urban areas
• Data: DEM (slope), water bodies (NDWI), road network
• Question: “Where can vehicles NOT go?”

Key Terrain

• Analysis: High ground, chokepoints, road intersections, bridges
• Data: DEM (local maxima), road/bridge GIS data
• Question: “What terrain features, if controlled, provide a decisive advantage?”

Avenues of Approach

• Analysis: Routes with gentle slopes, good trafficability, available cover
• Data: DEM (slope < 15 deg corridors), land cover (avoid water, dense forest)
• Question: “What are the most likely movement corridors?”

def simple_ocoka_analysis(elevation, pixel_size_m=30):
    """
    Basic OCOKA terrain analysis from DEM.
    Returns dict of analysis layers.
    """
    slope = compute_slope(elevation, pixel_size_m)
 
    # Obstacles: steep terrain
    obstacles = slope > 30  # degrees
 
    # Trafficable: gentle slopes suitable for vehicles
    trafficable = slope < 15
 
    # Key terrain: local high points (higher than all neighbors within radius)
    from scipy.ndimage import maximum_filter
    radius_pixels = int(500 / pixel_size_m)  # 500m radius
    local_max = maximum_filter(elevation, size=2*radius_pixels+1)
    key_terrain = (elevation == local_max) & (slope < 20)
 
    # Dead ground: areas not visible from hilltops (require viewshed for each)
    # Simplified: areas in deep valleys
    from scipy.ndimage import minimum_filter
    local_min = minimum_filter(elevation, size=2*radius_pixels+1)
    valley_depth = elevation - local_min
    deep_valleys = valley_depth < 10  # within 10m of local minimum
 
    return {
        "slope": slope,
        "obstacles": obstacles,
        "trafficable": trafficable,
        "key_terrain": key_terrain,
        "valleys": deep_valleys,
    }

---

Geolocation from Imagery

Shadow Analysis: Measuring Object Height

If you know the sun position (azimuth and elevation) and can measure shadow length in an image, you can compute object height.

object_height = shadow_length * tan(sun_elevation)

import numpy as np
from datetime import datetime
 
def estimate_sun_position(lat, lon, dt):
    """
    Estimate sun azimuth and elevation for a given location and time.
    Simplified solar position calculation.
    """
    day_of_year = dt.timetuple().tm_yday
    hour = dt.hour + dt.minute / 60
 
    # Declination
    declination = 23.45 * np.sin(np.radians(360 / 365 * (day_of_year - 81)))
 
    # Hour angle
    solar_noon = 12 - lon / 15  # approximate (ignores equation of time)
    hour_angle = 15 * (hour - solar_noon)
 
    lat_rad = np.radians(lat)
    dec_rad = np.radians(declination)
    ha_rad = np.radians(hour_angle)
 
    # Elevation
    sin_elev = (np.sin(lat_rad) * np.sin(dec_rad) +
                np.cos(lat_rad) * np.cos(dec_rad) * np.cos(ha_rad))
    elevation = np.degrees(np.arcsin(np.clip(sin_elev, -1, 1)))
 
    # Azimuth
    cos_az = ((np.sin(dec_rad) - np.sin(lat_rad) * sin_elev) /
              (np.cos(lat_rad) * np.cos(np.radians(elevation)) + 1e-10))
    azimuth = np.degrees(np.arccos(np.clip(cos_az, -1, 1)))
    if hour_angle > 0:
        azimuth = 360 - azimuth
 
    return {"elevation_deg": elevation, "azimuth_deg": azimuth}
 
 
def height_from_shadow(shadow_length_m, sun_elevation_deg):
    """
    Compute object height from shadow length and sun elevation.
    shadow_length_m: measured shadow length on ground (in meters)
    sun_elevation_deg: sun elevation above horizon
    """
    return shadow_length_m * np.tan(np.radians(sun_elevation_deg))
 
 
def shadow_length_from_pixels(shadow_pixels, pixel_size_m):
    """Convert shadow measurement from pixels to meters."""
    return shadow_pixels * pixel_size_m
 
 
# Example: measure building height from Sentinel-2
# Image taken at 10:30 UTC over Tallinn on June 15
sun = estimate_sun_position(59.44, 24.75, datetime(2024, 6, 15, 10, 30))
print(f"Sun elevation: {sun['elevation_deg']:.1f} deg")
print(f"Sun azimuth: {sun['azimuth_deg']:.1f} deg")
 
# If we measure a shadow of 3 pixels in a 10m image
shadow_m = shadow_length_from_pixels(3, 10)
height = height_from_shadow(shadow_m, sun["elevation_deg"])
print(f"Shadow: {shadow_m}m -> Object height: {height:.1f}m")

Mountain/Horizon Matching

Geolocate a ground-level photo by matching the visible mountain silhouette against the DEM.

def generate_horizon_profile(elevation, observer_row, observer_col,
                              pixel_size_m=30, azimuth_range=(0, 360),
                              n_azimuths=360):
    """
    Generate expected horizon profile from a DEM location.
    Returns array of elevation angles for each azimuth.
    Used for matching against ground-level photographs.
    """
    rows, cols = elevation.shape
    observer_elev = elevation[observer_row, observer_col] + 1.7  # eye height
 
    azimuths = np.linspace(azimuth_range[0], azimuth_range[1], n_azimuths)
    horizon = np.zeros(n_azimuths)
 
    for i, az in enumerate(azimuths):
        az_rad = np.radians(az)
        dx = np.sin(az_rad)  # east component
        dy = -np.cos(az_rad)  # north component (row decreases northward)
 
        max_angle = -90  # minimum possible
 
        for step in range(1, max(rows, cols)):
            col = int(observer_col + dx * step)
            row = int(observer_row + dy * step)
 
            if row < 0 or row >= rows or col < 0 or col >= cols:
                break
 
            dist_m = step * pixel_size_m
            target_elev = elevation[row, col]
            angle = np.degrees(np.arctan2(target_elev - observer_elev, dist_m))
 
            if angle > max_angle:
                max_angle = angle
 
        horizon[i] = max(max_angle, 0)
 
    return azimuths, horizon

---

Exercises

Exercise 1: Compute Viewshed from a Hilltop

1. Download SRTM DEM for an area you know (your region)
2. Choose a high point (hilltop, tower location)
3. Compute viewshed with observer height = 10m (tower)
4. What percentage of the surrounding 10km is visible?
5. Where are the major blind spots?

Exercise 2: Determine Object Height from Shadow

1. Find a satellite image where building shadows are visible (morning or evening image preferred — longer shadows)
2. Determine the image acquisition time (from metadata)
3. Compute sun position for that time and location
4. Measure shadow length in pixels, convert to meters
5. Estimate building height — verify against known data if possible

Exercise 3: Terrain Assessment for Observation Post

1. Given a 10x10 km area DEM, identify the 3 best locations for an observation post
2. Criteria: high ground, covers the most area, has cover (forest nearby)
3. Compute viewshed from each candidate position
4. Present as a comparison: which position has the best observation coverage?

---

Self-Test Questions

1. What is the difference between DSM and DTM? Which does SRTM provide?
2. A slope of 25 degrees — is it passable for wheeled vehicles? Tracked vehicles?
3. How does viewshed analysis help in communication planning?
4. At sun elevation 30 degrees, a shadow is 20 meters long. How tall is the object?
5. Why is OCOKA analysis relevant to both offensive and defensive operations?

---

See also: Satellite Fundamentals | Change Detection | Case Study - Monitoring Military Installations Next: Multi-Source Intelligence Fusion