python refresher for spatial analysis

Warm-up before we touch a single raster or shapefile

Before we open a single satellite image or elevation model, we need to make sure everyone speaks the same Python language. This session is not about GIS itself. It is about the small set of Python tools that every single session in this course will lean on. Nothing here is random. Every example was picked because you will see the exact same pattern again later, when you calculate NDVI, slope, camera parameters, or filter road networks.

How to use this notebook

Read the explanation, run the example, then try the exercise before looking at the solution. Solutions are hidden inside a collapsible box. Click to reveal them only after you have tried on your own.

Variables, Lists, and Dictionaries

In spatial work, almost everything is a bunch of numbers with a label attached to it: pixel values, band names, coordinate systems, ground control point coordinates, weights used in a model. Lists and dictionaries are the two basic Python structures used to organize that kind of information.

Variables

A variable is simply a name that points to a value. The value can be a number, text, or something more complex. Let’s start with something familiar from surveying: a single ground control point (GCP) used in photogrammetry to align an aerial photo with real world coordinates.

gcp_name = "GCP-01"
gcp_x = 512340.25
gcp_y = 3891220.10
gcp_elevation = 812.4

print(gcp_name, "is located at elevation", gcp_elevation, "meters")

GCP-01 is located at elevation 812.4 meters

# Example 2: a camera parameter used in photogrammetry
focal_length_mm = 35
sensor_width_mm = 23.5
image_width_px = 6000

# Ground sample distance depends on all three parameters
flying_height_m = 1200
gsd_cm = (sensor_width_mm * flying_height_m * 100) / (focal_length_mm * image_width_px)
print("Ground sample distance:", round(gsd_cm, 2), "cm/pixel")

Ground sample distance: 13.43 cm/pixel

Lists

A list stores an ordered collection of values. This is the natural way to represent a series of readings, such as elevation values collected at several ground control points, or NDVI values from a set of field plots.

# Elevation readings (in meters) from 6 ground control points
elevations = [812, 790, 845, 803, 860, 798]

print("Number of GCPs:", len(elevations))
print("Highest point:", max(elevations))
print("Lowest point:", min(elevations))
print("Total:", sum(elevations))

Number of GCPs: 6
Highest point: 860
Lowest point: 790
Total: 4908

Lists come with several built-in methods that you will use constantly. Below are the ones that matter most for this course.

band_names = ["red", "nir"]

# add a single item
band_names.append("swir")
print(band_names)

# add several items at once
band_names.extend(["blue", "green"])
print(band_names)

# remove a specific value
band_names.remove("swir")
print(band_names)

# sort alphabetically
band_names.sort()
print(band_names)

# check whether a value exists
print("red" in band_names)

# find the position of a value
print(band_names.index("nir"))

['red', 'nir', 'swir']
['red', 'nir', 'swir', 'blue', 'green']
['red', 'nir', 'blue', 'green']
['blue', 'green', 'nir', 'red']
True
2

# Slicing: getting a subset of a list
flight_line_ids = ["L01", "L02", "L03", "L04", "L05", "L06"]

print(flight_line_ids[0])      # first item
print(flight_line_ids[-1])     # last item
print(flight_line_ids[1:4])    # items 1 to 3
print(flight_line_ids[::2])    # every second item

L01
L06
['L02', 'L03', 'L04']
['L01', 'L03', 'L05']

Looping through a list

elevations = [812, 790, 845, 803, 860, 798]

for index, value in enumerate(elevations):
    print(f"GCP {index + 1}: {value} m")

GCP 1: 812 m
GCP 2: 790 m
GCP 3: 845 m
GCP 4: 803 m
GCP 5: 860 m
GCP 6: 798 m

# A second example, from remote sensing: cloud cover percentage per scene
cloud_cover = [12, 45, 3, 67, 8, 30]
scene_ids = ["S01", "S02", "S03", "S04", "S05", "S06"]

usable_scenes = []
for scene_id, cover in zip(scene_ids, cloud_cover):
    if cover < 20:
        usable_scenes.append(scene_id)

