Chapter 2: The Geospatial Python Toolbox¶

Part I — Scientific Python for Researchers

Learning Objectives¶

By the end of this chapter you will be able to:

Store and manipulate geographic values using Python's core data types
Write conditional logic to classify locations based on their coordinates
Organize multiple records with lists, tuples, and dictionaries
Automate repetitive work with for loops and your own functions
Combine these building blocks into a small, working rule-based model
Recognize when a vectorized NumPy operation should replace a Python loop

In [1]:

Copied!





# Standard imports — add chapter-specific imports below
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Standard imports — add chapter-specific imports below
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Data as Coordinates¶

In the world of computing, data refers to any collection of values or information that can be processed, stored, and transformed by a computer. This data can represent a wide variety of information, including text, images, and coordinates. In reality, any information processed by a computer must eventually be translated into numbers, as computers remain, at their core, sophisticated machines for numerical processing. In geography, everything starts with a location. To a computer, a location is simply a set of numbers stored in variables. A variable acts as a labeled container for data. For example, if we want to calculate the distance between two points, we first need to store their coordinates as floats (decimal numbers).

To illustrate this, let's represent the locations of two iconic cities: Rio de Janeiro and Paris. The following code demonstrates how we assign these numerical values to specific variables:

In [2]:

Copied!





# Assignment: Placing the coordinates of Rio and Paris into variables
# Rio de Janeiro, Brazil
x1 = -43.1729
y1 = -22.9068

# Paris, France
x2 = 2.3522
y2 = 48.8566

# Now the computer has these 'values' stored in 'labeled containers'
# Assignment: Placing the coordinates of Rio and Paris into variables
# Rio de Janeiro, Brazil
x1 = -43.1729
y1 = -22.9068

# Paris, France
x2 = 2.3522
y2 = 48.8566

# Now the computer has these 'values' stored in 'labeled containers'

Did you know?

In Python, the equals sign (=) does not represent mathematical equality. Instead, it is an assignment operator. Think of it as an instruction that says: 'Take the value on the right and store it inside the labeled container on the left.'

In [3]:

Copied!





# Every variable has a type. Python infers it from the value you assign.
print(type(x1))   # <class 'float'>
print(type(y1))

# Coordinates are just numbers, so ordinary arithmetic applies to them
delta_lon = x2 - x1
delta_lat = y2 - y1

print(f"Difference in longitude: {delta_lon:.2f} degrees")
print(f"Difference in latitude: {delta_lat:.2f} degrees")
# Every variable has a type. Python infers it from the value you assign.
print(type(x1))   # <class 'float'>
print(type(y1))

# Coordinates are just numbers, so ordinary arithmetic applies to them
delta_lon = x2 - x1
delta_lat = y2 - y1

print(f"Difference in longitude: {delta_lon:.2f} degrees")
print(f"Difference in latitude: {delta_lat:.2f} degrees")

<class 'float'>
<class 'float'>
Difference in longitude: 45.53 degrees
Difference in latitude: 71.76 degrees

Making Decisions¶

A script that only stores numbers is not very useful on its own — the value of code comes from what it decides to do with those numbers. Python, like every programming language, lets you branch your logic with if, elif ("else if"), and else, combined with comparison operators (<, >, <=, >=, ==, !=) and logical operators (and, or, not).

Geographically, a latitude tells you which hemisphere a point falls in: positive values are north of the Equator, negative values are south of it. Longitude works the same way relative to the Prime Meridian. This gives us a first, tiny geographic classifier.

In [4]:

Copied!





# Classifying Rio de Janeiro by hemisphere using its latitude and longitude
if y1 > 0:
    ns = "Northern"
elif y1 < 0:
    ns = "Southern"
else:
    ns = "on the Equator"

if x1 > 0:
    ew = "Eastern"
elif x1 < 0:
    ew = "Western"
else:
    ew = "on the Prime Meridian"

print(f"Rio de Janeiro ({x1}, {y1}) lies in the {ns} and {ew} Hemisphere.")

# The same test written with a single boolean expression
is_southern_and_western = (y1 < 0) and (x1 < 0)
print("Southern and Western?", is_southern_and_western)
# Classifying Rio de Janeiro by hemisphere using its latitude and longitude
if y1 > 0:
    ns = "Northern"
elif y1 < 0:
    ns = "Southern"
else:
    ns = "on the Equator"

if x1 > 0:
    ew = "Eastern"
elif x1 < 0:
    ew = "Western"
else:
    ew = "on the Prime Meridian"

print(f"Rio de Janeiro ({x1}, {y1}) lies in the {ns} and {ew} Hemisphere.")

# The same test written with a single boolean expression
is_southern_and_western = (y1 < 0) and (x1 < 0)
print("Southern and Western?", is_southern_and_western)

Rio de Janeiro (-43.1729, -22.9068) lies in the Southern and Western Hemisphere.
Southern and Western? True

Collections¶

