Chapter 7 Drone Remote Pilot: Geographic Information System#

Introduction: GIS
Simulation: Geographic and Map Projection
Simulation: DEM analysis and Visualization
Self-Assessment

1. Introduction#

🧭 What is GIS?#

Geographic Information System (GIS) is a framework for capturing, storing, analyzing, and visualizing spatial and geographic data.
It enables users to explore patterns, relationships, and trends through location-based insights using maps, layers, and queries.

🗂️ Data Models in GIS#

📦 Raster Data Model#

Represents the world as a grid of cells or pixels
Each cell holds a value (e.g., elevation, temperature, land cover)
Common in remote sensing and continuous data representation
Example: Satellite imagery, DEMs

🧩 Vector Data Model#

Represents spatial features using points, lines, and polygons
Points: Locations (e.g., wells, trees)
Lines: Roads, rivers
Polygons: Areas (e.g., land parcels, lakes)
Ideal for discrete and detailed geographic objects

🌐 Geographic vs. Map Projection#

🌍 Geographic Coordinate System (GCS)#

Based on latitude and longitude
Spherical model of Earth
Example: WGS84 (World Geodetic System 1984)

🗺️ Map Projection#

Transforms 3D Earth surface into a 2D map
Preserves one or more properties (e.g., area, shape, distance)
Common Projections:
- UTM (Universal Transverse Mercator): Divides Earth into zones, preserves shapes and distances
- Mercator: Great for navigation, distorts poles
- Albers Equal-Area: Used for thematic maps

🧮 Raster Analysis Using Map Algebra#

Map Algebra enables mathematical operations on raster layers.
Each raster is treated as a matrix of values and operations are applied cell-by-cell.

🧪 Types of Operations#

Local: Cell-by-cell (e.g., A + B, log(C))
Focal: Neighborhood-based (e.g., moving average, slope)
Zonal: Aggregates values by zones (e.g., mean elevation per watershed)
Global: Entire raster analysis (e.g., total area, histogram)

📊 Example#

# Using raster calculator logic
ElevationCorrected = ElevationRaster + OffsetRaster
VegetationMask = NDVI > 0.3

References#

[Bolstad, 2019]suitable for civil engineering applications — especially in remote sensing, drone-based mapping, and spatial analysis. The book introduces GIS concepts, Raster and Vector Data analysis, Spatial interpolation and Surface modeling, Field Data integration, and remote sensing foundations.

2. Simulation#

🌍 Interactive Projection Comparison with ipywidgets & Cartopy#

This Python code enables users to compare two cartographic projections by overlaying them in a visually distinct way using matplotlib, Cartopy, and ipywidgets.

🔧 Key Libraries#

ipywidgets: For creating interactive dropdown menus
matplotlib.pyplot: For rendering maps
cartopy.crs: Defines coordinate reference systems (CRS)
cartopy.feature: Adds geographic features like coastlines, borders, and states

🗺️ What Does It Do?#

Users select two map projections from dropdowns.
The script dynamically generates a map showing:
- Projection A in blue (α = 0.6)
- Projection B in red (α = 0.4)
Geographic features (states, borders, coastline, land/ocean) are displayed with transparency to emphasize differences.
The map focuses on the continental US with an extent of [-130, -65, 23, 50].

🎛️ How It Works#

Dropdown Menus
- Users pick two projections from an expanded set (e.g., Mercator, Albers, Lambert, Orthographic).
update_map Function
- Responds to dropdown changes.
- Clears previous output.
- Creates ax with Projection A.
- Overlays ax2 with Projection B.
- Adds features to visualize boundary differences.
Interactivity
- observe(update_map, names="value"): Attaches event listeners to dropdowns.
- display(VBox([...])): Renders dropdowns and map output vertically.
- update_map(): Initializes the map display on first load.

🎨 Example Use Case#

Compare how state boundaries appear when projected using:

WGS84 (EPSG:4326) vs. Lambert Conformal Conic
Orthographic vs. Azimuthal Equidistant

This helps geospatial analysts and students understand how projections affect shape, area, and relative distances on maps.

💡 Bonus Suggestion#

You could enhance this visualization by:

Adding tooltips or hover descriptions for each projection
Including country boundaries for international comparisons
Exporting selected projection data for further GIS analysis

import ipywidgets as widgets
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 🌐 Expanded projection options
projection_options = {
    "WGS84 (EPSG:4326)": ccrs.PlateCarree(),
    "Mercator": ccrs.Mercator(),
    "Transverse Mercator": ccrs.TransverseMercator(central_longitude=-96),
    "UTM Zone 15N": ccrs.UTM(zone=15),
    "Albers Equal Area": ccrs.AlbersEqualArea(central_longitude=-96, central_latitude=23),
    "Lambert Conformal Conic": ccrs.LambertConformal(central_longitude=-96, central_latitude=23),
    "Orthographic": ccrs.Orthographic(central_longitude=-96, central_latitude=39),
    "Azimuthal Equidistant": ccrs.AzimuthalEquidistant(central_longitude=-96, central_latitude=39)
}

# 🎛️ Dropdown selectors
dropdown1 = widgets.Dropdown(
    options=list(projection_options.keys()),
    value="WGS84 (EPSG:4326)",
    description='Projection A:'
)
dropdown2 = widgets.Dropdown(
    options=list(projection_options.keys()),
    value="Albers Equal Area",
    description='Projection B:'
)

output = widgets.Output()