print("Usable scenes (cloud cover below 20 percent):", usable_scenes)

Usable scenes (cloud cover below 20 percent): ['S01', 'S03', 'S05']

Dictionaries

A dictionary stores key and value pairs. Instead of creating several separate variables, a dictionary keeps related information together under one name. In Session 3, rasterio will return raster metadata in exactly this format.

raster_metadata = {
    "filename": "sentinel.tif",
    "width": 512,
    "height": 512,
    "crs": "EPSG:32639",
    "bands": 4
}

print("This raster uses CRS:", raster_metadata["crs"])
print("Dimensions:", raster_metadata["width"], "x", raster_metadata["height"])

This raster uses CRS: EPSG:32639
Dimensions: 512 x 512

# Dictionary methods that matter most
print(raster_metadata.keys())
print(raster_metadata.values())
print(raster_metadata.items())

# safely get a value that might not exist
print(raster_metadata.get("nodata", "not specified"))

# update or add a key
raster_metadata["nodata"] = -9999
print(raster_metadata)

dict_keys(['filename', 'width', 'height', 'crs', 'bands'])
dict_values(['sentinel.tif', 512, 512, 'EPSG:32639', 4])
dict_items([('filename', 'sentinel.tif'), ('width', 512), ('height', 512), ('crs', 'EPSG:32639'), ('bands', 4)])
not specified
{'filename': 'sentinel.tif', 'width': 512, 'height': 512, 'crs': 'EPSG:32639', 'bands': 4, 'nodata': -9999}

# A nested dictionary example: metadata for several photogrammetry cameras
cameras = {
    "camera_A": {"focal_length_mm": 35, "sensor_width_mm": 23.5},
    "camera_B": {"focal_length_mm": 50, "sensor_width_mm": 36.0}
}

for camera_name, specs in cameras.items():
    print(f"{camera_name}: focal length {specs['focal_length_mm']} mm")

camera_A: focal length 35 mm
camera_B: focal length 50 mm

Note: dictionaries versus lists

Use a list when order matters and items are of the same kind, such as a series of elevation readings. Use a dictionary when you need to label each value, such as raster metadata or camera specifications.

Exercise 1.1

Create a dictionary called image_info with keys sensor, resolution_m, and cloud_cover_percent. Fill in realistic values. Then print a sentence reporting the sensor name and its resolution using an f-string.

Show solution

image_info = {
    "sensor": "Sentinel-2",
    "resolution_m": 10,
    "cloud_cover_percent": 5
}

print(f"{image_info['sensor']} has a resolution of {image_info['resolution_m']} m")

Sentinel-2 has a resolution of 10 m

Exercise 1.2

Given elevations = [812, 790, 845, 803, 860, 798], write a loop that prints only the ground control points with elevation above 800 m, along with their index in the list.

Show solution

elevations = [812, 790, 845, 803, 860, 798]

for index, value in enumerate(elevations):
    if value > 800:
        print(f"GCP {index + 1}: {value} m")

GCP 1: 812 m
GCP 3: 845 m
GCP 4: 803 m
GCP 5: 860 m

Exercise 1.3

Starting from band_names = ["blue", "green", "red"], add "nir" and "swir" using two different list methods, then sort the final list and print it.

Show solution

band_names = ["blue", "green", "red"]
band_names.append("nir")
band_names.extend(["swir"])
band_names.sort()
print(band_names)

['blue', 'green', 'nir', 'red', 'swir']

Quick quiz

Question 1

Which method would you use to safely retrieve a value from a dictionary without raising an error if the key does not exist?

A. dictionary.keys()

B. dictionary.get()

C. dictionary.values()

D. dictionary.items()

Answer

B. The .get() method returns a default value instead of raising an error when the key is missing, which is exactly what raster_metadata.get("nodata", "not specified") did above.

Functions

Almost every session in this course ends with writing one function, such as calculate_ndvi, calculate_slope, or calculate_final_score. Becoming comfortable with functions now will remove a lot of confusion later.

A basic function