Two cities are easy to keep track of with four separate variables. Twenty cities are not. Python's built-in collections let you group related values together instead of inventing a new variable name for every one of them.

A list ([...]) holds an ordered, changeable sequence of items — useful for a column of coordinates.
A tuple ((...)) is like a list but immutable — a natural fit for a single, fixed coordinate pair (longitude, latitude).
A dictionary ({...}) maps a key to a value — perfect for looking up a city's coordinates by name.

These three structures are the pure-Python ancestors of the GeoSeries and GeoDataFrame objects you will meet in Chapter 7 — understanding them here makes those structures far less mysterious later.

In [5]:

Copied!





# A list of (longitude, latitude) tuples
coordinates = [
    (-43.1729, -22.9068),  # Rio de Janeiro
    (2.3522, 48.8566),     # Paris
    (-44.3028, -2.5297),   # São Luís, Brazil
]

# A dictionary mapping city names to their coordinate tuple
cities = {
    "Rio de Janeiro": (-43.1729, -22.9068),
    "Paris": (2.3522, 48.8566),
    "São Luís": (-44.3028, -2.5297),
    "Tokyo": (139.6917, 35.6895),
    "Cairo": (31.2357, 30.0444),
}

print(cities["São Luís"])       # Look up a single value by key
print(list(cities.keys()))      # All the keys
print(len(cities), "cities stored")
# A list of (longitude, latitude) tuples
coordinates = [
    (-43.1729, -22.9068),  # Rio de Janeiro
    (2.3522, 48.8566),     # Paris
    (-44.3028, -2.5297),   # São Luís, Brazil
]

# A dictionary mapping city names to their coordinate tuple
cities = {
    "Rio de Janeiro": (-43.1729, -22.9068),
    "Paris": (2.3522, 48.8566),
    "São Luís": (-44.3028, -2.5297),
    "Tokyo": (139.6917, 35.6895),
    "Cairo": (31.2357, 30.0444),
}

print(cities["São Luís"])       # Look up a single value by key
print(list(cities.keys()))      # All the keys
print(len(cities), "cities stored")

(-44.3028, -2.5297)
['Rio de Janeiro', 'Paris', 'São Luís', 'Tokyo', 'Cairo']
5 cities stored

Automation¶

Copy-pasting the hemisphere logic from before for all five cities in cities would work, but it does not scale — and it is exactly the kind of repetitive task computers exist to remove from you. Two tools handle this:

A for loop repeats a block of code once for every item in a collection.
A function (def) packages a block of code under a name, so you can reuse it without retyping it.

Combining the two turns the one-off classification from the previous section into something you can apply to any number of cities, forever, with a single line of code.

In [6]:

Copied!





def classify_hemisphere(lon, lat):
    """Return a short description of which hemisphere (lon, lat) falls in."""
    ns = "Northern" if lat > 0 else "Southern" if lat < 0 else "Equatorial"
    ew = "Eastern" if lon > 0 else "Western" if lon < 0 else "Prime Meridian"
    return f"{ns}/{ew}"

# Apply the function to every city with a for loop
for name, (lon, lat) in cities.items():
    print(f"{name:15s} -> {classify_hemisphere(lon, lat)}")
def classify_hemisphere(lon, lat):
    """Return a short description of which hemisphere (lon, lat) falls in."""
    ns = "Northern" if lat > 0 else "Southern" if lat < 0 else "Equatorial"
    ew = "Eastern" if lon > 0 else "Western" if lon < 0 else "Prime Meridian"
    return f"{ns}/{ew}"

# Apply the function to every city with a for loop
for name, (lon, lat) in cities.items():
    print(f"{name:15s} -> {classify_hemisphere(lon, lat)}")

Rio de Janeiro  -> Southern/Western
Paris           -> Northern/Eastern
São Luís        -> Southern/Western
Tokyo           -> Northern/Eastern
Cairo           -> Northern/Eastern

Building Models¶

You now have every ingredient a spatial model needs at its smallest scale: variables to hold state, conditionals to encode rules, collections to hold many observations, and functions to make the rules reusable. A model, at its simplest, is nothing more than a function (or a small set of them) that turns input data into a decision or a derived quantity, applied systematically across a dataset.

As a first, deliberately toy example, let's build a function that measures the straight-line ("as the crow flies") distance between two coordinate pairs, using the Pythagorean theorem, and use it to find which city in our dictionary is closest to Rio de Janeiro. This is not how you should measure real geographic distances — the Earth is curved, and Chapter 7 will introduce geopandas methods that account for that — but it is exactly the reasoning every one of the discrete spatial models in Part III of this book will build on: define a rule as a function, then apply it across many locations.

In [7]:

Copied!





import math

def flat_distance(coord1, coord2):
    """Straight-line distance between two (lon, lat) pairs, in degrees."""
    (x1, y1), (x2, y2) = coord1, coord2
    return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)

reference = cities["Rio de Janeiro"]