# 🗺️ Update map function
def update_map(change=None):
    output.clear_output()
    
    proj_a = dropdown1.value
    proj_b = dropdown2.value
    crs_a = projection_options[proj_a]
    crs_b = projection_options[proj_b]
    
    with output:
        print(f"\n🟦 {proj_a} (blue, α=0.6)\n🟥 {proj_b} (red, α=0.4)")
        
        fig = plt.figure(figsize=(10, 8))
        ax = plt.axes(projection=crs_a)
        ax.set_extent([-130, -65, 23, 50], crs=ccrs.PlateCarree())

        # Base: Projection A (semi-transparent)
        ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='lightgray', alpha=0.3)
        ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='lightblue', alpha=0.3)
        ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='blue', linewidth=1.2, zorder=1, alpha=0.6)
        ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='--', alpha=0.6)
        ax.add_feature(cfeature.COASTLINE.with_scale('50m'), alpha=0.6)

        # Overlay: Projection B (lighter transparent)
        ax2 = plt.axes(projection=crs_b)
        ax2.set_extent([-130, -65, 23, 50], crs=ccrs.PlateCarree())
        ax2.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='red', linewidth=0.8, zorder=2, alpha=0.4)

        ax.set_title(f'Overlay of State Boundaries\n{proj_a} vs. {proj_b}', fontsize=14)
        plt.show()

# 🔄 Interaction
dropdown1.observe(update_map, names="value")
dropdown2.observe(update_map, names="value")
display(widgets.VBox([dropdown1, dropdown2, output]))
update_map()

3. Simulation#

🌄 DEM Analysis and Visualization Tool#

This interactive dashboard enables users to explore a synthetic Digital Elevation Model (DEM) with dynamic controls for:

Grid resolution adjustment via slider for initial detail level
Aggregation factor to simulate spatial coarsening and resampling
Analysis types including:
- 🏞️ Raw DEM visualization
- ↗️ Slope computation
- 🟢 Elevation masking (>120m)
- 🧭 Elevation classification into bands

Users can intuitively observe how aggregation and resolution affect terrain representation and derivative spatial analyses. Grid transformations and visualizations update in real time to enhance geospatial understanding.

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from scipy.ndimage import generic_filter, zoom

# 🎚️ Sliders for resolution and aggregation level
resolution_slider = widgets.IntSlider(
    value=100, min=50, max=300, step=25,
    description='Initial Grid Resolution:',
    style={'description_width': 'initial'},
    continuous_update=False
)

aggregation_slider = widgets.IntSlider(
    value=1, min=1, max=5, step=1,
    description='Aggregation Factor:',
    style={'description_width': 'initial'},
    continuous_update=False
)

# 📊 Dropdown for analysis type
analysis_dropdown = widgets.Dropdown(
    options=['Raw DEM', 'Slope', 'Elevation Mask (>120m)', 'Classified Elevation'],
    description='Analysis Type:',
    value='Raw DEM',
    style={'description_width': 'initial'}
)

output = widgets.Output()

# ⚙️ Slope calculator
def compute_slope(raster):
    def slope_func(window):
        center = window[len(window)//2]
        diffs = np.abs(window - center)
        return np.mean(diffs)
    return generic_filter(raster, slope_func, size=3)

# 🔄 Interactive update function
def update_analysis(change=None):
    output.clear_output()
    res = resolution_slider.value
    agg = aggregation_slider.value
    method = analysis_dropdown.value

    # 🔧 Create raw synthetic DEM
    np.random.seed(42)
    base = np.random.normal(loc=100, scale=20, size=(res, res))
    hill = np.exp(-((np.linspace(-1,1,res)**2)[:,None] + (np.linspace(-1,1,res)**2)[None,:])) * 50
    dem = base + hill

    # 🔁 Apply aggregation by resizing with averaging
    if agg > 1:
        scale = 1 / agg
        dem = zoom(dem, zoom=scale, order=1)  # bilinear downsampling
        dem = zoom(dem, zoom=agg, order=1)    # upscale to original size for display

    with output:
        plt.figure(figsize=(6, 5))
        
        if method == 'Raw DEM':
            plt.imshow(dem, cmap='terrain')
            plt.title(f'🏞️ Aggregated DEM (Factor {agg})')
        elif method == 'Slope':
            slope = compute_slope(dem)
            plt.imshow(slope, cmap='gray')
            plt.title(f'↗️ Slope from Aggregated DEM (Factor {agg})')
        elif method == 'Elevation Mask (>120m)':
            mask = dem > 120
            plt.imshow(mask, cmap='Greens')
            plt.title(f'🟢 Elevation Mask >120m (Factor {agg})')
        elif method == 'Classified Elevation':
            classes = np.zeros_like(dem)
            classes[dem < 90] = 1
            classes[(dem >= 90) & (dem < 110)] = 2
            classes[dem >= 110] = 3
            plt.imshow(classes, cmap='viridis')
            plt.title(f'🧭 Classified DEM (Factor {agg})')

        plt.colorbar(label='Value')
        plt.axis('off')
        plt.show()

# 🧩 Hook up interactivity
for widget in [resolution_slider, aggregation_slider, analysis_dropdown]:
    widget.observe(update_analysis, names='value')

# 🚀 Launch dashboard
display(widgets.VBox([resolution_slider, aggregation_slider, analysis_dropdown, output]))
update_analysis()

Chapter 7 Drone Remote Pilot: Geographic Information System

Contents