def celsius_to_fahrenheit(celsius):
    """Convert a temperature from Celsius to Fahrenheit."""
    fahrenheit = (celsius * 9 / 5) + 32
    return fahrenheit

print(celsius_to_fahrenheit(21.5))

70.7

Functions with default parameters

def classify_elevation(elevation_m, threshold=800):
    """Return 'high' or 'low' depending on a threshold."""
    if elevation_m > threshold:
        return "high"
    return "low"

print(classify_elevation(860))
print(classify_elevation(790))
print(classify_elevation(860, threshold=900))

high
low
low

Note: default parameters

threshold=800 means use 800 unless the caller specifies otherwise. This exact pattern appears again in Session 8, where the suitability function has default weights that a student can later override.

Functions that return multiple values

def image_stats(pixel_values):
    """Return the minimum, maximum, and average of a list of pixel values."""
    minimum = min(pixel_values)
    maximum = max(pixel_values)
    average = sum(pixel_values) / len(pixel_values)
    return minimum, maximum, average

pixels = [120, 98, 145, 132, 110, 156]
low, high, avg = image_stats(pixels)
print("Min:", low, "Max:", high, "Average:", round(avg, 2))

Min: 98 Max: 156 Average: 126.83

Functions that call other functions

def dn_to_reflectance(dn_value, gain=0.0001, offset=0.0):
    """Convert a digital number to reflectance using a linear model."""
    return dn_value * gain + offset

def classify_reflectance(reflectance, threshold=0.2):
    """Classify a reflectance value as vegetation or non-vegetation."""
    if reflectance > threshold:
        return "vegetation"
    return "non-vegetation"

digital_number = 2500
reflectance_value = dn_to_reflectance(digital_number)
label = classify_reflectance(reflectance_value)

print(f"DN {digital_number} corresponds to reflectance {reflectance_value:.3f}, classified as {label}")

DN 2500 corresponds to reflectance 0.250, classified as vegetation

Applying a function to a whole list

def classify_all(elevation_list, threshold=800):
    """Classify every elevation value in a list."""
    results = []
    for value in elevation_list:
        results.append(classify_elevation(value, threshold))
    return results

elevations = [812, 790, 845, 803, 860, 798]
print(classify_all(elevations))

['high', 'low', 'high', 'high', 'high', 'low']

Exercise 2.1

Write a function meters_to_feet(meters) that converts meters to feet, using the fact that 1 meter equals 3.281 feet. Test it with the flying height used earlier in this session.

Show solution

def meters_to_feet(meters):
    """Convert meters to feet."""
    return meters * 3.281

print(meters_to_feet(1200))

3937.2000000000003

Exercise 2.2

Write a function average_elevation(elevation_list) that returns the average of a list of elevation values, then test it on the elevations list used above.

Show solution

def average_elevation(elevation_list):
    """Return the average value of a list of numbers."""
    return sum(elevation_list) / len(elevation_list)

elevations = [812, 790, 845, 803, 860, 798]
print(average_elevation(elevations))

818.0

Exercise 2.3

Write a function gsd_from_flight(focal_length_mm, sensor_width_mm, image_width_px, flying_height_m) that returns the ground sample distance in centimeters, based on the formula used at the start of this session. Test it with two different flying heights.

Show solution

def gsd_from_flight(focal_length_mm, sensor_width_mm, image_width_px, flying_height_m):
    """Return ground sample distance in cm per pixel."""
    return (sensor_width_mm * flying_height_m * 100) / (focal_length_mm * image_width_px)

print(gsd_from_flight(35, 23.5, 6000, 1200))
print(gsd_from_flight(35, 23.5, 6000, 800))

13.428571428571429
8.952380952380953

Quick quiz

Question 2

What does a function return if there is no explicit return statement inside it?

A. An error, always

B. None

C. 0

D. The last variable defined inside the function

Answer

B. Python functions return None by default when there is no explicit return statement.

From Lists to NumPy Arrays