# Build a small report: distance from Rio to every other city
report = []
for name, coord in cities.items():
    if name == "Rio de Janeiro":
        continue
    report.append((name, flat_distance(reference, coord)))

# Sort the report by distance, closest first
report.sort(key=lambda item: item[1])

for name, dist in report:
    print(f"Rio de Janeiro -> {name:15s} {dist:8.2f} degrees")

closest = report[0][0]
print(f"\nClosest city to Rio de Janeiro (by this toy model): {closest}")
import math

def flat_distance(coord1, coord2):
    """Straight-line distance between two (lon, lat) pairs, in degrees."""
    (x1, y1), (x2, y2) = coord1, coord2
    return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)

reference = cities["Rio de Janeiro"]

# Build a small report: distance from Rio to every other city
report = []
for name, coord in cities.items():
    if name == "Rio de Janeiro":
        continue
    report.append((name, flat_distance(reference, coord)))

# Sort the report by distance, closest first
report.sort(key=lambda item: item[1])

for name, dist in report:
    print(f"Rio de Janeiro -> {name:15s} {dist:8.2f} degrees")

closest = report[0][0]
print(f"\nClosest city to Rio de Janeiro (by this toy model): {closest}")

Rio de Janeiro -> São Luís           20.41 degrees
Rio de Janeiro -> Paris              84.99 degrees
Rio de Janeiro -> Cairo              91.33 degrees
Rio de Janeiro -> Tokyo             192.02 degrees

Closest city to Rio de Janeiro (by this toy model): São Luís

Vectorized Operations¶

The for loop in the previous section is perfectly readable for five cities. For five million grid cells in a raster — which is exactly the kind of data Part II and Part III will hand you — a Python-level loop becomes a bottleneck, because each iteration carries its own interpreter overhead.

numpy solves this by letting you express the same calculation as a single operation over an entire array at once. Instead of looping city by city, you store every longitude in one array and every latitude in another, and let NumPy apply the arithmetic element-by-element in optimized C code. This is called vectorization, and it is the single most important performance idea in the rest of this book — Chapter 20 comes back to it in depth once real raster grids are on the table.

In [8]:

Copied!





# The same "distance from Rio" calculation, vectorized with NumPy
names = np.array([n for n in cities if n != "Rio de Janeiro"])
lons = np.array([cities[n][0] for n in names])
lats = np.array([cities[n][1] for n in names])

ref_lon, ref_lat = reference

# One expression computes the distance to every city simultaneously
distances = np.sqrt((lons - ref_lon) ** 2 + (lats - ref_lat) ** 2)

for name, dist in zip(names, distances):
    print(f"{name:15s} {dist:8.2f} degrees")

# No explicit Python loop was needed to do the arithmetic — NumPy did it internally
print("\nClosest city (vectorized):", names[np.argmin(distances)])
# The same "distance from Rio" calculation, vectorized with NumPy
names = np.array([n for n in cities if n != "Rio de Janeiro"])
lons = np.array([cities[n][0] for n in names])
lats = np.array([cities[n][1] for n in names])

ref_lon, ref_lat = reference

# One expression computes the distance to every city simultaneously
distances = np.sqrt((lons - ref_lon) ** 2 + (lats - ref_lat) ** 2)

for name, dist in zip(names, distances):
    print(f"{name:15s} {dist:8.2f} degrees")

# No explicit Python loop was needed to do the arithmetic — NumPy did it internally
print("\nClosest city (vectorized):", names[np.argmin(distances)])

Paris              84.99 degrees
São Luís           20.41 degrees
Tokyo             192.02 degrees
Cairo              91.33 degrees

Closest city (vectorized): São Luís

Exercises¶

Hemisphere counter. Using the cities dictionary and a for loop, count how many cities lie in the Southern Hemisphere and how many lie in the Northern Hemisphere.
Extend the toolbox. Add three more cities of your choice to cities (use real coordinates — an easy way to find them is to search "[city name] latitude longitude"). Re-run the distance report from Building Models and check whether the closest city to Rio de Janeiro changes.
A second rule. Write a function classify_distance(dist) that returns "near" if dist < 20, "moderate" if 20 <= dist < 60, and "far" otherwise. Apply it to the report list from Building Models using a for loop.
Vectorize your own rule. Rewrite Exercise 3 using NumPy arrays and a boolean mask (np.where or direct comparisons) instead of a loop.

In [9]:

Copied!

# Your code here
# Your code here

Summary¶

Key concepts introduced¶

Variables and data types (float, str, bool) as the basic containers for geographic values
Comparison and logical operators, and if / elif / else for encoding geographic rules
Lists, tuples, and dictionaries for grouping many observations together
for loops and functions (def) for automating repeated work
A first rule-based "model" built purely from these ingredients
Vectorization with NumPy as the scalable alternative to a Python-level loop

The rest of Part I builds directly on this toolbox: Chapter 3 replaces the hand-written cities dictionary with a pandas.DataFrame, and the collections, conditionals, and functions introduced here reappear, in vectorized form, throughout Part II and Part III.