This is the single most important concept in this session. A raster, whether it is a satellite image or an elevation model, is nothing more than a large NumPy array of numbers. Every calculation performed in later sessions, including NDVI, slope, aspect, and the final suitability score, is simply NumPy math applied to arrays.

The problem with plain lists

elevations = [812, 790, 845, 803, 860, 798]

# converting every value to feet forces a manual loop with a plain list
elevations_ft = []
for value in elevations:
    elevations_ft.append(value * 3.281)

print(elevations_ft)

[2664.172, 2591.9900000000002, 2772.445, 2634.643, 2821.6600000000003, 2618.2380000000003]

The NumPy way

import numpy as np

elevations = np.array([812, 790, 845, 803, 860, 798])

# multiplication now applies to every value at once, no loop needed
elevations_ft = elevations * 3.281
print(elevations_ft)

[2664.172 2591.99  2772.445 2634.643 2821.66  2618.238]

Note: vectorized operations

When an operation applies to a whole array at once, without an explicit loop, this is called a vectorized operation. This is not only shorter code, it is dramatically faster, which matters when the array is an image with millions of pixels.

Creating arrays with built-in functions

zeros_array = np.zeros(5)
ones_array = np.ones(5)
sequence_array = np.arange(0, 50, 10)
even_spacing = np.linspace(0, 1, 5)

print("Zeros:", zeros_array)
print("Ones:", ones_array)
print("Arange:", sequence_array)
print("Linspace:", even_spacing)

Zeros: [0. 0. 0. 0. 0.]
Ones: [1. 1. 1. 1. 1.]
Arange: [ 0 10 20 30 40]
Linspace: [0.   0.25 0.5  0.75 1.  ]

These functions are commonly used to create placeholder arrays before filling them with real raster values, or to generate evenly spaced coordinates along a flight line.

Filtering with conditions

elevations = np.array([812, 790, 845, 803, 860, 798])

# boolean mask: True or False for each value
above_800 = elevations > 800
print(above_800)

# use the mask to pull out the matching values
print(elevations[above_800])

# equivalent one line version
print(elevations[elevations > 800])

[ True False  True  True  True False]
[812 845 803 860]
[812 845 803 860]

Note: this is exactly how we will pick suitable pixels

In Session 5, code very similar to slope_array[slope_array < 10] will be used to find every pixel with a gentle enough slope for solar panel installation. It is the same pattern used above.

Basic statistics

ndvi_values = np.array([0.1, 0.65, 0.72, 0.05, 0.55, 0.68])

print("Mean:", np.mean(ndvi_values))
print("Median:", np.median(ndvi_values))
print("Standard deviation:", np.std(ndvi_values))
print("Minimum:", np.min(ndvi_values))
print("Maximum:", np.max(ndvi_values))
print("Sum:", np.sum(ndvi_values))
print("Count above 0.5:", np.sum(ndvi_values > 0.5))

Mean: 0.4583333333333334
Median: 0.6000000000000001
Standard deviation: 0.2762496857715659
Minimum: 0.05
Maximum: 0.72
Sum: 2.7500000000000004
Count above 0.5: 4

Two dimensional arrays

A raster is a grid of numbers, which in NumPy is represented as a two dimensional array. Imagine the grid below as a small elevation model.

dem = np.array([
    [812, 815, 818, 820, 825],
    [805, 808, 812, 818, 822],
    [795, 798, 802, 808, 815],
    [780, 785, 790, 796, 805],
    [760, 770, 775, 782, 790]
])

print("Shape of this DEM:", dem.shape)
print("Value at row 2, column 3:", dem[2, 3])

Shape of this DEM: (5, 5)
Value at row 2, column 3: 808

Note: rows versus columns

In NumPy, and in every raster used later, the first index is the row (the y direction, top to bottom) and the second index is the column (the x direction, left to right). This is one of the most common beginner mistakes, so keep it in mind starting from Session 3.

Slicing a two dimensional array

print(dem[0:2, 0:2])   # top-left 2x2 corner
print(dem[0, :])        # entire first row
print(dem[:, -1])       # entire last column
print(dem[-1, :])       # entire last row

[[812 815]
 [805 808]]
[812 815 818 820 825]
[825 822 815 805 790]
[760 770 775 782 790]

Applying conditions to a two dimensional array

high_pixels = dem > 800
print(high_pixels)

print("Count of pixels above 800 m:", np.sum(high_pixels))
print("Mean elevation of the whole grid:", np.mean(dem))
print("Mean elevation of only the high pixels:", np.mean(dem[high_pixels]))

[[ True  True  True  True  True]
 [ True  True  True  True  True]
 [False False  True  True  True]
 [False False False False  True]
 [False False False False False]]
Count of pixels above 800 m: 14
Mean elevation of the whole grid: 800.24
Mean elevation of only the high pixels: 813.2142857142857

Reshaping and combining arrays

flat_values = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])

grid = flat_values.reshape(3, 3)
print(grid)

back_to_flat = grid.flatten()
print(back_to_flat)

[[1 2 3]
 [4 5 6]
 [7 8 9]]
[1 2 3 4 5 6 7 8 9]

# A remote sensing example: combining two bands into one array
red_band = np.array([[80, 85], [90, 95]])
nir_band = np.array([[150, 160], [140, 130]])

stacked = np.stack([red_band, nir_band])
print("Stacked shape:", stacked.shape)
print(stacked)

Stacked shape: (2, 2, 2)
[[[ 80  85]
  [ 90  95]]

 [[150 160]
  [140 130]]]

Exercise 3.1

Using the dem array defined above, find the minimum and maximum elevation in the grid, and count how many pixels have elevation below 790 m.

Show solution

print("Min:", np.min(dem))
print("Max:", np.max(dem))
print("Pixels below 790:", np.sum(dem < 790))

Min: 760
Max: 825
Pixels below 790: 6

Exercise 3.2

Create a NumPy array from the NDVI values [0.12, 0.58, 0.61, 0.09, 0.70, 0.33]. Print only the values that represent healthy vegetation, using a threshold of 0.5, and print how many values passed the threshold.

Show solution

ndvi = np.array([0.12, 0.58, 0.61, 0.09, 0.70, 0.33])
healthy = ndvi[ndvi > 0.5]
print(healthy)
print("Count:", len(healthy))

[0.58 0.61 0.7 ]
Count: 3

Exercise 3.3

Create a one dimensional array with np.arange containing the numbers from 0 to 15, reshape it into a 4 by 4 grid, then print only the last column of that grid.

Show solution

values = np.arange(16)
grid = values.reshape(4, 4)
print(grid)
print(grid[:, -1])

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]
[ 3  7 11 15]

Quick quiz

Question 3

Why are NumPy arrays preferred over plain Python lists when working with raster data?

A. NumPy arrays use less disk space than lists

B. NumPy arrays support fast, vectorized math without explicit loops

C. Lists cannot store numbers, only text

D. NumPy arrays automatically carry a coordinate system

Answer

B. The real power of NumPy for this course is vectorized operations: applying math to millions of pixels at once with no loop required. Option D is tempting but incorrect. A plain NumPy array has no geographic information attached; that is what a library such as rasterio adds later.

Working with Files and Paths

Starting from Session 1, the project will have a real folder structure, including data, data/outputs, and modules. A reliable way of building file paths is needed that works the same way on Windows, Mac, and Linux.

The pathlib way

from pathlib import Path

data_folder = Path("data")
outputs_folder = data_folder / "outputs"
ndvi_path = outputs_folder / "ndvi.tif"

print(ndvi_path)
print("Does this file exist:", ndvi_path.exists())

data/outputs/ndvi.tif
Does this file exist: True

Note: the slash operator

The / symbol joins path pieces together when used with pathlib. It automatically uses the correct separator for the operating system, so the same code works on any computer.

Useful path attributes and methods

image_path = Path("data/sentinel_scene_01.tif")

print("Name:", image_path.name)
print("Stem (name without extension):", image_path.stem)
print("Suffix (extension):", image_path.suffix)
print("Parent folder:", image_path.parent)
print("Is this an absolute path:", image_path.is_absolute())

Name: sentinel_scene_01.tif
Stem (name without extension): sentinel_scene_01
Suffix (extension): .tif
Parent folder: data
Is this an absolute path: False

Listing and creating folders

data_folder = Path("data")

if data_folder.exists():
    for file in data_folder.glob("*.tif"):
        print(file.name)
else:
    print("Folder 'data' does not exist yet in this environment")

dem.tif
sentinel.tif

outputs_folder = Path("data/outputs")
outputs_folder.mkdir(parents=True, exist_ok=True)
print("Folder ready:", outputs_folder)

Folder ready: data/outputs

Note: exist_ok=True

Without this argument, Python raises an error if the folder already exists. With it, Python continues without complaining, which is useful when a notebook is run more than once.

Checking file type

sample_paths = [
    Path("data/sentinel.tif"),
    Path("data/roads.shp"),
    Path("data/notes.txt")
]

for path in sample_paths:
    print(path.name, "->", path.suffix)

sentinel.tif -> .tif
roads.shp -> .shp
notes.txt -> .txt

Exercise 4.1

Using pathlib, build a path pointing to data/roads.shp without typing the slash character manually, then print its stem and suffix.

Show solution

from pathlib import Path

roads_path = Path("data") / "roads.shp"
print(roads_path)
print(roads_path.stem)
print(roads_path.suffix)

data/roads.shp
roads
.shp

Building and Importing Your Own Module

Every session produces a .py file inside a modules folder, such as data_loader.py or raster_utils.py. By the end of the course, main.py and app.py import and combine all of them. This section practices the same pattern with something simple.

Creating a small module

Save the code below in a file named unit_converter.py, in the same folder as this notebook.

# unit_converter.py

def meters_to_feet(meters):
    """Convert meters to feet."""
    return meters * 3.281

def celsius_to_fahrenheit(celsius):
    """Convert Celsius to Fahrenheit."""
    return (celsius * 9 / 5) + 32

def classify_elevation(elevation_m, threshold=800):
    """Return 'high' or 'low' depending on a threshold."""
    if elevation_m > threshold:
        return "high"
    return "low"

Importing the whole module

import unit_converter

print(unit_converter.meters_to_feet(1200))
print(unit_converter.celsius_to_fahrenheit(21.5))

Importing specific functions

from unit_converter import meters_to_feet, classify_elevation

print(meters_to_feet(500))
print(classify_elevation(860))

Note: where the module file should live

Python looks for unit_converter.py in the same folder as the notebook, or somewhere on its search path. If a ModuleNotFoundError appears, the most common cause is the file being in the wrong folder. Starting Session 1, this will be organized properly with a dedicated modules folder.

Exercise 5.1

Add a new function dn_to_reflectance(dn_value, gain=0.0001, offset=0.0) into unit_converter.py, matching the one defined earlier in this session. Import it and test it with two different digital number values.

Show solution

def dn_to_reflectance(dn_value, gain=0.0001, offset=0.0):
    """Convert a digital number to reflectance using a linear model."""
    return dn_value * gain + offset

print(dn_to_reflectance(2500))
print(dn_to_reflectance(1800))

0.25
0.18000000000000002

Quick quiz

Question 4

If the statement from unit_converter import meters_to_feet is used, how is the function called afterward?

A. unit_converter.meters_to_feet(500)

B. meters_to_feet(500)

C. import.meters_to_feet(500)

D. It cannot be called this way, this import style is only for classes

Answer

B. When a specific function is imported with from module import function_name, it is called directly by name, with no module prefix needed. Option A would be correct only if import unit_converter had been used instead.

Summary

This session covered every Python tool the rest of the course will lean on: lists and dictionaries for storing metadata and small collections of values, functions as the shape every processing step will take, NumPy arrays as the real structure behind every raster, pathlib for navigating the project folders, and custom modules as the pattern the entire ten session project is built on.

Coming up in Session 01

A real satellite image will be opened with rasterio, showing for the first time that it really is just a NumPy array with extra geographic information attached.