Module Python
Module Python
Imagine you have a very literal-minded assistant who can perform tasks at incredible
speed, but only if you give him instructions in a specific, unambiguous language. This
assistant has no common sense. If you say “draw a map of my city,” he will stare at you
blankly. But if you say:
1. Create a blank canvas 800 pixels wide and 600 pixels tall.
2. For each road in this list: draw a line from its start coordinate to its end coordinate,
with color black and width 2.
3. For each park in this list: draw a green polygon defined by its boundary points.
Programming is the art of writing such step-by-step instructions for a computer. The
instructions are written in a programming language—a formal language with strict rules
(syntax) that the computer can translate into actions. The written instructions are
called code. Running those instructions is called executing or running the program.
A program is not magic; it’s a recipe. And like a recipe, if you miss a step or write something
ambiguous, the result might be inedible—or the computer might refuse to run it at all and
show you an error.
You might ask: “If I want to make a map, why not just open QGIS and click around?” That is
a valid approach for one-off tasks. But imagine you need to:
• Convert coordinates for millions of GPS points and check each against 500 different
boundary files.
• Automatically generate custom maps for every city in a country, each with the same
style, and update them weekly.
Doing that by hand would take years. Programming allows you to teach the computer to do
repetitive, complex, and large-scale tasks exactly once, and then repeat them flawlessly
forever, with minor changes only when needed. It is the difference between digging a canal
with a spoon and using an excavator. The excavator requires more initial training, but once
you know how, you can reshape the world.
So, programming empowers you to automate, to scale, and to work with data that is too
vast for any human to inspect manually.
Python is a programming language created by Guido van Rossum and first released in 1991.
Its guiding principle is readability: code should look almost like English, so that it is easy to
write and, more importantly, easy to read later by others or by your future self. Python
emphasizes a clean, uncluttered syntax—no excessive punctuation like { } ; found in many
other languages. Instead, it uses indentation (spaces at the beginning of a line) to define
blocks of code, which forces the programmer to write neatly.
Python is also interpreted, meaning you can write one line and immediately run it to see
the result, without a separate “compile” step. This interactive feedback makes learning
faster and experimentation easy.
• Case-sensitive: House and house are two different names. “Beijing” and “beijing”
are not the same.
• Dynamically typed: You don’t have to declare whether a variable will hold a number
or text; the computer figures it out when you assign it. This makes initial coding
faster but demands care.
• Whitespace matters: Indentation is not just for looks; it defines structure. We will
learn exactly how to use it.
• Free and open-source: Python costs nothing and runs on Windows, macOS, and
Linux.
3. Why Python for Geospatial? The Ecosystem
Geospatial work involves handling coordinates, maps, satellite imagery, elevation models,
and spatial relationships. Python has become the lingua franca of geospatial programming
because of a convergence of powerful, open-source libraries that build upon decades of
compiled C/C++ libraries (GDAL, GEOS, PROJ). You don't need to know those low-level
libraries; Python wraps them nicely.
• Shapely – Create and manipulate points, lines, and polygons. Check if a point is
inside a park.
• GeoPandas – Work with vector data (shapefiles, GeoJSON) as if they were Excel
tables with a geometry column. Filter, merge, and spatial-join with simple
commands.
• Rasterio – Read, write, and process raster data (GeoTIFF, satellite imagery) as arrays
of numbers.
These tools are used by government agencies, environmental scientists, urban planners,
disaster response teams, and tech companies. Learning Python for geospatial means you
can join this community and solve real-world problems that affect billions of people.
But before we can fly with these libraries, we must learn to walk with pure Python. The
fundamentals—variables, loops, functions—are the foundation upon which everything else
rests. This module will take you there step by step.
Think of Python as your kitchen. You have basic ingredients (numbers, text) and tools
(knives, pans). Programming is writing a recipe (the script). A library is like a pre-made
sauce or a specialized appliance you can bring into your kitchen to avoid making everything
from scratch. But you still need to know how to chop an onion (variables, loops) before you
can cook a feast. In this module, we will first master the knife skills, then we'll bring in the
blender and the sous-vide machine.
We will now set up your environment. I will assume you are using a Windows PC, as it's the
most common beginner platform, but I'll provide macOS notes where necessary. If you're
on Linux, you likely already know these steps, but I can assist separately.
2. The site will suggest the latest stable version (e.g., Python 3.12.x). Click the big
yellow button Download Python 3.12.x (or whatever number is current). Do not
download Python 2; that is obsolete.
• Tick the box at the bottom that says “Add Python 3.x to PATH”. This is extremely
important. It allows you to run Python from any command prompt.
• Then click “Install Now” (or “Customize installation” if you want to change
location, but the default is fine).
On macOS, Python is often pre-installed, but it’s best to install the latest official version
from the Python website. Open the downloaded .pkg file and follow the prompts. You’ll also
need to install the command line tools if prompted; we can cover that later if needed.
We need to open a terminal (also called command prompt). On Windows, press Windows
Key + R, type cmd, and press Enter. On macOS, open Terminal from Applications > Utilities.
bash
python --version
Press Enter. You should see something like Python 3.12.3. If you see a version number
starting with 3, congratulations! Python is installed.
Now type:
bash
python
And press Enter. This launches the Python interactive interpreter—a mode where you can
type Python code line by line and see results immediately. You’ll see a prompt like >>>.
Type:
python
print("Hello, World!")
Press Enter. The computer will print Hello, World! below your command. This is your first
line of code!
To exit the interpreter, type exit() and press Enter, or just close the terminal window.
You can write Python in Notepad, but a code editor provides syntax highlighting, error
detection, and a smoother experience. I recommend Visual Studio Code (VS Code)—free,
powerful, and used by millions.
2. Install it with the default options. On Windows, you might want to tick “Add to PATH”
and “Register Code as an editor for supported file types” for convenience.
3. Once installed, open VS Code. On the left sidebar, click the Extensions icon (looks
like four squares). Search for “Python” and install the official Microsoft Python
extension. This gives you IntelliSense (code completion), debugging, and more.
1. Open VS Code.
2. Click File > Open Folder... and choose a folder where you will keep all your learning
projects. For example, create a new folder on your Desktop
called python_geospatial_basics and select it.
3. In the Explorer sidebar (top icon looks like two documents), right-click inside the
folder area and choose New File. Name it lesson0_1.py. The .py extension tells VS
Code and the computer that this is a Python file.
python
6. To run this script, click the small play button (▶) at the top-right of the editor window
(if it’s there). Alternatively, open the terminal inside VS Code by clicking Terminal >
New Terminal from the menu. In the terminal panel at the bottom, type:
bash
python lesson0_1.py
If you get an error like python is not recognized, it means Python is not in PATH. Re-run the
Python installer, ensure “Add Python to PATH” is checked, or restart your computer. If
problems persist, tell me and I’ll help.
In the code print("I am learning Python for geospatial!"), you used several symbols. Let’s
decode them without deep theory, just for awareness:
• print : a built-in function – a piece of code that does a specific task (here, display
text on the screen). We’ll learn about functions in detail later.
• ( and ) : parentheses. They enclose the information you give to the function
(called arguments). print needs to know what to print, so you put the text inside.
• " " : double quotes. They mark the beginning and end of a text string (a sequence of
characters). Text must be enclosed in quotes so Python knows it’s not a variable
name or command.
• Whitespace : the space after print is optional but conventional; it improves
readability. The lack of space inside quotes is part of the text.
Notice there are no # symbols, no import, no def. We are building from the absolute
minimum.
You might think, “Printing a sentence is trivial; how will this help me map the world?” Every
complex geospatial script you will ever write is built from these basic blocks.
The print function will be your primary debugging tool for years—you will
insert print(variable_name) to peek inside your program’s brain. Understanding how to run
a script and see output is the gateway to everything else. Without this, you cannot test any
code. So we start here to establish a safe sandbox where you can make mistakes without
breaking anything.
1. File names: Use descriptive names, all lowercase, with underscores if needed
(my_first_script.py). Avoid spaces and special characters. The .py extension is
mandatory.
2. Save before running: Always save your file. Most editors remind you, but it’s a good
habit.
3. One statement per line: Python generally expects one instruction per line. Do not
write print("Hi") print("Bye") on the same line (it’s possible with a semicolon but
considered bad style).
4. Case matters: Print is not the same as print. Capital letters change the meaning.
5. Comments: Lines starting with # are ignored by Python. They are notes to yourself.
We’ll use them later.
6. Indentation: For now, every line of code should start at the very left edge, no
spaces. Later, we’ll add indentation when needed, but incorrectly indenting code
that shouldn’t be indented will cause errors.
9. Problem: Convert Requirements into Writing Syntax (End-of-Lesson Exercise)
Now it is your turn to write a small script, to practice the actions you just learned.
Problem Statement:
Create a new Python file named hello_geo_beginner.py. In this file, write code that does the
following:
2. Print a blank line (hint: you can call print() with nothing inside the parentheses).
3. Print the following lines, each on its own separate output line:
4. After the last print, include a comment on its own line that says # End of my first
script. (Remember, comments start with # and are ignored when running.)
Requirements:
• The output must match exactly what is asked (including the blank line). The date can
be anything you like.
• The script must run without errors when you execute it from the terminal (or VS
Code play button).
Why this exercise: You’ll practice creating a new file, writing multiple print statements,
seeing the sequential execution, and adding a comment. The blank line teaches
that print() with no arguments outputs an empty line. The multi-line output demonstrates
that each print advances to a new line.
Self-check: After running, your console output should look something like:
text
My name is Li Wei.
I am learning to code so I can automate maps and spatial analysis.
(Note: after the last line, you won't see the comment because comments are silent.)
Please write this script now. If you encounter any error, read the error message carefully; it
usually tells you the line number and the problem. If you cannot solve it, note exactly what
the error says, and I will help you interpret it. Once you have successfully run the script, we
will proceed to Lesson 1.1 where we will learn about variables and store your first
coordinate.
Take your time, and remember: every master was once a beginner who refused to give up.
Imagine you have many small boxes, each can hold exactly one thing: a number, a piece of
text, a coordinate. To find a box later, you stick a label on it and write a name. In Python,
these labeled boxes are called variables.
A variable is a name that refers to a value stored in the computer's memory. You create a
variable by writing a name, then an equals sign =, then the value you want to store. This is
called an assignment statement.
python
city_name = "Beijing"
latitude = 39.909
population = 21540000
Here, city_name, latitude, and population are variables. The equals sign does not mean
equality like in mathematics; it means “take the value on the right and store it under the
name on the left.” Think of it as an arrow pointing left: Beijing → city_name.
In geospatial analysis, you constantly work with entities that have attributes: a city has a
name, a location, a population; a sensor station has an ID, a latitude, a longitude, an
elevation. Variables allow you to give meaningful names to these pieces of data, so your
code reads like a description of the world, not a soup of numbers.
Without variables, every time you need the latitude of Beijing, you'd have to write 39.909. If
that value ever changes (say you get a more precise measurement), you would have to find
and replace every occurrence—a tedious and error-prone process. With a variable, you
change the value in one place, and the update propagates everywhere the variable is used.
Variables also let you write generic code. Instead of writing a script that only works for
Beijing, you write a script that uses variables like city_lat, city_lon, and then you can apply it
to any city just by changing the assignment. This is the beginning of reusable, scalable
programming.
I've seen beginners hard-code values everywhere, then struggle to update their scripts. I've
also seen them try to use a variable before it's assigned, which raises a NameError and
crashes the program. Learning to declare, assign, and use variables properly is like learning
to use nouns in a language. Without nouns, you cannot form sentences. Without variables,
you cannot form a program that models the real world.
Python has strict rules for variable names, and strong community conventions. Let's learn
both.
• Variable names can contain letters (a-z, A-Z), digits (0-9), and underscores (_).
• They cannot be Python keywords (like if, else, for, import, def, etc.). You don't need
to memorize the keyword list; your editor will highlight them, and you'll learn them
over time.
• If a name consists of multiple words, separate them with underscores (_). This style
is called snake_case and is the official Python convention.
Example: city_name, utm_zone, max_speed.
• Avoid single-letter names except in very specific contexts (like x, y for coordinates in
a small loop) or in mathematical formulas. Even then,
prefer lon, lat over l, L, x, y when meaning is spatial.
• Use descriptive names that tell you what the value is, not how it's
used. population is good; pop is ambiguous (pop what?). temperature_celsius is
better than temp.
You've seen _ inside variable names like city_name. The underscore is just a visual
separator; it has no special power by itself. It's a character that Python allows in names,
and by convention, we use it to join words because spaces are not allowed. So city name is
invalid, but city_name is the Pythonic way.
Let's write code to see variables in action. Create a new file lesson1_1.py.
python
city = "Beijing"
city = "Shanghai"
A variable can be updated. The old value is forgotten (unless another variable also refers to
it). This is why we call them variables—they can vary.
We can even assign a value of one type to a variable previously holding a different type, but
it's usually a bad idea because it confuses readers. For clarity, keep a variable's type
consistent.
Python allows you to assign multiple variables in a single line, which is very handy for
coordinates:
python
This assigns 116.397 to lon and 39.909 to lat simultaneously. It works as if you wrote:
python
lon = 116.397
lat = 39.909
This idiom is common when working with point data. It's clean and readable.
1. Always name variables descriptively. The name is the first documentation of your
code.
4. Avoid using built-in names (like print, list, str, type) as variable names—it can
override the original function and cause confusing errors.
5. Assign before use. If you try to print(unknown_var), Python will raise NameError:
name 'unknown_var' is not defined. Always make sure a variable exists (via
assignment) before you reference it.
7. Comment your intent if a variable's purpose isn't obvious from its name.
Now it is your turn to put variables into practice. Create a new file
called variables_practice.py.
Problem Statement:
You are collecting data for three weather stations located in different parts of a country. For
each station, you need to store its name, latitude (in decimal degrees), longitude (in
decimal degrees), and elevation (in meters above sea level). You must then print a
formatted report.
1. Create six variables for station 1: st1_name, st1_lat, st1_lon, st1_elev. Assign them
appropriate values (you can use dummy data, e.g., "Alpha", 34.05, -118.25, 100).
2. For station 2, use multiple assignment in one line to set st2_name, st2_lat, st2_lon,
st2_elev to values like "Beta", 40.71, -74.00, 10.
3. For station 3, assign variables st3_name, etc., with values like "Gamma", 51.51, -
0.13, 35.
6. For each station, print a single line that reads exactly like:
Use the variables you created; do not re-type the values. You must embed the variable
values inside the string. The simplest way is to pass multiple arguments
to print: print("Station", st1_name, ": lat=", st1_lat, ...). (Later we'll learn f-strings for more
elegant formatting, but for now, use this method.)
7. After printing the three station lines, print a footer line: --- End of Report ---.
Requirements:
• You must not type the actual station values inside the print statements; you must
reference the variable names.
• The output must look clean and professional, with one blank line before the header,
and each station on its own line.
• Include a comment at the top of your file explaining what the script does.
Self-check: After running, your console should show something like (with your chosen
values):
text
Please write this script. If you encounter any errors, especially NameError, check if you've
misspelled a variable name or forgotten to assign it. Remember: Python is case-sensitive,
so st1_lat is different from st1_Lat. Once you succeed, you will have cemented the skill of
naming and using variables—the bedrock of all future code.
Take your time, and when you are ready, we will move on to data types in detail, where we'll
explore numbers, strings, and the subtle dangers of floating-point representation. Proceed.
• How much memory it occupies (for now we won’t worry about memory, but it’s good
to know).
Python is dynamically typed, meaning you do not have to declare a type when creating a
variable. Python infers the type from the value you assign. However, Python is also strongly
typed—it will not automatically convert between incompatible types unless you explicitly
ask. For example, you cannot add the number 39.9 to the string "Beijing" without an explicit
conversion. This prevents many hidden bugs.
There are several built-in data types in Python. We will master the three most
fundamental: integers, floating-point numbers, and strings. We’ll also touch
on booleans for true/false logic.
2. Integers (int)
Integers are whole numbers, positive, negative, or zero. They have unlimited precision in
Python (unlike some languages where integers have a fixed size).
Examples in geospatial:
Operations on Integers
You can perform arithmetic: +, -, *, /, // (floor division), % (modulo,
remainder), ** (exponent).
Important: The division / always returns a float, even if the result is a whole number.
python
a = 10
b=3
print(a % b) # 1 (remainder)
Why floor division and modulo? In geospatial, you might need to split a grid of tiles: if you
have 1045 points and each tile holds 100, 1045 // 100 gives you the number of full tiles,
and 1045 % 100 gives the remainder on the last tile.
You can convert a string or float to an integer using int(), but it will truncate (not round) a
float and fail on non-numeric strings.
python
elev_str = "150"
Floats represent real numbers with a decimal point. In Python, they are implemented as
IEEE 754 double-precision (64-bit) numbers. They can represent numbers as large as
~1.8×10^308 and as small as ~2.2×10^-308 with about 15-17 decimal digits of precision.
Examples in geospatial:
Because floats are stored in binary, many decimal fractions cannot be represented exactly.
This leads to the famous gotcha:
python
In geospatial, if you compute a coordinate by summing many small offsets, tiny errors can
accumulate. Never test float equality with ==; instead, use a tolerance. We'll revisit this
when we write coordinate operations.
Scientific Notation
For very large or small numbers, you can use scientific notation: 1.5e6 means 1.5 × 10^6
(1,500,000). 2.54e-4 means 0.000254. This is useful for Earth-related constants like the
equatorial circumference in meters: EARTH_CIRCUM_M = 4.0075e7.
python
lat_str = "39.909"
4. Strings (str)
Strings are sequences of characters, enclosed in single quotes '...' or double quotes "...".
Both are equivalent; just be consistent. Use double quotes if your string contains an
apostrophe: "It's a map".
Examples in geospatial:
python
python
python
len("Beijing") # 7
• Accessing characters: You can index a string (like a list) to get a single character,
but we’ll cover indexing later. For now, just know it’s possible.
python
lat = 39.909
Without str(), "Latitude: " + lat would raise a TypeError because you can’t concatenate a
string and a float.
5. Booleans (bool)
Booleans have only two values: True or False (capitalization matters). They are the result of
comparisons and logical operations, and they are fundamental for decision-making.
Examples in geospatial:
== equal to 5 == 5 → True
Crucial note for floats: Because of precision issues, avoid == for floats; use tolerance-
based checks. For integers and strings, == is safe.
Logical Operators
python
is_valid = True
has_data = False
You can always discover the type of any value using the built-in type() function.
python
This is invaluable when debugging: if your script behaves unexpectedly, check the type of
your variables.
1. Use integers for counts, indices, and whole numbers. Don’t use floats for things
like number of points.
3. Use strings for text, labels, and identifiers that might have leading zeros. A
station ID like "001" is a string, not an integer, because 001 as an integer would
become 1.
4. Convert types explicitly using int(), float(), str(). Don’t rely on implicit conversion.
python
7. Use type() liberally while learning to understand what you’re working with.
python
center_lat = 39.9042
center_lon = 116.4074
# Simple Euclidean in degree space (not accurate, just for type demonstration)
# Build report
print("Site:", site_name)
Notice how we used every type naturally. The float precision issue isn’t visible here
because we’re not comparing floats for equality.
Problem Statement:
You are managing three weather stations. For each, you must store:
• Latitude (float)
• Longitude (float)
1. Assigns appropriate values for two stations using descriptive variable names in
snake_case, e.g., sta1_code, sta1_lat, etc.
2. Computes the simple planar distance in degrees between the two stations (using
the same Euclidean formula as above). Store in a variable delta_deg.
5. At the end, print a blank line, then print the type of the dist_km variable using type().
Self-check: Ensure the script runs without errors. Example output (with your chosen
values):
text
When you finish, you will have a solid grasp of the basic building blocks: variables,
assignment, and the fundamental data types. We'll then move into Lesson 1.3, where we'll
learn how to control the flow of your program—making decisions with if statements. Write
your script now.
A program without decisions is like a train on a single track—it can only go one way.
An if statement introduces a branch: a point where the program examines a condition
(something that is either True or False) and then decides which block of code to execute
next.
python
if condition:
statement1
statement2
The condition must be an expression that evaluates to a boolean (True/False). After the
condition, you write a colon :. Then, on the next line, you indent (usually 4 spaces) all the
lines that belong to the if block. When you return to the previous indentation level, you
signal the end of the block.
This indentation is not decoration—it is Python’s way of defining code structure. In many
other languages, you use curly braces { }; in Python, you use consistent whitespace. This
forces you to write visually clean code.
python
print("Hot day")
else:
If you have multiple exclusive conditions, chain if with elif (short for "else if"):
python
category = "Clear"
else:
category = "Overcast"
Only the first True block executes; the rest are skipped. The else is optional and catches
anything that didn’t match.
Why This Matters for Geospatial
By mastering if, you can teach your code to react intelligently to the data.
I’ve seen beginners write scripts without conditions that blindly process invalid data. For
example, they compute a distance for a coordinate (200, 100) that is clearly outside Earth’s
range. With a simple if check, you can stop the script or correct it early. An if statement is
your first line of defense against garbage-in, garbage-out.
Also, indentation errors are the most common frustration for beginners. Python will give
you an IndentationError if your spaces are inconsistent. Understanding exactly why
indentation matters will save you hours of staring at red error messages.
3.1 Basic if
python
altitude = 1200
print("Check complete.")
Output:
text
3.2 if-else
python
precipitation = 0.2
if precipitation > 0:
print("Rain detected")
else:
print("No rain")
python
magnitude = 5.7
if magnitude < 4:
alert = "Minor"
alert = "Moderate"
alert = "Strong"
else:
alert = "Major"
python
sensor_ok = True
reading = 45.3
if sensor_ok:
if reading > 100:
else:
print("Reading normal")
else:
print("Sensor malfunction")
You can nest as deep as needed, but deep nesting reduces readability. Often you can
combine conditions with and/or to flatten the structure.
python
lat = 35.0
lon = 135.0
hemisphere = "Northeast"
hemisphere = "Northwest"
# etc.
and requires both conditions to be True. or requires at least one. not flips a boolean.
python
python
if x > 0:
Every if body must be indented at least one space more than the if line. The standard is 4
spaces. Never mix tabs and spaces.
python
python
a = 0.1 + 0.2
python
python
if x < 0:
pass # Placeholder
else:
print("Non-negative")
2. Use parentheses for clarity when combining many and/or, but they are not
required if the order of operations (not, and, or) is clear.
3. Limit nesting to 2 levels. If deeper, consider breaking into functions (we’ll learn
later).
4. Include an else only if it handles a real case. Don’t add a useless else: pass.
5. Use elif for mutually exclusive conditions instead of a series of if statements that
might all be evaluated. elif stops after the first match.
6. Always indent 4 spaces. Configure your editor to insert spaces when you press Tab
(VS Code does this by default).
Let’s apply conditionals to a common task: checking whether a point lies within a
rectangular region (bounding box). This is a simplification of spatial filtering.
Create bounding_box_check.py:
python
# Define bounding box for study area (min_lon, min_lat, max_lon, max_lat)
min_lon = 116.2
max_lon = 116.6
min_lat = 39.8
max_lat = 40.0
# A point to test
point_lon = 116.4074
point_lat = 39.9042
if min_lon <= point_lon <= max_lon and min_lat <= point_lat <= max_lat:
else:
Python allows chained comparisons like min_lon <= point_lon <= max_lon, which is clean
and readable.
Problem Statement:
Write a script that categorizes weather conditions and issues an alert if dangerous.
1. Assign test values to the three variables. (Try different values to test all branches.)
o Else if wind_speed > 15 or precipitation > 1.0, set alert = "YELLOW - Caution"
4. Print a report:
text
Use your variables in the print statement (convert numbers to strings or use multiple
arguments).
Requirements:
• Test your script with at least two different sets of values to ensure all branches work.
Why this exercise: You practice multi-condition logic, nesting, and the critical safety that
conditionals bring to geospatial scripts. Weather categorization directly mirrors how you’d
classify land cover, hazard zones, or any spatial feature based on attributes.
Please write the script now. Debug any IndentationError or SyntaxError by checking your
colons and consistent spaces. Once you are satisfied, we will proceed to loops—the
engine that lets you repeat actions over thousands of features.
A loop is a block of code that runs multiple times. Python provides two primary loop
constructs:
• for loop: Iterates over a sequence (like a list, a string, or a range of numbers) and
executes the indented block once for each element in the sequence. You know in
advance (or can determine) how many times the loop will run.
• while loop: Repeats a block as long as a condition remains True. It may run zero
times, a fixed number of times, or indefinitely if the condition never becomes False.
In geospatial work, for loops are ubiquitous: loop over every point in a GPS track to
calculate speed; loop over every polygon in a land-use layer to compute area; loop over
every pixel in a raster window to apply a filter. while loops are less common for simple
iteration but are useful when reading data until a sentinel value is found, or when
implementing iterative algorithms that converge (like spatial interpolation).
The key idea: Do not repeat yourself (DRY). If you find yourself copying and pasting the
same code with minor changes (e.g., for station1, station2, station3), you should be using a
loop.
Without loops, you would be forced to process geodata manually or write absurdly long
scripts. Imagine checking if each of 10,000 points is inside a protected area. A loop lets you
write the check once and feed it the list of points. Loops are the difference between a
hand-drawn map and a computer-generated one. They unlock the ability to handle
datasets that are larger than your patience.
python
# indented block
statement1
statement2
• The indented block is the body of the loop; it runs once per element.
python
print(f"{temp}°C is HIGH")
else:
print(f"{temp}°C is normal")
This will print a line for each temperature in the list. Notice how temp is just a variable
name we choose; it could be t, reading, value. Choose descriptive names.
Often you need to loop a specific number of times, or over indices. range(start, stop,
step) generates a sequence of integers.
python
# Loop 5 times (i = 0, 1, 2, 3, 4)
for i in range(5):
print(f"Iteration {i}")
range(5) starts at 0 and stops before 5. You can provide a start and step:
python
python
coord_str = "116.4074,39.9042"
for ch in coord_str:
print(ch)
This prints each character (including comma and dot) on a new line.
When you have parallel lists (e.g., lons and lats), use zip to iterate them together:
python
zip pairs elements position-wise and stops when the shortest list ends.
If you need both the index (position) and the value while looping, use enumerate:
python
stations = ["Alpha", "Beta", "Gamma"]
4.1 Syntax
python
while condition:
# body
The loop checks the condition before each iteration. If it's True, runs the body; then checks
again. If the condition never becomes False, you get an infinite loop (press Ctrl+C to stop).
python
countdown = 5
print("Liftoff!")
You might use a while loop to read lines from a file until an empty line is found, or to run a
spatial interpolation until the change between iterations is below a threshold. For now,
we’ll focus on for loops as the primary tool.
5. Loop Control: break, continue, pass
• break: Immediately exits the innermost loop. Useful for early termination when a
condition is met.
python
if val < 0:
break
print(val)
• continue: Skips the rest of the current iteration and moves to the next element.
python
if val < 0:
print(val)
# Prints 1, 5, 7, 2
• pass: Does nothing; it's a placeholder. Can be used when a loop body is
syntactically required but you haven't written the logic yet.
python
6. Common Pitfalls
python
for i in range(3):
pass
• Infinite while loop due to a condition that never becomes False. Always ensure the
loop variable is updated inside the body.
• Modifying a list while iterating over it can produce unexpected results. Don't
remove items from a list while looping over the same list. Instead, create a new list.
1. Use for loops for iterating over a known collection or range. This is the default
choice for most spatial data processing.
2. Use while loops when the number of iterations is unknown and depends on a
condition computed inside the loop.
3. Choose descriptive loop variable names: for station in stations:, not for x in
stations:.
4. Keep the loop body short and focused. If it grows large, consider extracting parts
into functions (we’ll learn soon).
5. Avoid deeply nested loops (e.g., a for inside another for inside another for). Deep
nesting is hard to read and can be slow. Often, you can flatten logic using zip or list
comprehensions (advanced).
6. When using range(len(...)) to loop over indices, think twice. Usually, direct
iteration for item in list is clearer. Use enumerate if you need the index.
8. Use break sparingly—it should be obvious why the loop is terminating early.
8. Geospatial Application: Generating a Grid of Points
Let’s apply for loops to create a regular grid of sampling points, a common task in field
survey design.
Create grid_generator.py:
python
min_lon = 116.0
max_lon = 116.2
min_lat = 39.8
max_lat = 40.0
step = 0.05
# Generate longitudes and latitudes using range and a while-like approach with
multiplication
for i in range(lon_count):
for j in range(lat_count):
print(f"({lon:.4f}, {lat:.4f})")
This uses a nested loop: for each longitude step, we loop over all latitude steps. The result
is a regular set of points covering the rectangle. Notice how we converted floating-point
range to integer counts to avoid floating-point accumulation in loop counters.
Problem Statement:
You have a list of weather station codes and parallel lists of their latest temperature
readings (°C) and humidity percentages.
python
Write a script that uses a loop to process each station and produce a formatted report.
Additionally, use a while loop to simulate a simple iterative correction.
Requirements:
1. Use a for loop with zip to iterate over stations, temps, humidity simultaneously.
Self-check:
The station report should show 5 lines with appropriate statuses based on the data.
The while loop output might be:
text
Calibrating sensor...
(Note: 100.5 is still > 100.0? Actually 100.5 > 100.0 is True, so it would subtract again to
get 99.5. You'll want to adjust your condition to stop exactly when <= target. Use raw_value
> target. Starting 102.5: subtract 1.0 -> 101.5 (print), still > 100 -> subtract -> 100.5 (print),
still > 100? 100.5 > 100.0 is True, so subtract -> 99.5 (print), loop stops because 99.5 <=
100. Final printed calibrated value = 99.5. That's correct. The output would be three prints,
not two. So adjust expectation: 101.5, 100.5, 99.5 then final "Calibrated value: 99.5".)
Write the script now. Experiment by changing the lists to see different outputs. When you're
ready, we will move on to organizing data with more complex containers like lists, tuples,
and dictionaries—the backbone of geospatial data structures.
A variable holding a single number or string is useful, but the world is made of collections.
A GPS track is not one point but many. A land parcel has an owner name, a soil type, an
area, and a polygon geometry. To represent these in code, we need data structures that
can hold multiple values and, optionally, label them.
Python provides three built-in container types that you will use every day:
• List (list): An ordered, mutable (changeable) sequence of items. Items are accessed
by integer index (position). Use it when you have a collection of similar items where
order matters, and you might need to add, remove, or change elements.
Understanding when to use each is crucial for writing clean, efficient, and bug-free
geospatial Python.
I’ve seen beginners store coordinates as two separate lists (lons = [], lats = []) and then
struggle to keep them synchronized. A simple zip helps, but the root problem is poor data
modeling. A better approach is a list of tuples: points = [(lon1, lat1), (lon2, lat2)]. Each tuple
binds the two numbers into a single concept—a point. If you then need to attach attributes,
a dictionary per point or a GeoJSON-like structure is natural. Choosing the right container
makes your code reflect the structure of the problem, reducing errors and improving
readability.
3. Lists in Depth
A list is created using square brackets [], with items separated by commas.
python
empty = []
Indexing starts at 0 for the first element. Negative indices count from the end.
python
print(codes[0]) # WX01
print(codes[2]) # WX03
You can extract a sublist using [start:stop:step]. start is inclusive, stop is exclusive.
python
nums = [0, 1, 2, 3, 4, 5]
Slicing is extremely useful for extracting windows from raster data (once we have arrays) or
splitting a trajectory.
python
• [Link]() sorts in place (modifies list). sorted(lst) returns a new sorted list.
We already did this. You can loop directly over elements, or use enumerate for index+value.
You can create a list of lists, e.g., representing a 3x3 grid of elevation values.
python
elevation_grid = [
4. Tuples in Depth
python
color = (255, 0, 0)
empty = ()
4.2 Immutability
python
point[0] = 117.0 # TypeError: 'tuple' object does not support item assignment
This makes tuples safe to use as dictionary keys (see later) and conveys that the sequence
is a fixed unit. If you need a modifiable coordinate, use a list [lon, lat]; but in geospatial,
points are often better as tuples because a coordinate’s identity shouldn’t change.
python
This is how we elegantly swap values or return multiple values from a function.
• Use a tuple for fixed collections where position carries meaning: (lon, lat), (r, g,
b), (min_x, min_y, max_x, max_y).
• Use a list for variable-length collections of similar items: a list of GPS points, a list
of station names.
5. Dictionaries in Depth
Dictionaries are written with curly braces {} containing key-value pairs separated by colons.
python
city = {
"name": "Beijing",
"lat": 39.909,
"lon": 116.397,
"population": 21540000,
"is_capital": True
empty_dict = {}
python
print(city["name"]) # "Beijing"
print(city["population"]) # 21540000
If the key doesn’t exist, you get a KeyError. To avoid this, use get() with a default:
python
python
python
del city["is_capital"]
• [Link]() – returns a view of (key, value) tuples. Very handy for looping.
python
print(f"{key}: {value}")
You can combine dictionaries and lists to model geospatial data similar to GeoJSON:
python
feature = {
"type": "Feature",
"geometry": {
"type": "Point",
},
"properties": {
"visitors_per_year": 16000000
}
print(feature["geometry"]["coordinates"]) # (116.4074, 39.9042)
This structure is exactly how real geospatial data is represented in Python before we use
GeoPandas. Mastering nested dicts and lists gives you complete control.
1. Choose list for ordered, mutable collections of similar items. Use [].
6. Prefer direct iteration over lists rather than indexing with range(len(...)) unless you
need the index.
7. Use [Link](key, default) to safely access dictionary values when the key might be
missing.
8. Never modify a list while iterating over it with a for loop; create a new list or iterate
over a copy.
Let's combine these containers to represent three weather stations and print a summary.
Create containers_geo_demo.py:
python
stations = [
"code": "WX01",
"location": (116.38, 39.92), # tuple (lon, lat)
"temp": 22.5,
"humidity": 55,
"active": True
},
"code": "WX02",
"temp": 25.0,
"humidity": 60,
"active": False
},
"code": "WX03",
"temp": 19.0,
"humidity": 72,
"active": True
code = station["code"]
temp = station["temp"]
hum = station["humidity"]
This script shows how a list of dictionaries (or GeoJSON features) naturally models a vector
dataset. We’ll build on this structure heavily.
Problem Statement:
You have recorded a simple GPS track as a list of tuples (lon, lat). You also have attribute
dictionaries for the start and end points.
python
track = [
(116.395, 39.910),
(116.400, 39.912),
(116.405, 39.913),
(116.410, 39.915),
(116.415, 39.916)
2. Prints the first and last points using indexing (negative index for last).
3. Loops over the track points using enumerate and prints each point with its index (1-
based) and coordinates formatted to 3 decimal places.
4. Computes the total number of segments (segments = len(track) - 1) and prints it.
6. Prints the metadata dictionary (just print the whole dict) and then separately prints
the "length_deg" value formatted to 4 decimal places.
Requirements:
• Use proper container types: track as list of tuples, info as dicts, metadata as dict.
• For the distance loop, use indexing or zip to iterate over consecutive points. I
recommend zip(track[:-1], track[1:]) which pairs each point with the next one.
• Use [Link] for square root (import math). The distance formula: sqrt((lon2-
lon1)**2 + (lat2-lat1)**2).
Why this exercise: You will practice list indexing, tuple unpacking, dictionary construction,
nested containers, and a common geospatial pattern—computing segment lengths along a
polyline. This directly prepares you for working with Shapely LineStrings later.
Self-check: With the given track, the approximate length should be about 0.02 degrees
(roughly 2.2 km depending on the exact coordinates). Print the metadata dictionary to
verify its structure.
Write the script now. As always, run it, fix any errors, and observe the output. Next lesson,
we will learn about functions—the final piece that allows you to wrap logic into reusable
blocks.
A function is like a recipe in a cookbook. It has a name (e.g., "boil an egg"), it expects
certain ingredients (water, eggs), it performs a sequence of steps, and it produces a result
(a boiled egg). In Python, you define a function using the def keyword, give it a name, and
inside the indented body you write the instructions. Once defined, you can call the
function by its name, providing the required inputs (arguments), and it will execute its code
and possibly return a value back to you.
python
# Indented body
return result
• Parentheses () hold zero or more parameters—variables that will receive the input
values when the function is called.
• return (optional) sends a value back to the caller. If no return is written, the function
returns None (a special Python object meaning "nothing").
When you call the function, you provide arguments (the actual values) inside parentheses:
python
sum_value = function_name(5, 3)
Here 5 and 3 are arguments passed to parameter1 and parameter2. The function runs,
returns the sum, and it gets assigned to sum_value.
I’ve seen beginners write scripts that are 500 lines long, with the same distance formula
repeated 20 times. When they discover a bug (e.g., they forgot to convert degrees to
radians), they must hunt down every occurrence, and they inevitably miss one. By using a
function, you fix the bug once, and every call uses the corrected code. Functions also
reduce cognitive load: once you’ve written haversine_distance and trust it, you no longer
need to think about the math; you just use it. This is the art of abstraction—hiding complex
details behind a simple interface.
python
def print_greeting():
To call it:
python
Parameters are placeholders. You can name them anything, but they should be descriptive.
python
python
def celsius_to_fahrenheit(celsius):
return fahrenheit
Usage:
python
temp_f = celsius_to_fahrenheit(25)
print(temp_f) # 77.0
The return statement immediately exits the function. Any code after return (in the same
block) is not executed.
python
Usage:
python
You can make some parameters optional by providing a default value. If the caller omits the
argument, the default is used.
python
print(f"{greeting}, {name}!")
Usage:
python
When calling a function, you can specify arguments by name, which makes the code more
readable and allows you to skip some default parameters.
python
describe_point(lat=39.9, lon=116.4)
Keyword arguments must come after positional arguments if you mix them.
Immediately after the def line, a triple-quoted string describes what the function does, its
parameters, and what it returns. This is not just a comment; it becomes part of the
function’s documentation accessible via help().
python
Returns:
"""
return a + b
Then in the interactive interpreter, help(add) displays this text. This is a professional
standard.
Variables defined inside a function are local to that function. They exist only while the
function is executing and cannot be accessed from outside. Variables defined outside any
function are global; they can be read inside a function, but if you try to assign to them,
Python treats it as a new local variable unless you explicitly declare global var_name (not
recommended for beginners). This isolation prevents accidental interference.
python
x = 10 # global
def modify():
print("Inside:", x)
modify()
print("Outside:", x) # still 10
If you need a value from outside, pass it as an argument. If you need a result, return it.
4. Geospatial Examples: Functions You’ll Actually Write
Even though not accurate for long distances, it’s useful for quick checks.
python
import math
dx = lon2 - lon1
dy = lat2 - lat1
Call it:
python
python
return (min_lon <= lon <= max_lon) and (min_lat <= lat <= max_lat)
Usage:
python
python
c = 2 * [Link]([Link](a))
return R * c
3. One function, one task – if a function does two things, split it.
4. Keep it short – ideally a function should fit on your screen (20-30 lines). If longer,
break into smaller helper functions.
5. Use docstrings – always, even for small functions. It saves you time later.
6. Limit parameters – more than 4-5 parameters makes the function hard to use.
Consider grouping related parameters into a tuple or dictionary.
7. Avoid side effects – a function should either compute and return something, or do
something (like print to console), but not both unless clearly intended. Functions
that return values are more testable and composable.
8. Don’t use global variables inside functions unless absolutely necessary. Prefer
passing arguments and returning results.
9. Use return to exit early if a condition isn’t met (e.g., invalid input).
Create a new file geo_functions.py. You will define a set of utility functions for basic
geospatial calculations and then use them in a small script.
Requirements:
1. Write a function midpoint(lon1, lat1, lon2, lat2) that returns the midpoint
coordinates (lon, lat) of two points. The midpoint is simply the average of longitudes
and average of latitudes.
o Docstring.
o Compute and print the Haversine distance between them (reuse from
lesson).
o Use a check: if the distance (Haversine) is greater than 1 km, print "Points are
far apart", else "Close together".
Why this exercise: You will practice defining multiple functions with different signatures,
using imports, writing docstrings, and calling them in a script. The geospatial context
makes the functions meaningful.
Important: In the if __name__ == '__main__': block, we finally use this pattern that you’ve
seen. Explain: if __name__ == '__main__': means “execute this block only if the script is run
directly, not imported as a module”. This allows you to reuse the functions in other scripts
without running the test code. We’ll discuss it more in future lessons, but for now, include it
exactly as shown.
text
Write the script now. Debug any errors by reading the traceback carefully. Once your
functions are working, you’ve reached a major milestone. After this, we will start working
with Python’s built-in modules and eventually external geospatial libraries. Proceed.
LESSON 1.7 — Working with Files: Reading, Writing, and Parsing Simple CSV
1. Concept & Theory: Persistence and Data Exchange
A file is a sequence of bytes stored on your disk under a name. Python treats files
as streams of data that you can read from or write to. To work with a file, you must:
1. Open it, specifying a path and a mode ('r' for reading, 'w' for writing, 'a' for
appending).
The most common text file formats for geospatial data are:
• CSV (Comma-Separated Values): Simple tables where each line is a row and
columns are separated by commas (or other delimiters). Example: GPS tracks,
weather data.
• GeoJSON: A text format for geographic features (we’ll parse manually later, then use
libraries).
Every serious geospatial project starts with data ingestion. You will read thousands of
coordinates from a CSV, process them, and export results to a new file. Without file
handling, your scripts cannot interact with the data ecosystem—they become isolated
toys. By learning file I/O now, you prepare to load real-world datasets in subsequent
lessons.
Variables live only as long as the program runs; they vanish when it ends. Files persist. If
you collect GPS points during fieldwork, you must save them to a file. If you want to share
your analysis results, you write them to a file. File handling is the gateway between your
transient code and the permanent, shareable world of data.
python
content = [Link]()
print(content)
• The path '[Link]' is relative to where your script is run (or absolute
like 'C:/data/[Link]'). Use forward slashes / or double backslashes \\ on
Windows.
Best for large files: iterate over the file object line by line (memory efficient).
python
If a file is not in the same folder as your script, you must provide the full path. You can use
the os module to build paths safely, but for now, manual strings work. Forward slashes
work on all platforms: 'C:/Users/Name/Desktop/[Link]'.
python
If you need to write multiple lines, you can also use [Link](lines) (with newlines
included).
Warning: Opening a file in 'w' mode destroys any existing content instantly. Double-check
your paths.
CSV format: first row is often a header (column names), subsequent rows are data.
text
code,lat,lon,temp
WX01,39.909,116.397,22.5
WX02,31.230,121.473,25.0
We can parse it manually to understand the process before using the csv module (which
we’ll learn with import in a future lesson).
python
lines = [Link]()
print("Columns:", header)
for line in lines[1:]: # skip header
parts = [Link]().split(',')
code = parts[0]
lat = float(parts[1])
lon = float(parts[2])
temp = float(parts[3])
This manual approach fails if fields contain commas (e.g., quoted text), which is why
the csv module exists, but it’s fine for simple numeric data.
python
data_rows = [
[Link]("code,lat,lon,temp\n") # header
line = f"{row[0]},{row[1]},{row[2]},{row[3]}\n"
[Link](line)
7. “How to Write” Rules for File Handling
4. When reading a file, always process line by line if the file could be large (e.g.,
millions of GPS points).
6. Handle missing files – for now, if the file doesn’t exist, Python
raises FileNotFoundError. In production, you’d wrap in try/except (later lesson).
7. Never hard-code sensitive file paths in scripts you share; use relative paths or
command-line arguments.
8. Use strip() on each line to remove trailing newline and spaces before parsing.
8. Geospatial Example: Reading GPS Track CSV and Writing Filtered Results
Suppose you have a file gps_track.csv with columns: lon,lat,elevation,time. You want to
read it, filter out points with elevation below 50 meters, and write a new CSV with the
remaining points.
python
valid_points = []
for line in f:
parts = [Link]().split(',')
lon = float(parts[0])
lat = float(parts[1])
elev = float(parts[2])
time = parts[3]
[Link]("lon,lat,elevation,time\n")
for pt in valid_points:
[Link](f"{pt[0]},{pt[1]},{pt[2]},{pt[3]}\n")
This script mirrors exactly what you’ll do later with GeoPandas but at the raw Python level.
Create a new file file_processing.py. You will read a CSV of weather stations, perform a
simple calculation, and output a report file.
Problem Statement:
Assume you have a file named weather_data.csv with the following content (create this file
in the same folder as your script):
text
station,lat,lon,rainfall_mm
Alpha,34.05,-118.25,12.5
Beta,40.71,-74.00,8.2
Gamma,51.51,-0.13,15.0
Delta,48.85,2.35,9.7
Epsilon,35.68,139.76,18.3
1. Opens and reads the CSV file manually (no csv module).
o Else: "Wet"
6. After writing, prints "Report generated with X stations." (X = number of data rows).
Requirements:
• Convert latitude and longitude to floats (for potential future use), but write them
back as strings.
Why this exercise: You’ll practice the complete file I/O cycle (read → process → write) in a
geospatial context, reinforcing loops, conditionals, string parsing, and proper file
management.
text
Alpha,34.05,-118.25,12.5,Moderate
Beta,40.71,-74.0,8.2,Dry
...
Write the script now. In the next lesson, we will finally unlock the power of Python’s
standard library with import, starting with the math module you briefly saw, and then
the csv module to simplify exactly this kind of task. Proceed.
A module is simply a Python file (.py) that contains variables, functions, and classes that
can be reused. When you write import math, you are telling Python: “Load the [Link] file
(which is part of the standard library), execute its code, and make its contents available to
me under the name math.”
Think of a module as a toolbox. Your current script is your workbench. You bring a toolbox
(import), open it, and use its tools by referencing the toolbox name: [Link](25). This
keeps your workbench uncluttered and lets you switch toolboxes as needed.
Python’s standard library is a collection of modules that come with every Python
installation. It includes modules for mathematics, file management, internet access,
date/time handling, and much more. The math module provides trigonometric, logarithmic,
and constant functions. The csv module provides robust tools for reading and writing CSV
files (handling quoted fields, different delimiters, etc.). Later, third-party modules (installed
via pip) will give us geospatial superpowers.
• import module_name as alias – imports with a short name (e.g., import numpy as
np).
I could ask you to write a cosine function from scratch using Taylor series. That would be
educational but impractical. The math module’s cos() is written in C, optimized, and tested
by thousands. Using it saves time and prevents bugs. The csv module correctly handles
edge cases (fields with commas, quotes, line breaks) that our manual split(',') could never
handle. Importing is the act of standing on the shoulders of giants. In geospatial, you will
import geopandas, not rewrite shapefile parsers. This lesson is the gateway to that world.
python
import math
radius = 5.0
print(area)
Here [Link] is a constant from the module. [Link](x, y) raises x to power y. The module
name acts as a namespace, preventing collisions with your own variables (e.g., you can still
have a variable called pi without conflict).
python
import math as m
print([Link]([Link](30))) # 0.5
This is useful for modules with long names or when a short alias is conventional
(e.g., import pandas as pd).
python
angle_deg = 60
print(cos(radians(angle_deg))) # 0.5
Now cos and radians are directly available without the math. prefix. This can make code
cleaner, but overuse might cause name conflicts.
python
print(sin(pi/2)) # works, but where did sin and pi come from? Confusing.
This dumps all names from math into your global namespace, which can silently override
your variables or built-ins. It’s acceptable in rare interactive sessions, but never in
production scripts.
If you have a file geo_utils.py in the same directory, you can do:
python
import geo_utils
geo_utils.haversine_distance(...)
Or
python
The math module provides everything you need for coordinate computations on a sphere.
python
import math
initial_bearing = math.atan2(x, y)
This function uses math.atan2 to get the correct quadrant, then normalizes to 0–360°.
We parsed CSV manually in the previous lesson. Now let the csv module do the heavy
lifting correctly.
python
import csv
reader = [Link](f)
print("Columns:", header)
[Link] returns an iterator that yields lists of strings. It automatically handles quoted
fields, embedded commas, etc. The next(reader) gives the first row, then the loop
consumes the rest.
Even more readable: each row becomes a dictionary with column names as keys.
python
lat = float(row['lat'])
lon = float(row['lon'])
temp = float(row['temp'])
DictReader automatically uses the first line as fieldnames. No need for next().
python
data = [
writer = [Link](f)
[Link](header)
[Link](data)
Important: Always use newline='' when opening files for CSV writing; otherwise, you might
get extra blank lines on Windows.
python
rows = [
[Link]()
[Link](rows)
python
if __name__ == '__main__':
main()
Now understand why. Every Python file has a special variable __name__. When you run a
file directly (e.g., python [Link]), Python sets __name__ to the string '__main__' for that
file only. If that same file is imported as a module (import script), its __name__ is set to the
module’s name ('script'), not '__main__'.
Thus, the if block lets you write code that should only run when the file is executed as the
main program, not when it’s imported. This allows you to include test code or demo code in
a module without it firing when someone imports your functions.
Example:
geo_tools.py:
python
# Test code
• In another script, import geo_tools → only defines midpoint, does not print anything.
1. Put all imports at the top of your file, right after the module docstring.
2. Order imports: standard library first, then third-party, then local modules. Separate
groups with a blank line.
4. Use aliases sparingly and only when they are conventional (import numpy as
np, import pandas as pd).
6. Be explicit: prefer from math import cos, radians only when you use them
frequently and there’s no risk of name clash.
8. Always use the __name__ guard in files that can be both scripts and modules.
Let’s rewrite our track processing using the proper modules. Assume gps_track.csv has
columns lon,lat,time.
python
import csv
import math
def haversine_distance(lon1, lat1, lon2, lat2):
R = 6371.0
a = [Link](dlat/2)**2 + [Link](lat1)*[Link](lat2)*[Link](dlon/2)**2
return R * 2 * [Link]([Link](a))
points = []
reader = [Link](f)
writer = [Link](f)
for i, pt in enumerate(points):
if i == 0:
else:
if __name__ == '__main__':
process_track('gps_track.csv', 'track_with_distances.csv')
Create a new file station_distances.py. This will combine csv reading, math functions, and
custom functions.
Problem Statement:
You have a CSV file [Link] with columns: id, lat, lon. Example content:
text
id,lat,lon
A,34.05,-118.25
B,40.71,-74.00
C,51.51,-0.13
D,35.68,139.76
o Print the two stations that are farthest apart and the distance.
o Write the distance matrix to a CSV file [Link] with the same tabular
format (IDs as first row and first column).
Requirements:
• The distance_matrix function must be pure: take a list of station dicts, return a list of
lists.
• The output CSV must have the station IDs as headers and first column.
• Follow all rules for imports, file I/O, and function design.
Why this exercise: You will integrate csv, math, nested loops, functions, and file writing.
This is a miniature version of a spatial analysis pipeline that calculates a distance matrix—
common in logistics, ecology, and urban planning.
Self-check: With the sample data, farthest pair might be A–D or B–C depending on
coordinates. Verify with an online distance calculator that your Haversine implementation
is correct. The output CSV should open in a spreadsheet and show distances.
Write the script now. When you finish, we will be fully prepared to step into the geospatial
libraries—first with Shapely for geometry objects. Proceed.
Up to now, we represented a coordinate as a tuple (lon, lat) and a line as a list of such
tuples. That worked for simple tasks, but imagine needing to:
• Polygon – an area bounded by an exterior ring (and optionally interior rings for holes).
Each geometry object stores its coordinates and provides methods (functions attached to
the object) that operate on that geometry: area, length, distance, buffer, intersection, etc.
They also provide predicates that return booleans: contains, intersects, touches, within.
Shapely is not a GIS; it doesn’t handle coordinate reference systems or attributes—it purely
manipulates abstract 2D Cartesian geometry. Later, geopandas will combine Shapely
geometries with data tables and CRS information. For now, we focus purely on the shapes.
Shapely operates in a flat, two-dimensional plane. When you feed it geographic longitude
and latitude, it will treat them as if they were X and Y on a plane. That means area
calculations are in square degrees (meaningless) and distances are in degrees (not
constant). This is fine for operations like containment checks where exact metric area isn’t
required, but be aware: for accurate measurements, you must first project to a planar
coordinate system (we’ll learn that with pyproj). The concepts we learn here apply
universally.
Without Shapely, you would spend a year implementing and debugging algorithms that
have been refined over decades. With Shapely, you write one line: [Link](point).
It’s not just about speed; it’s about correctness. GEOS is used by PostGIS, QGIS, and
thousands of organizations worldwide. By using Shapely, you tap into that reliability.
Additionally, Shapely geometries are the native objects that geopandas stores in
its geometry column, so mastering Shapely is mandatory for vector data analysis.
bash
python
We also import the module itself for advanced operations later: import shapely.
The simplest geometry: a single coordinate pair. You can create a Point from a tuple or by
passing coordinates directly.
python
# From a tuple
p2 = Point(121.4737, 31.2304)
Shapely points can have 2 or 3 dimensions (Z). There is also a “measured” dimension (M)
but we’ll ignore it.
4.2 Inspecting a Point
python
p = Point(116.4074, 39.9042)
print(p.geom_type) # 'Point'
print(p.is_empty) # False
The coords attribute is a CoordinateSequence object; you can convert it to a list or tuple to
see the raw values.
• distance(other) – returns the minimum distance to another geometry (in the same
planar units).
python
p1 = Point(0, 0)
p2 = Point(3, 4)
5.1 Creation
python
line2 = LineString(points)
python
print(line.geom_type) # 'LineString'
print([Link]) # (minx, miny, maxx, maxy) -> (0.0, 0.0, 3.0, 4.0)
[Link] gives the planar length. If the coordinates are degrees, the length is in degrees.
You can access the boundary points: [Link][0] gives the start, [Link][-1] the end.
However, to get interior Point objects, you can iterate over the coords and build Points, or
use [Link] for advanced splitting.
• interpolate(distance) – returns a Point at a given distance along the line from the
start.
• project(point) – returns the distance along the line to the nearest point to the given
point.
python
midpoint = [Link]([Link] / 2)
A Polygon represents a filled area defined by an exterior ring (a closed LinearRing, which is
like a LineString but explicitly closed) and zero or more interior rings (holes).
The exterior ring must be closed (first and last point identical). Shapely will close it if you
don’t, but it’s good practice.
python
# Square: (0,0) -> (1,0) -> (1,1) -> (0,1) -> (0,0)
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
If you omit the final point, Shapely automatically closes it, but explicitly closing ensures
clarity.
Pass the exterior ring as the first argument and a list of interior rings as the second.
python
exterior = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]
hole = [(3, 3), (7, 3), (7, 7), (3, 7), (3, 3)] # closed
python
print([Link]) # 1.0 (in square units)
The centroid is the geometric center (not always inside the polygon for concave shapes).
For a guaranteed interior point, use polygon.representative_point().
python
p = Point(0.5, 0.5)
print([Link](p)) # True
1. Import only what you need: from [Link] import Point, LineString,
Polygon. This avoids namespace clutter.
2. Use tuples for coordinates – Shapely accepts lists or tuples, but tuples convey
immutability of the coordinate pair.
3. Close your polygons explicitly – though not strictly required, it prevents confusion.
4. Be aware of coordinate order – Shapely uses Cartesian (x, y) = (lon, lat) for
geographic. Never swap.
6. Geometry objects are immutable – once created, you cannot change coordinates.
Create a new geometry instead.
7. Use geom_type for debugging – printing a Point shows POINT (x y), which is helpful.
Let’s write a script shapely_intro.py that defines some points, lines, and polygons,
performs basic operations, and prints results.
python
bbox = Polygon([
(115, 30), (125, 30), (125, 40), (115, 40), (115, 30)
])
buffer_zone = [Link](0.5)
intersection = buffer_zone.intersection(bbox)
Run this script. You’ll get numeric outputs. Notice that all operations are in degree space,
which is distorted. In later lessons, we’ll project to meters.
Problem Statement:
• Polygon P: a triangle with vertices (2, 2), (2, 8), (8, 2) representing a protected zone.
3. Check if Point A is on the line L (use [Link](point) < 1e-9 as tolerance; there is
no on_line predicate, but distance works). Print the result.
6. Buffer Point A by a radius of 3. Compute the area of intersection between this buffer
and Polygon P.
7. Create a new polygon Q that is the union of Polygon P and a square polygon (4, 4) ->
(7, 4) -> (7, 7) -> (4, 7) -> (4, 4). Print the area of Q.
Requirements:
Why this exercise: You will practice creating multiple geometry types, using spatial
predicates, set operations (intersection, union), and buffer. These are the core operations
in any spatial analysis workflow.
Self-check:
• A is on L (distance 0) → True.
• Buffer A radius 3: circle. Intersection area with triangle? Visually, triangle from (2,2)
to (2,8) to (8,2). Buffer around (0,0) radius 3 touches triangle? The triangle's closest
point to origin is (2,2) distance sqrt(8) ≈ 2.828, so buffer 3 will overlap. The
intersection area can be computed. I'll check with Shapely later, but the student will
get some value.
• Union area Q: triangle P (area 18) + square (3x3=9) overlapping; union area will be
less than 27. The square (4,4)-(7,7) partially overlaps triangle? Triangle occupies x
from 2 to 8, y from 2 to 8. Square is fully inside the bounding box but does it intersect
triangle? The triangle's hypotenuse is x+y=10. For square points: (4,4) -> 8 < 10
inside triangle? Actually triangle interior: x+y >? For points (2,2) is min, interior is x >
2, y > 2? Wait, the triangle vertices: (2,2), (2,8), (8,2). The triangle is the set of points
with x >= 2, y >= 2, and x + y <= 10 (below the hypotenuse). So square (4,4) is x+y=8
≤10 inside, (7,4) x+y=11 >10 outside, (7,7) 14 outside, (4,7) 11 outside. So only part
of the square is inside. Union area = 18 (triangle) + area of square not overlapped.
Overlap area could be computed. But we just print whatever Shapely gives.
In computational geometry, the relationships between two planar shapes are defined by
the DE-9IM model (Dimensionally Extended 9-Intersection Model). This model classifies
the intersection of the interiors, boundaries, and exteriors of two geometries. Shapely
implements the most common predicates derived from DE-9IM, giving you an intuitive
vocabulary.
Geometries are topologically equal (same shape and Are two representations of the
equals
location) same city boundary identical?
intersects The geometries have at least one point in common Does a road pass through a par
The geometries share interior points but neither Do two forest stands partially
overlaps
contains the other (same dimension) overlap?
disjoint The geometries have no points in common Is a point far from any building?
Each predicate is a method of the first geometry object, taking the second as
argument: [Link](point). They return a boolean. Understanding these allows you
to filter, select, and analyze spatial features with precision.
Why This Matters: In GIS, attribute queries (like population > 100000) can be done with
Pandas. Spatial queries (like point inside polygon) require Shapely predicates. They are the
building blocks of spatial joins, overlay analysis, and map algebra.
Many beginners use intersects when they really need within. A point on the boundary of a
polygon is not within (because the interior of the point does not lie entirely in the interior of
the polygon), but it does intersect. Misusing predicates leads to off-by-one errors in feature
counts, misclassification, and silent logical flaws. Learning the exact semantics of each
predicate makes your code robust and defensible.
Furthermore, performance: when checking many points against a complex polygon, naive
predicate calls can be slow. Shapely provides prepared geometries – a one-time
transformation that creates an optimized spatial index inside the polygon, making
repeated contains or intersects calls orders of magnitude faster. We’ll learn this pattern as
well.
We’ll build a script predicates_demo.py to explore each predicate. First, we need some
sample geometries.
python
outer = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]
hole = [(4, 4), (6, 4), (6, 6), (4, 6), (4, 4)]
# Points
line_cross = LineString([(-1, 5), (5, 5), (11, 5)]) # passes through hole
poly_overlap = Polygon([(8, 8), (12, 8), (12, 12), (8, 12), (8, 8)])
3.1 equals
Two geometries are equal if they represent the exact same set of points (same shape and
location, regardless of vertex order or representation details). Use equals or ==.
Shapely’s == uses equals.
python
p1 = Point(0, 0)
p2 = Point(0.0, 0.0)
print([Link](p2)) # True
print([Link](poly2)) # True
I'll use the standard definition: [Link](B) means B is completely inside A, possibly
touching the boundary. For points, a point on the boundary is not considered contained. I’ll
test with actual Shapely: I recall that a point on a polygon boundary does not
satisfy contains, but does intersect. I'll explain accordingly.
So in code:
python
print([Link](interior_point)) # True
print([Link](corner_point)) # False
print([Link](far_point)) # False
3.3 intersects
[Link](B) → True if the geometries have at least one point in common. This is the most
general "touches or crosses or overlaps or contains". Even boundary contact counts.
python
print([Link](interior_point)) # True
print([Link](far_point)) # False
print([Link](line_cross)) # True
intersects is the safest predicate when you just want to know if they share any point. But be
careful: it doesn't distinguish between a point inside and a point on the edge.
3.4 touches
[Link](B) → True if the geometries share at least one boundary point but their interiors
do NOT intersect. For polygon/point: the point must be on the boundary of the polygon (and
not inside). For line/polygon: the line touches if it shares boundary but doesn’t go inside.
python
print([Link](border_point)) # True
print([Link](corner_point)) # True
print([Link](far_point)) # False
3.5 crosses
[Link](B) → True if the intersection is not empty and the dimension of the intersection is
less than the maximum dimension of the two geometries, and the geometries are not
subsets of each other. Typically: a line crosses a polygon if part of the line is inside and part
is outside.
python
3.6 overlaps
[Link](B) → True if the geometries are of the same dimension and they share interior
points but neither contains the other. For two polygons: they partially overlap.
python
print([Link](poly_overlap)) # True (poly_overlap overlaps region partially)
3.7 disjoint
python
print([Link](far_point)) # True
python
union = [Link](poly_overlap)
intersection = [Link](poly_overlap)
print([Link])
diff = [Link](poly_overlap)
print([Link])
These operations return new Shapely geometries, possibly GeometryCollection if the result
is heterogeneous (e.g., polygon with a hole and a separate polygon).
When you need to test many points (or small geometries) against a complex polygon,
creating a prepared version of the polygon pre-computes a spatial index. Then
repeated contains, covers, or intersects calls run much faster.
python
prepared_region = prep(region)
for pt in list_of_points:
if prepared_region.contains(pt):
# do something
The prepared object has only contains, covers, and intersects methods. Use it in loops over
large datasets.
1. Choose the right predicate – use within/contains when you need to know if a
feature is fully inside; use intersects for any contact.
Let’s apply predicates to categorize weather stations relative to a study area polygon.
python
study_area = Polygon([(100, 20), (120, 20), (120, 40), (100, 40), (100, 20)])
# Stations
stations = {
prep_area = prep(study_area)
if prep_area.contains(point):
status = "Inside"
elif study_area.touches(point):
status = "On boundary"
elif study_area.intersects(point):
else:
status = "Outside"
print(f"{name}: {status}")
Problem Statement:
You are given a polygon representing a nature reserve, and a set of features (points and
lines) that represent sensors and trails. You need to classify them and compute spatial
summary.
Reserve polygon (list of exterior coordinates): [(0,0), (10,0), (10,10), (0,10), (0,0)] (no holes).
Features:
• Sensor A: Point(2, 2)
2. For each feature, determines and prints its relationship to the reserve using the
most specific predicate (choose among within/contains for points, crosses for
lines, touches, etc.). Print a message like "Sensor A is inside the reserve."
5. Creates a 0.5-unit buffer around Sensor A. Compute and print the area of
intersection between this buffer and the reserve.
6. Uses prepared geometry for the reserve when testing the sensors in a loop.
Requirements:
Why this exercise: You’ll integrate all the day’s concepts—predicates, classification,
measurements, and preparation—into a practical spatial query tool. This is a microcosm of
a GIS “select by location” workflow.
Self-check:
• Trail 1: within? Does a line that lies completely inside a polygon satisfy within?
Yes, [Link](Polygon) returns True if the entire line is inside (interior or
boundary? the line's interior must be inside the interior of polygon? Actually, line
within polygon: True if the line’s interior and boundary are completely inside the
polygon’s interior or boundary. Trail 1 is entirely inside, not touching boundary
except at endpoints? (1,1) is interior, (9,9) interior. So within = True. So Trail 1 is
"within". So within count for lines? I will just classify lines separately: "Trail 1 is
completely inside", "Trail 2 crosses", "Trail 3 touches boundary". The problem asks to
print relationship. So output accordingly.
• Trail 2 crosses.
• Buffer around A radius 0.5: circle area = pi*0.25 ≈0.7854. Intersection with reserve:
since buffer is fully inside (center at 2,2, radius 0.5, min x=1.5, max x=2.5, still
inside), intersection area = buffer area ≈0.7854. Print to a few decimals.
This exercise will solidify your command. Write the script now. After this, we will move to
the cornerstone of vector data analysis: GeoPandas. Proceed.
Recall our earlier work: we stored station attributes in dictionaries and coordinates in
separate tuples. If we wanted to filter stations by temperature and then check which are
inside a polygon, we had to juggle two separate structures—the list of attributes and the list
of Shapely objects. This quickly becomes messy and error-prone.
• pandas DataFrame – A 2-dimensional, tabular data structure with rows and labeled
columns, capable of holding mixed data types (numbers, strings, booleans). It’s like
an Excel spreadsheet inside Python.
A GeoDataFrame is a pandas DataFrame with a geometry column. It inherits all the power
of pandas for data manipulation (filtering, grouping, aggregating, merging) and adds spatial
operations (spatial joins, intersection, buffering, distance calculations) that operate on the
geometry column.
Key concepts:
• Geometry Column: The active geometry column (you can have multiple, but one is
active). You can access it with [Link]. It’s a GeoSeries, which is a
pandas Series of Shapely geometries with CRS information.
Without GeoPandas, every spatial analysis would require custom loops over lists of
Shapely objects and manual attribute management. GeoPandas provides a unified,
high-performance framework that handles millions of features efficiently (using underlying
GEOS and pygeos/rtree for spatial indexing). It also seamlessly reads and writes dozens of
vector file formats (Shapefile, GeoJSON, GeoPackage, PostGIS, etc.). In professional
geospatial Python, GeoPandas is as fundamental as a hammer to a carpenter.
I’ve seen researchers write 500-line scripts to load, filter, and intersect two shapefiles,
when GeoPandas can do it in 10 lines. Not only is the GeoPandas code shorter, it’s also
more readable and less likely to contain bugs. By mastering GeoPandas, you can stop
worrying about the mechanics of data handling and focus on the spatial questions that
interest you. You also become part of a massive community where solutions are widely
shared.
GeoPandas depends on Shapely, Fiona, PyProj, and other packages. The recommended
installation via conda (as we set up in Lesson 0.1) ensures all dependencies are
compatible. If you haven’t installed GeoPandas yet:
bash
python
python
import pandas as pd
python
A GeoDataFrame can be created in several ways. The most instructive is to build one from a
dictionary or a list of dictionaries, where one column contains Shapely geometry objects.
python
# Build dictionary
data = {
"code": station_codes,
"temperature": temperatures,
print(gdf)
Output:
text
Notice:
• The CRS is set to EPSG:4326 (WGS84 latitude/longitude). Always set the CRS when
creating a GeoDataFrame; without it, spatial operations that depend on CRS
(like to_crs) will complain.
python
features = [
]
gdf = [Link](features, crs="EPSG:4326")
The geometry column is automatically detected if its name is exactly 'geometry'. If your
geometry column has a different name (e.g., 'location'), you can specify it with
the geometry parameter:
python
GeoPandas can read most vector formats using the read_file() function, which is powered
by Fiona.
python
gdf = gpd.read_file("path/to/[Link]")
python
gdf = gpd.read_file("[Link]
python
gdf.to_file("[Link]", driver="GeoJSON")
For Shapefile, the driver is automatically selected, but a directory will be created. Better to
use GeoPackage for modern projects.
python
# Combine conditions
Suppose we have a polygon study_area (a Shapely Polygon). We can select only the points
that lie within it.
python
points_inside = gdf[within_mask]
2. Name your geometry column 'geometry' unless there’s a compelling reason not to.
It avoids having to specify geometry= every time.
4. Use pandas filtering for attribute queries; use predicate methods for spatial
queries.
5. Vectorized operations are faster than looping; use them whenever possible.
6. When plotting (we’ll cover later), the geometry column is automatically used.
Let’s write a script geopandas_intro.py that uses a sample dataset (Natural Earth) or a
simple CSV with coordinates. For the exercise, we will create a GeoDataFrame from a
manually defined list, simulating a real-world workflow.
python
data = {
print(gdf_cities.head())
print("Big cities:")
print(big_cities[["city", "pop_millions"]])
east_china_poly = Polygon([
])
mask = gdf_cities.within(east_china_poly)
print(gdf_cities[mask][["city"]])
# Write to GeoJSON
big_cities.to_file("big_cities.geojson", driver="GeoJSON")
Problem Statement:
You have a dataset of weather stations (attributes: id, elevation_m), and their coordinates
are given as lists. You also have a polygon representing a mountainous area of interest. You
need to:
text
ids = ["S1","S2","S3","S4","S5"]
2. Define a polygon mountain_zone with vertices: (115.9, 39.7), (116.7, 39.7), (116.7,
40.1), (115.9, 40.1), (115.9, 39.7). (Rough rectangle).
4. Filter stations with elevation > 1000 meters. Print the count and the filtered
GeoDataFrame (just id and elevation_m columns).
5. Spatially filter stations that are inside the mountain_zone. Print the resulting station
IDs.
6. Create a buffer of 0.1 degrees around each station in the full GeoDataFrame (you
can assign back or create a new column buffer_geom). Then compute the total area
(in square degrees) of all buffers combined? Actually, compute the area of each
buffer and print the average buffer area.
Requirements:
Why this exercise: You’ll practice GeoDataFrame creation, CRS assignment, attribute
filtering, spatial filtering, geometry operations (buffer, area), and file output. These form the
core of any vector data processing script.
Self-check:
• Stations inside mountain_zone: Check which ones fall within the polygon. Lons
115.9–116.7, lats 39.7–40.1. S1 (116.1,39.8) inside; S2 (116.4,40.0) inside (border?
(116.4,40.0) is exactly on northern boundary? It's 40.0, which is on the line,
but within might return False for boundary. We'll see; S2 might be considered not
inside. S3 (115.8,39.6) outside; S4 (116.3,39.9) inside; S5 (116.6,40.2) outside (40.2
> 40.1). So expected inside: S1 and S4 (S2 maybe not). Accept whatever Shapely
gives.
• Buffer area of 0.1 deg: area = π * 0.1² ≈ 0.0314 square degrees. Average buffer area
should be that per station.
Proceed to write the script. In the next lesson, we’ll learn about coordinate reference
systems in depth and how to project data for accurate measurements.I'm pleased with the
content developed so far. The lesson covers GeoDataFrame creation, reading/writing, and
basic filtering. The exercise reinforces these concepts. I should now provide the final
response to the [Link] now elevate our geospatial coding to the next level: working
with vector data that has both geometry and attributes. You already know how to create
Shapely geometries. You know how to handle lists, dictionaries, and files. Now we bring
them together into a single, powerful data structure called a GeoDataFrame, provided by
the library GeoPandas. This is the cornerstone of vector analysis in the Python geospatial
ecosystem.
Recall our earlier work: we stored station attributes in dictionaries and coordinates in
separate tuples. If we wanted to filter stations by temperature and then check which are
inside a polygon, we had to juggle two separate structures—the list of attributes and the list
of Shapely objects. This quickly becomes messy and error-prone.
• pandas DataFrame – A 2-dimensional, tabular data structure with rows and labeled
columns, capable of holding mixed data types (numbers, strings, booleans). It’s like
an Excel spreadsheet inside Python.
A GeoDataFrame is a pandas DataFrame with a geometry column. It inherits all the power
of pandas for data manipulation (filtering, grouping, aggregating, merging) and adds spatial
operations (spatial joins, intersection, buffering, distance calculations) that operate on the
geometry column.
Key concepts:
• Geometry Column: The active geometry column (you can have multiple, but one is
active). You can access it with [Link]. It’s a GeoSeries, which is a
pandas Series of Shapely geometries with CRS information.
Without GeoPandas, every spatial analysis would require custom loops over lists of
Shapely objects and manual attribute management. GeoPandas provides a unified,
high-performance framework that handles millions of features efficiently (using underlying
GEOS and pygeos/rtree for spatial indexing). It also seamlessly reads and writes dozens of
vector file formats (Shapefile, GeoJSON, GeoPackage, PostGIS, etc.). In professional
geospatial Python, GeoPandas is as fundamental as a hammer to a carpenter.
I’ve seen researchers write 500-line scripts to load, filter, and intersect two shapefiles,
when GeoPandas can do it in 10 lines. Not only is the GeoPandas code shorter, it’s also
more readable and less likely to contain bugs. By mastering GeoPandas, you can stop
worrying about the mechanics of data handling and focus on the spatial questions that
interest you. You also become part of a massive community where solutions are widely
shared.
GeoPandas depends on Shapely, Fiona, PyProj, and other packages. The recommended
installation via conda (as we set up in Lesson 0.1) ensures all dependencies are
compatible. If you haven’t installed GeoPandas yet:
bash
python
python
import pandas as pd
python
A GeoDataFrame can be created in several ways. The most instructive is to build one from a
dictionary or a list of dictionaries, where one column contains Shapely geometry objects.
python
# Build dictionary
data = {
"code": station_codes,
"temperature": temperatures,
print(gdf)
Output:
text
Notice:
• The CRS is set to EPSG:4326 (WGS84 latitude/longitude). Always set the CRS when
creating a GeoDataFrame; without it, spatial operations that depend on CRS
(like to_crs) will complain.
python
features = [
The geometry column is automatically detected if its name is exactly 'geometry'. If your
geometry column has a different name (e.g., 'location'), you can specify it with
the geometry parameter:
python
GeoPandas can read most vector formats using the read_file() function, which is powered
by Fiona.
python
gdf = gpd.read_file("path/to/[Link]")
python
gdf = gpd.read_file("[Link]
python
gdf.to_file("[Link]", driver="GeoJSON")
For Shapefile, the driver is automatically selected, but a directory will be created. Better to
use GeoPackage for modern projects.
python
# Combine conditions
Suppose we have a polygon study_area (a Shapely Polygon). We can select only the points
that lie within it.
python
points_inside = gdf[within_mask]
2. Name your geometry column 'geometry' unless there’s a compelling reason not to.
It avoids having to specify geometry= every time.
3. Use gpd.read_file() and gdf.to_file() for I/O – simple and robust.
4. Use pandas filtering for attribute queries; use predicate methods for spatial
queries.
5. Vectorized operations are faster than looping; use them whenever possible.
6. When plotting (we’ll cover later), the geometry column is automatically used.
Let’s write a script geopandas_intro.py that uses a sample dataset (Natural Earth) or a
simple CSV with coordinates. For the exercise, we will create a GeoDataFrame from a
manually defined list, simulating a real-world workflow.
python
data = {
print(gdf_cities.head())
print("Big cities:")
print(big_cities[["city", "pop_millions"]])
east_china_poly = Polygon([
])
mask = gdf_cities.within(east_china_poly)
print(gdf_cities[mask][["city"]])
# Write to GeoJSON
big_cities.to_file("big_cities.geojson", driver="GeoJSON")
Problem Statement:
You have a dataset of weather stations (attributes: id, elevation_m), and their coordinates
are given as lists. You also have a polygon representing a mountainous area of interest. You
need to:
text
ids = ["S1","S2","S3","S4","S5"]
2. Define a polygon mountain_zone with vertices: (115.9, 39.7), (116.7, 39.7), (116.7,
40.1), (115.9, 40.1), (115.9, 39.7). (Rough rectangle).
4. Filter stations with elevation > 1000 meters. Print the count and the filtered
GeoDataFrame (just id and elevation_m columns).
5. Spatially filter stations that are inside the mountain_zone. Print the resulting station
IDs.
6. Create a buffer of 0.1 degrees around each station in the full GeoDataFrame (you
can assign back or create a new column buffer_geom). Then compute the total area
(in square degrees) of all buffers combined? Actually, compute the area of each
buffer and print the average buffer area.
Requirements:
Self-check:
• Stations inside mountain_zone: Check which ones fall within the polygon. Lons
115.9–116.7, lats 39.7–40.1. S1 (116.1,39.8) inside; S2 (116.4,40.0) inside (border?
(116.4,40.0) is exactly on northern boundary? It's 40.0, which is on the line,
but within might return False for boundary. We'll see; S2 might be considered not
inside. S3 (115.8,39.6) outside; S4 (116.3,39.9) inside; S5 (116.6,40.2) outside (40.2
> 40.1). So expected inside: S1 and S4 (S2 maybe not). Accept whatever Shapely
gives.
• Buffer area of 0.1 deg: area = π * 0.1² ≈ 0.0314 square degrees. Average buffer area
should be that per station.
Proceed to write the script. In the next lesson, we’ll learn about coordinate reference
systems in depth and how to project data for accurate measurements.
Here is the complete script, annotated line-by-line so you understand the reasoning behind
each stroke.
python
#!/usr/bin/env python3
"""
"""
import geopandas as gpd
def main():
# Build a list of Point geometries, ensuring lon first (x) then lat (y)
mountain_zone = Polygon([
(115.9, 39.7),
(116.7, 39.7),
(116.7, 40.1),
(115.9, 40.1),
])
# 3. Print the full GeoDataFrame
print("Full GeoDataFrame:")
print(stations_gdf)
print()
print(high_stations[["id", "elevation_m"]])
print()
inside_mask = stations_gdf.within(mountain_zone)
inside_stations = stations_gdf[inside_mask]
print(inside_stations[["id"]])
print()
buffers = stations_gdf.buffer(0.1)
avg_area = buffer_areas.mean()
print()
# 7. Save high-elevation stations to GeoJSON
high_stations.to_file("high_stations.geojson", driver="GeoJSON")
print("Exported high_stations.geojson")
if __name__ == "__main__":
main()
• CRS: Always "EPSG:4326" for lat/lon data. Without it, the GeoDataFrame would be
“crs-less” and many operations would warn or fail.
If your script behaves similarly, excellent. If you had a different structure but achieved the
same outputs, that’s fine as long as you used vectorized methods.
In a regular table join (like in SQL or pandas merge), you combine two tables based on a
common attribute column: for example, join a table of city names to a table of city
populations using the city ID. That’s an attribute join.
A spatial join does the same thing, but the “key” is location. It takes two GeoDataFrames
and combines rows that have a specified spatial relationship. For each geometry in the left
GeoDataFrame, it looks for all geometries in the right GeoDataFrame that satisfy the
predicate (e.g., intersects, within, contains) and attaches the attributes of the matching
right row(s) to the left row. If multiple right rows match, the left row is duplicated
(one-to-many join). If no right row matches, you can choose to keep or drop the left row.
GeoPandas implements spatial joins via the [Link]() function. Under the hood, it builds
a spatial index (R-tree) on the right GeoDataFrame for efficiency, so it scales well to large
datasets.
Earlier, to assign a district name to each station, you would have to loop over every station,
then loop over every district polygon, test within, and if true, copy the name. For 10,000
stations and 500 districts, that’s 5 million predicate tests, and you’d also have to manually
align attributes. With sjoin, you write:
python
It’s a single line, it’s fast (using spatial index), and it’s correct. It handles multiple matches
gracefully. The sjoin function is your primary tool for enriching data with location-based
attributes.
python
The result is a new GeoDataFrame containing the columns of both left and right
DataFrames (if how="inner" or "left"), with the right’s geometry column by default dropped
to avoid duplication. The left’s geometry is retained.
python
# Points (cities)
cities_gdf = [Link]({
}, crs="EPSG:4326")
# Polygons (districts)
districts_gdf = [Link]({
"geometry": [Polygon([(0,0),(2,0),(2,2),(0,2)]),
Polygon([(3,3),(5,3),(5,5),(3,5)])]
}, crs="EPSG:4326")
# Spatial join: which district is each city in?
print(joined)
Output:
text
0 A POINT (1 1) 0 North
1 B POINT (3 3) 1 South
City C is outside both polygons, so it’s dropped (inner join). If we used how="left", city C
would remain with district=NaN.
Notice the index_right column – it indicates the index of the matching row in the right
GeoDataFrame. You can suppress it with lsuffix and rsuffix parameters if column names
clash, but the default behavior is safe.
If a left geometry matches multiple right geometries (e.g., a point lies on the boundary of
two polygons, depending on the predicate), the left row will be duplicated for each match.
This is by design. Use how and the predicate carefully to avoid unwanted duplicates.
sjoin automatically builds a spatial index on the right GeoDataFrame. If you are performing
multiple joins with the same right DataFrame (e.g., inside a loop), precompute the index
yourself using right_df.sindex and then use .query() for each geometry, but sjoin is already
optimized for one-shot joins.
1. Ensure both GeoDataFrames have a CRS – if CRSs differ, sjoin will raise an error.
Reproject one to match the other using to_crs().
3. Beware of boundary cases – points on polygon edges may not satisfy "within";
use "intersects" if you consider touching as inside.
4. Handle multiple matches – if you need a one-to-one join, you can later aggregate
the duplicate rows using dissolve or groupby (next lesson).
6. Always inspect the result – [Link] and [Link]() to verify expected row
counts.
python
# States (simplified)
states = [Link]({
"geometry": [
}, crs="EPSG:4326")
# Earthquake points
quakes = [Link]({
}, crs="EPSG:4326")
# Spatial join
print(quake_state[["mag", "state"]])
Output:
text
mag state
0 4.5 California
1 5.2 Nevada
2 3.8 Oregon
Problem Statement:
You are given two GeoDataFrames representing field survey plots and species sightings.
• P1: (2,3), P2: (6,2), P3: (2,8), P4: (8,4), P5: (3,5)
Note: point P5 lies exactly on the boundary between Plot 1 and Plot 3 at (3,5) (the boundary
line at y=5).
1. Create the two GeoDataFrames (with appropriate CRS, e.g., "EPSG:32610" for a
local UTM, but for simplicity just use a Cartesian CRS like "EPSG:3857" or even an
undefined CRS; we’ll set crs="EPSG:3857" as placeholder).
3. Notice what happens to P5. It is on the boundary, so "within" will exclude it. Now
perform the join using predicate "intersects" with how="left". Print the result. Does
P5 appear? How many rows? Discuss in a comment.
4. Calculate the number of sightings per plot. After joining with intersects, use
the groupby method (a core pandas operation) to count
sightings: [Link]("plot_id").size(). Print these counts.
Requirements:
Why this exercise: You’ll practice the spatial join with different predicates, observe
boundary effects, perform a simple aggregation, and export the result. Counting points per
polygon is one of the most common geospatial queries.
Self-check:
• Inner within join: P1 (inside Plot1), P2 (inside Plot2), P3 (inside Plot3), P4? (8,4) is
inside Plot2? Plot2 spans x=5–10, y=0–5. (8,4) is inside. So inner join gives 4 rows. P5
excluded.
• Left intersects join: P5 now appears, with polygon? It intersects both Plot1 and Plot3
because it's on their shared boundary. Since sjoin duplicates, P5 will appear twice!
That's important to note. So number of rows = previous 4 + 2 for P5 = 6.
Write the script now. Next, we’ll explore dissolve and aggregate operations to merge
geometries within groups.
Solution to Exercise: spatial_join_lab.py
Here is the complete script with extensive comments to illuminate every decision.
python
#!/usr/bin/env python3
"""
"""
def main():
plots_data = [
sightings_data = [
{"sighting_id": "P5", "geometry": Point(3, 5)}, # on boundary between Plot1 and Plot3
print(inner_within[["sighting_id", "plot_id"]])
print(left_intersects[["sighting_id", "plot_id"]])
# P5 now appears, but twice! It intersects both Plot1 and Plot3 at the shared edge.
counts = left_intersects.groupby("plot_id").size()
print(counts)
print()
# 6. Write joined result to GeoPackage
print("Exported sightings_per_plot.gpkg")
if __name__ == "__main__":
main()
• The group counts using intersects are: Plot1=2, Plot2=2, Plot3=2. This shows how
boundary cases can inflate counts if not handled intentionally.
This exercise reinforces a critical rule: Always inspect your join results, especially when
boundaries matter. If you need a one-to-one relationship, you might need to adjust
predicates or post-process duplicates.
Imagine you have a GeoDataFrame of land-use parcels, each with a land_use category
(forest, urban, water). You want to create a single merged polygon for each category—the
union of all parcels of that type. This operation is called dissolve. It performs two actions
simultaneously:
2. Aggregates the geometry column by merging all geometries within each group into
a single geometry (usually using unary_union, which is a fast union of all shapes).
Optionally, you can also aggregate other columns (e.g., sum area, average value).
The result is a new GeoDataFrame with one row per unique group value, with the geometry
being the union of all features in that group. It’s the spatial equivalent of a SQL GROUP
BY with a spatial aggregate function.
Before dissolve, you would have to loop over categories, subset the GeoDataFrame, and
call unary_union yourself—then build a new GeoDataFrame from
scratch. dissolve encapsulates this pattern into a single, readable, and optimized line. It
also handles the aggregation of other columns elegantly.
python
• by – The column(s) to group by. Can be a single column name (string) or a list of
column names. If by=None (default), it dissolves all geometries into one single
feature.
• aggfunc – How to aggregate other (non-geometry) columns. Default is 'first' (take the
first value in the group). You can pass a dictionary mapping column names to
aggregation functions, e.g., {'population':'sum', 'area':'mean'}.
• as_index – If True (default), the group column becomes the index of the resulting
GeoDataFrame. Set as_index=False to keep it as a column.
The geometry column is aggregated using unary_union (the union of all geometries in the
group). This operation dissolves internal boundaries, merging adjacent or overlapping
shapes.
Example: Dissolve parcels by land use
python
# Sample parcels
gdf = [Link]({
'geometry': [
Polygon([(0,0),(1,0),(1,1),(0,1)]),
Polygon([(1,0),(2,0),(2,1),(1,1)]),
Polygon([(0,1),(1,1),(1,2),(0,2)]),
Polygon([(1,1),(2,1),(2,2),(1,2)]),
Polygon([(2,0),(3,0),(3,1),(2,1)])
}, crs="EPSG:3857")
print(dissolved)
Output:
text
geometry area_ha
land_use
Notice:
• The geometry column contains the merged polygon(s). Adjacent forest parcels
merged into a single larger polygon.
If you don’t need attribute aggregation, just use [Link](by='land_use') and it will
discard other columns (or keep the first value).
python
all_merged = [Link]()
This creates a single-row GeoDataFrame with the union of all geometries. Useful for
creating a study area boundary from individual tiles.
You can pass a dictionary to aggfunc to specify different aggregation functions per column.
python
Available functions: 'sum', 'mean', 'min', 'max', 'first', 'last', 'count', or any custom function.
python
7. Performance Considerations
dissolve uses unary_union under the hood, which is efficient but can be slow for very
complex geometries with millions of vertices. For large datasets, consider simplifying
geometries first ([Link](tolerance)) if appropriate. Otherwise, it’s fine for most tasks.
8. “How to Write” Rules for Dissolve
1. Always check the CRS before dissolving; if you need area/length calculations later,
make sure the CRS is projected (meters), not geographic (degrees).
4. Reset index if needed – After dissolve, the group column becomes the index.
Use reset_index() or as_index=False to convert back to a column.
5. Watch for None geometries – If any geometry in a group is invalid or None, the
union may fail or produce unexpected results. Clean your data first.
python
# muni = gpd.read_file("[Link]")
muni = [Link]({
'geometry': [
Polygon([(0,0),(1,0),(1,1)]),
Polygon([(1,0),(2,0),(2,1),(1,1)]),
Polygon([(0,1),(1,1),(1,2)]),
Polygon([(1,1),(2,1),(2,2)]),
Polygon([(2,0),(3,0),(3,1)])
Output shows one row per province, with merged geometry, total population, and area in
km².
Problem Statement:
3. Dissolve by species:
6. Create a second dissolve without by (i.e., merge all stands into one). Print its area in
hectares.
Requirements:
Why this exercise: You’ll practice dissolving by attribute, aggregating numeric columns,
computing area in metric units, and exporting. This is a building block for land management
and ecological analysis.
Self-check:
• Original: 5 stands.
• Dissolved by species: should yield 3 rows: oak (stands A, C, E), pine (B), birch (D).
o Oak total volume: 1200+1500+900 = 3600. Area: union of three squares each
2x2 = 4 m²? Wait, coordinates 0-2,2-4 etc. Each stand is 2x2 = 4 sq meters,
but they are not all adjacent: A (0,0-2,2), C (0,2-2,4) touching A, E (4,0-6,2)
separate. So union area = area of A+C + area of E = (4+4? Actually A and C
merge into a 2x4 rectangle = 8, plus E 4 = 12 m². In hectares: 12 / 10000 =
0.0012 ha.
So area numbers will be small because units are small squares. That's fine; the point is the
process.
Proceed to write the script. Next, we’ll learn overlay operations (intersection, union,
identity) to combine multiple layers.
python
#!/usr/bin/env python3
"""
Computes total volume and area per species, and exports the result.
"""
def main():
stands_data = [
print()
species_dissolved = [Link](
by="species",
aggfunc={"volume_m3": "sum"},
as_index=False
print()
print()
species_dissolved.to_file("species_units.geojson", driver="GeoJSON")
print("Exported species_units.geojson")
if __name__ == "__main__":
main()
• Full Union: Calling dissolve() without by merges every geometry into a single row.
Note that .area on the resulting GeoDataFrame returns a GeoSeries, so we
index [0] to get the scalar value.
• Clean Output: Only relevant columns are printed, and the final file export is in an
open format (GeoJSON).
This pattern is reusable for any aggregation task where you need to combine features by a
categorical attribute.
In geoprocessing, overlay operations take two polygon (or mixed) layers and produce a
new layer where the boundaries are split and recombined based on the spatial
relationship. Think of it as a set of Boolean operations on space: given two layers A and B,
you can compute:
• Symmetric Difference ((A − B) ∪ (B − A)): Areas covered by only one of the two.
• Identity (A identity B): All features of A are retained, but split by overlapping features
of B—essentially an intersection that keeps the full extent of A.
GeoPandas provides these operations via the overlay() function, which is conceptually
similar to sjoin but outputs a new GeoDataFrame with geometries that are the result of the
spatial operation, and with attributes combined from both parent layers. This is how you do
“cookie-cutting,” “site selection,” or “impact area” analysis.
Unlike a simple dissolve that merges within the same layer, overlay combines two different
layers. You might have a layer of flood zones and a layer of property parcels; an
intersection will give you the portions of each parcel that are within the flood zone, with
attributes of both. A union gives you all areas, split into unique combinations of attributes.
Identity can be used to attach the flood zone information to parcels while keeping the full
parcel shapes (splitting them if they cross zone boundaries).
These operations are fundamental to vector analysis and are analogous to map algebra for
rasters, but for discrete vector features.
Overlay operations are the engine behind site suitability modeling, risk assessment (e.g.,
population in flood zones), and administrative unit analysis.
3. Syntax: [Link]()
python
python
• 'intersection' – keeps only the overlapping areas, with attributes from both.
• 'union' – keeps all areas, creating new polygons for all combinations.
Important: Both GeoDataFrames must have the same CRS, and usually you want polygon
layers. Points and lines can be used with some operations (like intersection of points with
polygons), but the output geometry type may change (e.g., intersecting a line with a
polygon yields a line or point). For simplicity, we focus on polygon-polygon overlays.
Example: Intersection of two polygons
python
counties = [Link]({
"geometry": [
}, crs="EPSG:3857")
flood_zone = [Link]({
"zone": "Flood",
}, crs="EPSG:3857")
print(intersect)
Output:
text
The resulting GeoDataFrame has columns from both parent layers, and the geometries are
the actual overlapping polygons. Here the flood zone crosses both counties, so we get two
rows.
Union Example
python
print(union)
Union produces all distinct regions created by the overlay. There will be areas that belong
only to County A, only to County B, only to the flood zone (outside counties, if any), and the
overlapping parts. Each row will have the combination of attributes; where attributes are
missing (e.g., a part of County A not in flood zone), the zone column will be None.
Difference Example
python
print(diff)
This returns the parts of the counties that are not in the flood zone, retaining the county
attributes.
Identity Example
python
print(identity)
Identity is like intersection but keeps the full extent of the first layer. It splits the first layer
by the second, so you get pieces of counties that are either inside or outside the flood zone.
The zone column will be populated for the intersected parts, and None for the
non-overlapping parts. This is extremely useful for enriching one layer with attributes from
another without losing area.
• Invalid geometries (self-intersections, etc.) can cause overlay to fail. Fix them
with .buffer(0) or .make_valid() (Shapely 2.0+).
1. Check CRS consistency – both layers must have the same CRS. Reproject if
needed.
2. Use how='identity' for attribute transfer – it retains the entire first layer, which is
often what you want.
3. Aggregate results – after overlay, you often need to dissolve by category and sum
areas.
We’ll demonstrate a classic municipal analysis: find the area of each land-use type within
each zoning district.
python
parcels = [Link]({
"land_use": ["Residential", "Commercial", "Park"],
"geometry": [
}, crs="EPSG:32633") # metric
# Zoning districts
zoning = [Link]({
"geometry": [
}, crs="EPSG:32633")
parcel_zone["area_sqm"] = parcel_zone.[Link]
Problem Statement:
You have two polygon layers representing soil types and protected areas. You need to
calculate the area of each soil type that falls within protected areas, and also create a
“combined” map of both layers.
3. Perform an intersection overlay to get soil polygons within protected areas. Print
the result with columns soil, protected, and area in square meters (as a new
column area_sqm).
4. Using the intersection result, dissolve by protected to get the total area of each
protected area (regardless of soil). Print these areas.
5. Perform a union overlay to create a combined layer of all soil and protected areas.
The result will have many small pieces with attributes from both. Compute the total
area of the union result (sum of all pieces) and print it.
6. Perform an identity overlay with soil as first layer and protected areas as second.
This will keep all soil areas, split by protected areas. Print the resulting area of the
soil portion that is outside any protected area (hint: filter where protected is None,
then sum areas).
Requirements:
Self-check:
• Intersection will produce pieces of soil squares overlapping the reserve and park.
• Park (5,1)-(9,5) overlaps Clay (4,0)-(8,4): intersection (5,1)-(8,4) area 3*3=9; overlaps
Loam? Loam (0,4)-(4,8) doesn't intersect Park. So Park area 9 m².
• Intersection result has rows: Silt-Reserve (4), Clay-Reserve (4), Loam-Reserve (4),
Clay-Park (9) → 4 rows.
• Union: total area = sum of all distinct regions. Initially soil total = 3*16=48 m² (each
4x4). Protected areas: Reserve 4x4=16, Park 4x4=16, but they overlap each other?
Park (5,1)-(9,5) and Reserve (2,2)-(6,6) overlap in (5,2)-(6,5) area 1*3=3. So union
area = 48+16+16 - (overlap of Reserve with soils? Actually overlay union will
correctly compute non-overlapping areas, sum will be ≤ total extent.
• Identity: splits soils by protected, outside parts have protected=None. Outside area
= total soil area - area of intersection = 48 - (12+9) = 27 m² (since intersection
covered 12+9=21, but check: reserve and park overlap in the clay region? Actually
park and reserve overlap each other, but intersection overlay we computed each
protected area separately, so the sum of intersection areas is 21. However, the two
protected areas overlap in the clay region: Clay-Reserve overlap 4 m², Clay-Park
overlap 9 m², but part of that clay area is counted in both? Wait, the park and
reserve are separate polygons, so their intersection with soils are independent. The
soil polygons are split, so there is no double counting; the identity will split the clay
soil into parts: one part only in reserve, one only in park, one in both (if they overlap).
Since reserve and park overlap each other, that overlapping region would be
assigned to both protected areas? Identity will split by the combined boundaries. So
the outside area computed as sum of areas where protected is None will be exactly
the soil area not covered by either protected area. This is 48 - area of union of
protected areas intersected with soils. It will be less than 27 if there is overlap. We'll
accept whatever geopandas calculates. The point is to practice.
python
#!/usr/bin/env python3
"""
"""
def main():
soil = [Link]({
"geometry": [
}, crs="EPSG:32610")
# 2. Create protected areas (Reserve and Park)
protected = [Link]({
"geometry": [
}, crs="EPSG:32610")
intersection["area_sqm"] = [Link]
print()
print(protected_total[["geometry", "area_sqm"]])
print()
union["area_sqm"] = [Link]
total_union_area = union["area_sqm"].sum()
print(f"Total union area: {total_union_area:.2f} m²")
print()
identity["area_sqm"] = [Link]
outside_protected = identity[identity["protected"].isna()]
outside_area = outside_protected["area_sqm"].sum()
print()
print("Exported soil_identity.gpkg")
if __name__ == "__main__":
main()
• Union creates a complex patchwork of all soil and protected boundaries. Summing
areas gives the total extent of the combined layers.
• Identity keeps the full soil layer but splits it. Filtering where protected is NaN yields
the unprotected soil areas.
• All calculations are in square meters because the CRS is projected (EPSG:32610),
which leads us directly into our next topic.
LESSON 4.1 — Coordinate Reference Systems and Projections: The Shape of the Earth
in Code
1. Concept & Theory: Why a Round Earth Looks Flat on Your Screen
The Earth is not a perfect sphere; it is an oblate spheroid, slightly flattened at the poles and
bulging at the equator. To represent its surface on a flat map or in a planar geometry engine
like Shapely, we must define a mathematical transformation from the curved,
three-dimensional Earth to a two-dimensional plane. This transformation is called a map
projection.
A Coordinate Reference System (CRS) is a complete set of rules that defines how
coordinates relate to actual locations on Earth. It consists of:
• A datum – the reference ellipsoid (model of the Earth’s shape) and its position
relative to the Earth’s center. WGS84 is a global datum; NAD83 is used in North
America.
In Python, the library PyProj (wrapped by GeoPandas’ to_crs()) handles all the complex
mathematics, translating between CRSs using the PROJ library (the same engine behind
QGIS, PostGIS, and many other GIS).
• To calculate meaningful distances, areas, or angles, you must use a projected CRS
appropriate for your region. Computing length in degrees leads to gross errors (1
degree of longitude varies from ~111 km at the equator to 0 km at the poles).
• Different data providers use different CRSs; to overlay them, they must be
transformed to a common CRS.
• Some operations (like buffers) in Shapely are purely planar; they will not work as
expected if you feed them lat/lon degrees expecting meters.
The Role of PyProj and GeoPandas
• PyProj is the low-level library that can transform individual coordinate pairs or
whole arrays. You can use it directly with [Link], but GeoPandas
simplifies this.
• GeoPandas stores the CRS in the .crs attribute and provides the .to_crs() method to
reproject entire GeoDataFrames.
2. “Why This?” – The Difference Between 100 Degrees and 100 Meters
Imagine you are analyzing deforestation within a 10-km radius of a ranger station. If your
data is in EPSG:4326 (WGS84 lat/lon), a buffer of 0.1 degrees will not be 10 km everywhere;
near the poles it will be a tiny fraction of that. Your analysis would be completely wrong.
Conversely, if you project to a local UTM zone (meters), a buffer of 10000 meters is exactly
10 km everywhere. The choice of CRS is not optional technicality; it is the foundation of
accurate spatial reasoning.
In GeoPandas, [Link] returns a [Link] object. You can print it to see its definition:
python
gdf = [Link](
crs="EPSG:4326"
print([Link])
# Output: EPSG:4326
• 326xx – UTM zones in the northern hemisphere (e.g., 32610 for zone 10N). Southern
hemisphere: 327xx.
You have already done this: [Link](data, crs="EPSG:4326"). You can also
assign to .crs after creation:
python
[Link] = "EPSG:4326"
But if the coordinates are already in that CRS, you are just telling GeoPandas the truth. If
you assign a wrong CRS, reprojection will be wrong.
The most common operation: convert from geographic to projected (or between
projections).
python
gdf_web = gdf.to_crs("EPSG:3857")
print(gdf_web.[Link]())
You can pass an EPSG code, a PROJ string, or a [Link] object. GeoPandas calls PyProj
under the hood to transform every coordinate in the geometry column.
Performance tip: Reprojection of large GeoDataFrames can be slow. If you only need to
compute distances or areas on the fly, consider converting just the coordinates and doing
the operation without persisting the full reprojected dataset (though for simplicity we’ll
project the whole GeoDataFrame).
• UTM zones are a safe default: each zone covers 6° of longitude. You can compute
the zone for a given longitude: zone = int((lon + 180) / 6) + 1. For northern
hemisphere, EPSG = 32600 + zone; southern = 32700 + zone.
For now, UTM is enough. We’ll use a helper function to select the right UTM zone.
Let’s write a script that demonstrates the danger of degree-based buffers and the fix using
projection.
python
buffer_deg = [Link](0.1) # 0.1 degree, expecting ~7.45 km radius? Actually it's large.
point_proj = point.to_crs("EPSG:32648")
buffer_deg_geo = buffer_deg.to_crs("EPSG:4326")
buffer_m_geo = buffer_m.to_crs("EPSG:4326")
You’ll see that the degree-based buffer is not a circle on the ground (its shape and area
depend on latitude), whereas the metric buffer, when reprojected back to geographic, is a
true circle on the ellipsoid (more accurately, a geodesic circle). The areas will differ.
1. Always know the CRS of your input data. If read from a file, [Link] tells you.
2. Assign CRS if missing but only if you are sure what it is. Never guess.
4. Use UTM zones for local or regional analysis. For global analysis, consider
equal-area projections (e.g., "+proj=cea").
5. Keep original geographic copy if you need to write to GeoJSON (which uses
WGS84) or visualize on web maps (Web Mercator).
6. Project once and reuse – don’t reproject the same data multiple times in a pipeline.
Problem Statement:
You are given a list of three weather stations in different parts of the world, with their
geographic coordinates (WGS84). You need to compute accurate distances between them
and buffer zones, using UTM projections.
2. Write a function get_utm_epsg(lon, lat) that returns the EPSG code for the UTM zone
containing the given point. Use the formula: zone = int((lon + 180) / 6) + 1, and if lat
>= 0, EPSG = 32600 + zone; else 32700 + zone. (Assume latitudes are between -80
and 84, where UTM is valid).
3. For each station, reproject it to its own UTM zone and create a buffer of 100 km
(100,000 meters). Store the buffered polygons in a list.
4. Convert each buffer back to EPSG:4326 and combine them into a single
GeoDataFrame buffers_geo. Print the area of each buffer in square kilometers
(compute area in the original projected CRS before reprojecting, or compute after?
Better to compute area in the projected CRS and then record it; you can store the
area as a column before reprojecting). Print each station name and its buffer area in
km².
5. Compute the accurate straight-line distance (in km) between Station A and Station
B. You can either:
o Reproject both to a common UTM zone (choose the zone of one of them, or a
zone in between, e.g., halfway meridian) and compute Euclidean distance.
o Or use the Haversine formula (we coded that earlier). Since we’re focusing on
projections, use the projection method: find a suitable UTM zone for the
midpoint of the two longitudes, reproject both points, and compute the
distance. Compare with the Haversine distance (import
your geo_functions or re-implement). Print both distances.
6. Write a comment discussing why the UTM distance might differ slightly from
Haversine and under what circumstances it’s acceptable.
Requirements:
Self-check:
• Buffer areas should be close to π * 100² ≈ 31,416 km² (since it’s a 100 km radius
circle).
• Distance between NYC and SF: Haversine gives ~4,130 km. UTM method using a
central zone should be close, but may vary due to projection distortion over long
distances. That’s okay; note the difference.
• The function get_utm_epsg for NYC (lon -73.9) = zone 18N (EPSG:32618); SF (lon -
122.4) = zone 10N (32610); Tokyo (lon 139.7) = zone 54N (32654).
Proceed to write the script. Next lesson, we’ll dive into raster data and how Python handles
imagery and elevation [Link] will now review the solution to the projection exercise,
solidifying your ability to work with coordinate reference systems dynamically. Then, we
shall cross the threshold from vector to raster—the world of pixels, grids, and continuous
fields. Raster data powers satellite imagery, digital elevation models, climate surfaces, and
land-cover classification. Mastering it will complete your geospatial Python foundation.
python
#!/usr/bin/env python3
"""
"""
import math
"""
Return the EPSG code of the UTM zone for the given longitude and latitude.
Args:
Returns:
"""
if lat >= 0:
else:
R = 6371.0
a = [Link](dlat/2)**2 + [Link](lat1)*[Link](lat2)*[Link](dlon/2)**2
return R * 2 * [Link]([Link](a))
def main():
data = {
buffer_areas = {}
buffers_geo_list = []
name = row["name"]
point_proj = point_geo.to_crs(f"EPSG:{epsg}")
buffer_proj = point_proj.buffer(100_000)
area_m2 = buffer_proj.area[0]
area_km2 = area_m2 / 1e6
buffer_areas[name] = area_km2
buffer_geo = buffer_proj.to_crs("EPSG:4326")
buffer_geo["name"] = name
buffers_geo_list.append(buffer_geo)
buffers_all = [Link](
[Link](buffers_geo_list, ignore_index=True),
crs="EPSG:4326"
print()
pt_nyc_proj = pt_nyc.to_crs(f"EPSG:{common_epsg}")
pt_sf_proj = pt_sf.to_crs(f"EPSG:{common_epsg}")
dx = pt_nyc_proj.x[0] - pt_sf_proj.x[0]
dy = pt_nyc_proj.y[0] - pt_sf_proj.y[0]
# Haversine distance
print("(Small differences are expected due to projection distortion over long distances.)")
buffers_all.to_file("station_buffers.geojson", driver="GeoJSON")
print("\nExported station_buffers.geojson")
if __name__ == "__main__":
import pandas as pd # for [Link]
main()
Key Points
• Buffering workflow: extract row → project → buffer → compute area in projected CRS
→ reproject back to 4326. Always compute area before reprojecting, because area in
4326 is in square degrees.
• Common UTM for distance: choosing the UTM of the midpoint minimizes distortion
for both points. The Haversine distance provides the “true” great-circle distance on
the sphere; the UTM distance is a planar approximation that works well within a
zone (6° width). Over the large longitudinal span of NYC–SF (~48°), UTM distorts
significantly; yet we still used the method to demonstrate the principle. In practice,
you’d use a equidistant projection or Haversine for such long distances.
LESSON 5.1 — Raster Data as Arrays: Rasterio, Numpy, and the Pixel World
Vector data (points, lines, polygons) represents discrete objects. Raster data represents
continuous fields as a regular grid of cells (pixels), where each cell stores a numeric value.
Examples:
• Satellite imagery: multi-band (Red, Green, Blue, NIR) where each pixel is a
brightness value.
• Georeferencing metadata (the CRS, and a transformation that maps pixel indices
to geographic coordinates).
• Nodata value (a sentinel indicating missing data, like -9999).
In Python, the primary library for raster I/O is rasterio. It reads raster data efficiently (often
using memory-mapped files) and exposes the pixel arrays as numpy arrays. Numpy is the
fundamental package for numerical computing in Python; it provides the ndarray object—a
fast, memory-efficient multi-dimensional array—along with vectorized mathematical
functions. Rasterio and numpy together allow you to perform pixel-level computations with
the speed of compiled C code.
Imagine you have a 10,000×10,000 pixel satellite image (100 million pixels). With pure
Python loops, multiplying each pixel by 2 would take seconds. Numpy performs the same
operation in a fraction of a second using optimized, vectorized operations (array * 2).
Rasterio understands the spatial context, so you can window-read only the portion you
need, reproject rasters, or even process in parallel.
• Extract pixel values at point locations (e.g., elevation for weather stations).
• Classify land cover from spectral indices (NDVI from Red and NIR bands).
If you haven’t yet: pip install rasterio or conda install -c conda-forge rasterio.
python
import rasterio
[Link] is an Affine object containing six coefficients: (a, b, c, d, e, f). The mapping
from pixel (row, col) to coordinate (x, y) is:
text
x = a * col + b * row + c
y = d * col + e * row + f
Typically, b and d are 0 for north-up rasters; a is pixel width; e is negative pixel height
(north-up).
The .read() method returns a [Link]. You can then apply numpy operations.
python
import numpy as np
# Compute NDVI from Red (band 3) and NIR (band 4) of a Landsat image
red = [Link](3).astype('float32')
nir = [Link](4).astype('float32')
You can create a new GeoTIFF with the same georeferencing by copying the metadata
(profile) and writing a numpy array.
python
[Link]([Link]('float32'), 1)
Important: The array shape must match (bands, rows, cols). For a single band, write a 2D
array; it’s automatically interpreted as band 1.
To read only a portion of a raster (e.g., a city block), use window parameter.
python
This avoids loading the whole image into memory, critical for large datasets.
1. Always use with [Link](...) as src: to ensure the file is properly closed.
4. Use appropriate data types – int16 for integer values, float32 for computations.
Convert with .astype().
5. Profile copying – when creating output rasters, always start from [Link] and
update only what’s needed.
6. Compression – add compress='lzw' (lossless) or 'deflate' for GeoTIFF to reduce file
size.
7. Window read – use Window for large files to keep memory low.
We’ll simulate a small example with fake data to demonstrate the full workflow.
python
import rasterio
import numpy as np
profile = {
'driver': 'GTiff',
'height': 100,
'width': 100,
'count': 3,
'dtype': 'uint16',
'crs': 'EPSG:32633',
[Link](42)
[Link](red, 1)
[Link](green, 2)
[Link](nir, 3)
# Now process it
red = [Link](1).astype('float32')
nir = [Link](3).astype('float32')
out_profile = [Link]
[Link]([Link]('float32'), 1)
Create a new file raster_basics.py. You will generate a synthetic elevation raster, compute
slope, and extract values at points.
Problem Statement:
2. Compute a slope raster (in degrees) from the DEM. Slope at each pixel can be
approximated by the magnitude of the gradient (using central differences) divided by
pixel size.
o Write the slope raster as [Link] (float32, same profile but nodata at borders
where gradient can’t be computed? You can set edge rows/cols to 0 or
nodata). For simplicity, compute slope only for interior pixels (rows 1 to -2,
cols 1 to -2) and keep edges as 0.
4. Extract the elevation and slope values at these points from the rasters you created.
You can use rasterio’s sample method or iterate over points using [Link](point.x,
point.y) to get row/col and then index the array. Print for each point: name, elevation,
slope.
5. Plotting (optional): if you have matplotlib installed, display the slope raster. (We’ll
learn plotting in detail later, but you may attempt a simple imshow).
Requirements:
Why this exercise: You will practice creating synthetic rasters, writing and reading them,
performing a common terrain analysis (slope), and extracting raster values at point
locations—a fundamental building block for zonal statistics and sampling.
Self-check:
• Slope should be highest on steep flanks of the sine/cosine waves (max maybe >
30°).
• Point extraction must return correct elevation and slope values; you can manually
compute the expected elevation from the formula at those coordinates.
Proceed to write the script. Next, we’ll explore raster-vector integration in depth: zonal
statistics, masking, and [Link] will now review the solution to the raster basics
exercise, which ties together numpy array creation, raster writing, slope computation, and
point extraction. Mastering these patterns will give you the confidence to handle any raster
dataset. After that, we will deepen your raster skills by learning raster-vector integration—
how to mask, clip, and compute statistics across zones—a cornerstone of geospatial
analysis.
python
#!/usr/bin/env python3
"""
raster_basics.py – Create a synthetic DEM, compute slope, and extract values at points.
"""
import numpy as np
import rasterio
def main():
pixel_size = 30.0
# Elevation formula
elevation = 1000.0 + 500.0 * [Link](2 * [Link] * col_idx / 100) * [Link](2 * [Link] * row_idx /
100)
elevation = [Link]('float32')
# Profile
profile = {
'driver': 'GTiff',
'height': height,
'width': width,
'count': 1,
'dtype': 'float32',
'crs': 'EPSG:32610',
'transform': transform,
'compress': 'lzw'
[Link](elevation, 1)
print("Created [Link]")
# dz_dx and dz_dy have shape (100, 98) and (98, 100); need to align to interior (98,98)
slope_deg = [Link](slope_rad)
slope_full = np.zeros_like(elevation)
slope_profile = [Link]()
slope_profile.update(dtype='float32', nodata=0)
with [Link]('[Link]', 'w', **slope_profile) as dst:
[Link](slope_full.astype('float32'), 1)
print("Created [Link]")
points_data = {
"geometry": [
Point(500150, 4200150),
Point(501500, 4201500),
Point(502000, 4202000)
pt = [Link]
r, c = [Link](pt.x, pt.y)
elev_val = [Link](1)[r, c]
slope_val = [Link](1)[r, c]
else:
if __name__ == "__main__":
main()
Key Points
• Transform creation: We used an Affine directly, with origin_y set to the top edge.
Since pixel_height is negative, origin_y must be the maximum y (north).
• Raster value extraction: [Link](x, y) returns (row, col). Then we index the
numpy array directly. This is very efficient.
Real-world analyses rarely stay within a single data model. You may need to:
• Clip a raster to the boundary of a study area (vector polygon) to reduce processing.
• Compute zonal statistics: sum, mean, maximum of raster values within each
polygon of a vector layer (e.g., average elevation per watershed, total rainfall per
administrative district).
• Rasterize vector features into a grid (e.g., burn roads into a raster for cost-distance
analysis).
These operations bridge the two worlds. They rely on the fact that both layers share a
common CRS, and that we can spatially relate pixel indices to coordinates.
Core functions:
• Manual indexing – as we did in the exercise, using [Link]() and numpy array
indexing for maximum control.
Imagine you are assessing flood risk. You have a DEM raster and a vector layer of building
footprints. To determine which buildings are at risk, you need the average elevation within
each footprint. Zonal statistics will give you that table directly. Masking can exclude areas
outside a watershed boundary, focusing a complex model on the region of interest. Without
these tools, you would write nested loops over polygons and pixels—slow and error-prone.
We will use the mask function from [Link]. It requires a dataset (open for reading), a
geometry (Shapely object or GeoJSON-like dict), and optionally crop=True to also reduce
the raster extent.
python
import rasterio
out_meta = [Link]()
out_meta.update({
"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform
})
[Link](out_image)
• [study_area] – a list of geometries; the raster will be masked to the union of all.
• crop=True – reduces the output extent to the bounding box of the geometry, saving
space.
• out_image is a 3D numpy array (bands, rows, cols) with masked pixels set to the
raster’s nodata value (or 0 if not defined). If you want a numpy masked array,
use masked=True in [Link]() earlier.
The rasterstats package (not always installed by default) provides a convenient high-level
interface. Install it with pip install rasterstats.
python
But you can also achieve the same manually using mask and numpy operations, which
gives you more control. We’ll practice both.
For educational clarity, we’ll implement a simple zonal mean extraction using mask and
then compute the mean of the resulting array.
python
def zonal_mean_vectorized(polygon, raster_path):
data = out_image[0]
4. Apply functions row-wise if the number of zones is moderate; for large zones
(>1000), consider using rasterstats which uses optimized routines.
5. Write intermediate rasters if you will repeat the extraction many times with
different polygons (pre-clip once).
6. Be aware of partial pixels – the mask function creates a boolean mask per pixel; a
pixel is included if its center is within the polygon (depending
on all_touched parameter). The default is all_touched=False (center point rule).
python
import numpy as np
# Load watersheds
watersheds = gpd.read_file('[Link]')
slope_crs = [Link]('[Link]').crs
watersheds = watersheds.to_crs(slope_crs)
data = out_img[0]
print(watersheds[['name', 'avg_slope']])
Problem Statement:
Using the DEM and slope rasters you created in the previous exercise
([Link] and [Link]), and a new vector layer of circular survey plots, you will perform
masking, point extraction, and zonal statistics.
3. Point extraction: For the centroid of each plot (the center point), extract the
elevation and slope from the original rasters. Print: plot name, elevation, slope.
4. Zonal statistics (manual): For each plot polygon, compute the mean
elevation and mean slope using the mask function and numpy. Add these as
columns mean_elev and mean_slope to the plots GeoDataFrame. Print the table
(without geometry) showing plot_id, mean_elev, mean_slope.
5. Export the enriched plots GeoDataFrame (with the new mean columns) to a
GeoPackage file plots_stats.gpkg.
Requirements:
Why this exercise: You will integrate the day’s skills: masking, point extraction, and
polygon zonal statistics. This workflow is the basis for ecological surveys, environmental
sampling, and precision agriculture.
Self-check:
• The masked DEM should be a small array covering only the three circles.
• Centroids: compute expected elevation from the sine-cosine formula we used. For
Plot 1 centroid (500500, 4200500), convert to pixel coordinates:
o Elevation formula at (col, row) = (16.67, 83.33) will be around 1000 + 500 *
sin(2π*16.67/100)*cos(2π*83.33/100). You can approximate and check
against script output.
• Mean values will be averages over the circle areas, which will differ from point
values.
Proceed to write the script. In the next lesson, we will move to more advanced raster
operations: reprojection, resampling, and windowed processing for large [Link] will
now examine the solution to the raster-vector integration exercise, consolidating your
ability to mask, extract, and compute zonal statistics—skills that form the bedrock of
environmental analysis. After that, we move into the next critical operation: raster
reprojection and resampling. You will learn how to warp a raster from one CRS to another,
align grids, and change pixel resolution, all while preserving data integrity.
python
#!/usr/bin/env python3
"""
"""
import numpy as np
import rasterio
def main():
centers = [
plots = [Link](
crs="EPSG:32610"
union_geom = plots.unary_union
out_meta = [Link]()
out_meta.update({
"height": dem_masked.shape[1],
"width": dem_masked.shape[2],
"transform": out_transform
})
[Link](dem_masked)
cent = [Link]
r, c = [Link](cent.x, cent.y)
# Check bounds
elev = [Link](1)[r, c]
slp = [Link](1)[r, c]
else:
print("\nZonal statistics:")
print("\nExported plots_stats.gpkg")
if __name__ == "__main__":
main()
Key Points
• Point extraction: Uses the centroid of each plot; we open both rasters together with
a compound with statement.
In an ideal world, all your raster data would share the same CRS, pixel size, and extent.
Reality is messier: a DEM might be in a local UTM projection while a satellite image is in a
different UTM zone or in geographic coordinates. To perform arithmetic between them (e.g.,
computing vegetation index from Landsat bands that are already orthorectified but in
WGS84), you must warp one raster to match the other’s grid.
Reprojection changes the CRS of a raster. Resampling changes the pixel size and/or the
grid alignment. These operations are combined when you warp a raster to a target CRS and
a target resolution. The process involves:
1. Computing the output grid (extent and cell positions) based on the target CRS and
transform.
2. For each output pixel, determining which input pixel(s) it corresponds to.
3. Applying a resampling method to compute the output value from the input values.
• nearest (default) – takes the value of the closest input pixel; good for categorical
data (land cover).
• bilinear – linear interpolation between 4 nearest pixels; good for continuous data
(elevation, temperature).
• cubic – cubic convolution using 16 pixels; slightly sharper but can overshoot.
• average – averages pixels that fall within the output pixel; good for downsampling to
avoid aliasing.
Why this matters: Reprojecting a raster allows you to align it with vector data or another
raster for pixel-by-pixel calculations. Resampling lets you reduce resolution for faster
processing or increase resolution (though that does not add information).
I once worked with a climate model that output global precipitation data on a 0.5° grid in
WGS84, and a hydrological model that required input on a 1-km Albers Equal-Area grid.
Without reprojection and resampling, the two could never talk.
The [Link] function is the Swiss Army knife that handles this translation
efficiently, in memory or on disk.
The function can be used to copy data from one open raster to another, or from an
in-memory numpy array to a file. The basic signature:
python
import rasterio
# Update metadata
kwargs = [Link]()
[Link]({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
reproject(
source=[Link](src, i),
destination=[Link](dst, i),
src_transform=[Link],
src_crs=[Link],
dst_transform=transform,
dst_crs=dst_crs,
resampling=[Link]
• For land-cover classifications, nearest preserves the original class values without
mixing.
Often you want to warp a raster to exactly match the grid of another raster (so they can be
stacked or compared pixel-wise). You can read the profile of the reference raster and use
its transform and dimensions directly.
python
dst_transform = [Link]
dst_crs = [Link]
dst_width = [Link]
dst_height = [Link]
6. Resample before other operations – it’s often efficient to resample to the target
resolution early in a pipeline.
Assume we have a DEM in UTM zone 10N (California) and we need it in zone 11N for an
overlapping study area.
python
import rasterio
dst_crs = 'EPSG:32611'
kwargs = [Link]()
reproject(
source=[Link](src, 1),
destination=[Link](dst, 1),
src_transform=[Link],
src_crs=[Link],
dst_transform=transform,
dst_crs=dst_crs,
resampling=[Link]
Problem Statement:
You have two synthetic rasters representing the same area but in different CRSs and
resolutions. Your task is to align them so they can be compared.
1. Create the first raster (ndvi_utm.tif) exactly as the NDVI result from the NDVI
exercise (or create a new one):
o Write as float32.
3. Reproject and resample the geographic NDVI raster to match the UTM NDVI
raster’s grid (CRS EPSG:32610, same transform and dimensions as ndvi_utm.tif).
Save the result as ndvi_geo_warped.tif. Use bilinear resampling.
4. Compare the two rasters pixel-wise: read both the original UTM NDVI and the
warped geographic NDVI, compute the difference array (diff = ndvi_utm -
ndvi_geo_warped), and print the maximum absolute difference. Also print the mean
absolute difference.
5. Extract values at three sample points (in UTM coordinates):
Requirements:
Why this exercise: You’ll practice the complete warping workflow: creating rasters with
different CRSs, aligning them, and evaluating the impact of resampling.
Self-check:
• The two original rasters have different resolutions and CRSs, so the warped one will
be interpolated, causing differences.
Proceed to write the script. In the next lesson, we’ll move to the final vector–raster
integration: point cloud processing and 3D [Link] will now review the solution to the
raster warping exercise, which demonstrates the crucial ability to align datasets from
different coordinate systems. Following that, we will step into the third dimension—
literally—by working with point cloud data. LiDAR and other 3D surveying techniques
produce millions of points with X, Y, Z, and often attributes like classification. Python gives
you the tools to read, filter, and analyze these massive datasets.
python
#!/usr/bin/env python3
"""
"""
import numpy as np
import rasterio
def main():
rng1 = [Link].default_rng(0)
rng2 = [Link].default_rng(1)
pixel_size = 30.0
profile_utm = {
[Link](ndvi_utm, 1)
print("Created ndvi_utm.tif")
profile_geo = {
[Link](ndvi_geo, 1)
print("Created ndvi_geo.tif")
dst_crs = 'EPSG:32610'
# However, to exactly match the UTM raster's transform and size, we'll force them.
transform = transform_utm
warped_width = width
warped_height = height
profile_warped = [Link]()
profile_warped.update({
'crs': dst_crs,
'transform': transform,
'width': warped_width,
'height': warped_height,
'compress': 'lzw'
})
reproject(
source=[Link](src, 1),
destination=[Link](dst, 1),
src_transform=[Link],
src_crs=[Link],
dst_transform=transform,
dst_crs=dst_crs,
resampling=[Link]
print("Created ndvi_geo_warped.tif")
# 4. Compare the two rasters
utm_arr = [Link](1)
warped_arr = [Link](1)
max_abs_diff = [Link]([Link](diff))
mean_abs_diff = [Link]([Link](diff))
print("\nPoint extraction:")
r, c = [Link](x, y)
val_utm = [Link](1)[r, c]
val_warp = [Link](1)[r, c]
if __name__ == "__main__":
main()
Key Observations
• Forced matching: Instead of using the default transform computed
by calculate_default_transform, we explicitly set the output transform and
dimensions to match the reference UTM raster exactly. This is crucial when you
need a pixel-perfect stack.
• Differences arise because the original geographic raster has a different pixel size
(0.001° ≈ 111 m) and the resampling (bilinear) interpolates values onto the new
30-m grid. The UTM raster was generated independently with different random
noise, so the differences reflect both the noise difference and the resampling
smoothing.
A point cloud is a collection of 3D points (X, Y, Z) that represent the surface of objects.
LiDAR (Light Detection and Ranging) systems on aircraft, drones, or terrestrial scanners
emit laser pulses and record the returns, generating millions of points with attributes:
Point clouds are used for creating digital elevation models (DEM), digital surface models
(DSM), canopy height models, building footprints, and powerline inspection.
The standard file format for LiDAR data is LAS (or its compressed version LAZ). In Python,
the library laspy provides efficient reading, writing, and manipulation of LAS/LAZ files. It
supports the LAS 1.2–1.4 specifications and can handle large files via memory-efficient
iterators.
Raster grids are uniform; point clouds are irregular. Many advanced geospatial tasks
require extracting information directly from the raw points before gridding. For example,
ground filtering (separating ground from non-ground points) is essential for producing
bare-earth DEMs. Canopy metrics (e.g., height percentiles) are computed from vegetation
points. By learning laspy, you can bypass pre-processed products and work with the raw
survey data.
python
import laspy
print([Link].point_count)
print([Link].point_format)
points = [Link]()
# Access arrays
x = points.x
y = points.y
z = points.z
classification = [Link]
laspy returns a LasData object, where each attribute (x, y, z, intensity, etc.) is a numpy
array. You can slice, filter, and analyze them using standard numpy operations.
4. Memory-Efficient Iteration
Loading millions of points into memory might exceed RAM. Use las.chunk_iterator to
process in chunks.
python
y = chunk.y
z = chunk.z
# process chunk
Using numpy boolean indexing, you can quickly select ground points, high-vegetation, etc.
python
ground_points = points[ground]
non_ground = points[~ground]
• 0 – Never classified
• 1 – Unassigned
• 2 – Ground
• 3 – Low Vegetation
• 4 – Medium Vegetation
• 5 – High Vegetation
• 6 – Building
• 8 – Reserved
• 9 – Water
• etc.
The simplest way is to use rasterio with a grid. We can iterate over ground points and bin
them into cells, averaging Z. Or use laspy’s create_dem or external libraries
like pdal (Python PDAL is also powerful). We’ll implement a basic binning for educational
purposes.
7. “How to Write” Rules for Point Clouds
1. Check the point format – LAS files store different attributes depending on the
format version; laspy handles this.
2. Use iterators for large files – never assume the whole file fits in memory.
Assume we have a LAS with ground and vegetation classes. We can compute the height
above ground for each point and then rasterize canopy height.
python
import numpy as np
import laspy
import rasterio
points = [Link]()
ground = points[[Link] == 2]
veg = points[[Link]([3,4,5])]
# Create a simple ground model by interpolating ground points (we'll just use min Z per grid
cell)
# Define grid
# Bin ground points: store minimum Z per cell (simplistic ground model)
ground_grid[row, col] = z
chm = np.zeros_like(ground_grid)
if 0 <= row < height and 0 <= col < width and ground_grid[row, col] != [Link]:
[Link]([Link]('float32'), 1)
(Note: In practice, you would use proper interpolation like Delaunay triangulation or inverse
distance weighting for the ground surface, and use libraries like [Link] or pyinterp.)
Create a new file point_cloud_basics.py. We will generate a synthetic point cloud (since
real LAS files may not be available) to practice reading, filtering, and creating a simple
DEM.
Problem Statement:
1. Generate a synthetic point cloud using numpy to simulate a flat ground with some
random noise and a few “buildings” (vertical pillars):
o Populate x, y, z, classification.
o Set the CRS (choose EPSG:32610 for UTM). Use [Link](...) to
add a WKT CRS if you want, but for simplicity, just note that the coordinates
are in that system.
o For each cell, compute the mean Z of the ground points falling into that cell.
Leave empty cells as [Link].
o Write the DEM to dem_synthetic.tif with CRS EPSG:32610. Use mean Z as the
cell value.
o Grid all points (ground + buildings) using maximum Z per cell (to capture
building tops).
o Write to dsm_synthetic.tif.
Requirements:
Why this exercise: You will simulate the full LiDAR workflow: point generation, LAS
writing/reading, classification filtering, and rasterization to DEM/DSM. This is the
foundation of 3D geospatial analysis.
Self-check:
• The synthetic LAS should open correctly in point cloud viewers (optional).
• DEM should show a nearly flat surface (mean ~0) with holes where buildings are (if
you only use ground points). Actually, ground points exist everywhere, so DEM will
cover the whole area with near-zero values.
Proceed to write the script. In the following lesson, we will explore web mapping and
deployment—turning your analysis into interactive [Link] will now review the
solution to the point cloud exercise. This lesson ties together synthetic data generation,
LAS file creation, classification filtering, and rasterization—a complete LiDAR simulation.
After that, we will shift from the analytical back end to the visual front end: interactive web
mapping. You will learn to turn your GeoDataFrames and rasters into dynamic, shareable
maps using Python alone.
python
#!/usr/bin/env python3
"""
"""
import numpy as np
import laspy
import rasterio
def main():
n_ground = 10_000
n_bldg1 = 1000
n_bldg2 = 1000
# Concatenate all
las_file.x = x
las_file.y = y
las_file.z = z
las_file.classification = classification
las_file.write("[Link]")
with [Link]("[Link]") as f:
las = [Link]()
ground_mask = [Link] == 2
building_mask = [Link] == 6
gx = las.x[ground_mask]
gy = las.y[ground_mask]
gz = las.z[ground_mask]
valid = (row >= 0) & (row < height) & (col >= 0) & (col < width)
# Write DEM
profile_dem = {
[Link]([Link]('float32'), 1)
print("Created dem_synthetic.tif")
valid_all = (row_all >= 0) & (row_all < height) & (col_all >= 0) & (col_all < width)
[Link]([Link]('float32'), 1)
print("Created dsm_synthetic.tif")
if __name__ == "__main__":
main()
Key Points
• DSM uses all points, taking the maximum Z per cell, so building tops are captured.
LESSON 7.1 — Interactive Web Mapping with folium and leafmap
Static maps printed on paper or saved as PNGs are useful for reports, but the modern
geospatial workflow often demands interactive, zoomable, clickable maps that can be
shared via a URL. Python can generate such maps using libraries that write HTML, CSS, and
JavaScript behind the scenes, leveraging the [Link] mapping library.
• folium is a lightweight library that creates Leaflet maps from Python. You add
GeoDataFrames as layers, customize pop-ups, add markers, and save as a
standalone HTML file or integrate into Jupyter notebooks.
• leafmap builds on top of folium and ipyleaflet, adding support for raster layers, COG
(Cloud Optimized GeoTIFF) streaming, and splitting maps for comparison. It’s a
more feature-rich tool for raster-vector combination.
Today, we focus on folium as the entry point, because it’s simple and directly linked to the
GeoDataFrame world you know.
You’ve performed complex analysis; now you must present it to stakeholders who may not
have Python. An interactive map lets them explore the data, click on features to see
attributes, and understand spatial patterns intuitively. This is often the final stage of a
project, and Python handles it seamlessly.
python
import folium
[Link]('[Link]')
3.2 Adding a GeoDataFrame (Points)
python
gdf = gpd.read_file('[Link]')
[Link](
gdf,
name='Stations',
tooltip=[Link](fields=['name', 'elevation'])
).add_to(m)
python
[Link](
geo_data=gdf,
data=gdf,
columns=['id', 'population'],
key_on='[Link]',
fill_color='YlOrRd',
legend_name='Population'
).add_to(m)
python
[Link]().add_to(m)
4. leafmap for Rasters
leafmap can directly add a Cloud Optimized GeoTIFF from a URL or local file
using add_raster(). This renders the raster on the map using a tile service.
python
import leafmap
m = [Link]()
m.to_html('raster_map.html')
1. Keep data layers reasonably sized – large GeoJSONs can slow down browsers.
Simplify geometries if needed.
5. For rasters, use COG format – leafmap can stream them efficiently.
Problem Statement:
Using the plots_stats.gpkg from the raster-vector lab (the survey plots with mean elevation
and slope), and the dem_synthetic.tif from the point cloud lesson, create an interactive
map that:
1. Centers on the extent of the plots (approximate coordinates from the previous
exercise). You can set location to [4201500, 501000] (UTM coordinates? Wait, folium
expects lat/lon. Reproject the plots to EPSG:4326 before mapping).
o Use [Link] with a style function that fills polygons with a color
based on mean_elev (e.g., green for low, red for high). Define a simple color
ramp manually: if mean_elev < ? (you decide thresholds based on your data
range). Or use [Link] with mean_elev column.
Requirements:
• Ensure the raster CRS is also EPSG:4326; if not, reproject the raster or use leafmap’s
ability to handle other CRSs (leafmap can reproject on the fly).
Why this exercise: You will bridge analysis and presentation, learning to export your
spatial results to an interactive medium. This is often the final deliverable in real-world
projects.
Self-check:
• The HTML file should open in a web browser showing an interactive map.
• DEM layer can be toggled; it may need to be warped to EPSG:3857 (leafmap can
handle that automatically with COG).
Proceed to write the script. In the final lesson, we’ll cover packaging your code into
reusable modules and a concluding [Link] have reached the final stage of this
foundational journey. The previous lesson taught you to present data; now, we learn to
package and share the tools themselves. A script is a tool for yourself;
a module or package is a tool for others (and your future self). This lesson covers the
structure, testing, and deployment of your geospatial Python code into something
reusable, robust, and professional.
A single .py file is a script or a module. When a project grows to multiple related modules
(e.g., [Link], [Link], [Link], [Link]), you organize them into a package—a
directory containing an __init__.py file and the modules. This structure allows clean
imports:
python
• Tests – a tests/ directory with scripts that verify the correctness of your functions.
Why this matters: You’ve written many functions across exercises. If you needed to
reuse haversine_distance in a new project, you’d have to copy-paste it. That leads to
divergence and bugs. By creating your own geo_utils package, you install it once and import
it anywhere. Professionally, all geospatial libraries you use (GeoPandas, rasterio, shapely)
are packages structured this way.
When you find yourself emailing a script to a colleague, you are distributing code in its most
fragile form. They may have slightly different paths, missing dependencies, or Python
versions. A package with a clear environment file ([Link]) and tests ensures that
your tool works as intended on any machine. If you later need to publish your work (e.g., on
GitHub, PyPI), this structure is mandatory.
3. Package Structure
geo_utils/
├── geo_utils/
│ ├── __init__.py
│ ├── [Link]
│ ├── crs_helpers.py
│ └── raster_utils.py
├── tests/
│ ├── test_distances.py
│ └── test_crs_helpers.py
├── [Link]
├── [Link]
└── [Link]
• __init__.py can be empty or import key functions to expose at the package level.
• tests/ contains files named test_*.py. Each file contains functions starting
with test_.
pytest is a testing framework. You write a function that calls your code and asserts an
expected result. If the assertion fails, the test fails.
python
# tests/test_distances.py
import math
def test_haversine_same_point():
dist = haversine_distance(0, 0, 0, 0)
def test_haversine_known_distance():
Run tests from the terminal: pytest tests/. You’ll see a green dot for each passing test.
CI services (like GitHub Actions) automatically run your tests every time you push code. A
configuration file (.github/workflows/[Link]) sets up a Python environment, installs
dependencies, and runs pytest. This ensures that contributions don’t break existing
functionality.
2. One function per major task – keep functions small and focused.
6. Include tests for the most critical functions (especially those involving coordinate
transformations and geometry operations).
8. Add a [Link] explaining what the package does and how to install it.
Create a package from the geospatial functions you’ve developed during this course. This
will be your capstone—a personal toolkit you can expand forever.
Task:
1. Create a directory structure as described above for a package
called geocraft (choose any name). You can do this in a new
folder geocraft_package/.
o __init__.py – import the key functions so they are accessible as from geocraft
import haversine_distance.
toml
[build-system]
requires = ["setuptools"]
build-backend = "setuptools.build_meta"
[project]
name = "geocraft"
version = "0.1.0"
requires-python = ">=3.9"
dependencies = ["numpy", "rasterio", "geopandas", "shapely"]
5. Create a [Link] explaining how to install and use the package (1–2
paragraphs).
6. (Optional) If you have Git installed, initialize a repository and make a commit.
7. Run pip install -e . from the root directory (where [Link] is) to install your
package in development mode. Then open a Python interpreter and test:
python
Why this final exercise: You’re not just writing code; you’re building a tool. This structure is
exactly how professionals deliver geospatial capabilities. It consolidates everything you’ve
learned: functions, CRS, files, testing, and documentation.
Self-check:
You have traveled from the very first print("Hello") to building your own geospatial Python
package. You now understand variables, loops, functions, file I/O, Shapely geometries,
GeoPandas spatial joins, raster I/O with rasterio, reprojection, point cloud processing, and
even interactive mapping. Most importantly, you have a mental model of how spatial data
is represented in Python and the rules of writing clean, professional code.
You are a geospatial analyst for a municipal government. A river flows through the city from
west to east. Recent climate projections indicate that extreme rainfall could raise the
river’s water level by 5 meters above its normal bank. Your task is to:
1. Delimit the flood inundation zone—the area that would be submerged under a
5-meter rise.
4. Produce an interactive map showing the river, flood zone, and at-risk buildings with
attributes.
This is not a toy problem. It integrates raster elevation models, vector building
footprints, spatial predicates, zonal statistics, and web map visualization—a complete
geospatial pipeline that mirrors what city planners and emergency services use.
We require three core datasets, all in the same projected CRS (metric units, e.g., UTM):
• DEM (Digital Elevation Model): A raster where each pixel’s value is the ground
elevation in meters. Resolution e.g. 5 m.
• Building footprints: A polygon vector layer with an id and perhaps floors or height.
• Population density raster (optional for advanced step): A raster of people per pixel,
same grid as DEM or different.
For educational purposes, we will synthesize these datasets in-memory. This guarantees
you can run the entire pipeline without external data, and you see how to create test data
for any geospatial problem.
o DEM: a gentle slope with some hills and a river valley carved into it.
o Buffer the river centerline by a distance (e.g., 100 m) to capture the river’s
floodplain.
o Mask this buffered polygon to only include cells where DEM <
threshold_elevation (threshold = river surface elevation + 5 m).
o Zonal sum of population density raster within each at-risk building polygon
(or simply sum over the whole flood zone).
o Add layers: river, flood zone (translucent blue), buildings (color by at-risk),
population overlay.
6. Export results
#!/usr/bin/env python3
"""
Demonstrates raster-vector integration, spatial joins, zonal stats, and web mapping.
"""
import numpy as np
import rasterio
import folium
• UTM CRS (e.g., EPSG:32610, WGS84/UTM zone 10N). Center coordinates roughly
(500000, 4200000) east/north.
How to write: Use numpy to create an array, then write a GeoTIFF. We’ll set the
transformation explicitly.
python
# DEM parameters
crs = "EPSG:32610"
# Create a geometry of the river line, then use `geopandas` with a buffer to mask.
# We'll rasterize the line into a binary mask, then compute distance transform.
all_touched=True)
# from_origin gives positive pixel_width and negative pixel_height if we use standard north-
up? Actually, from_origin takes west, north, xsize, ysize (positive ysize means south-down?
We'll use Affine explicitly to avoid confusion.
# Correction: river_raster is 1 where line exists, 0 elsewhere. We need distance from pixels
to the line.
elevation -= river_depth
# Save DEM
dem_profile = {
'compress': 'lzw'
[Link]([Link]('float32'), 1)
We generate random rectangular buildings, slightly clustered near the river, to simulate a
small town.
python
[Link](42)
n_buildings = 150
building_geoms = []
building_ids = []
for i in range(n_buildings):
center = [Link](town_centers)
building_geoms.append(rect)
building_ids.append(f"B{i:04d}")
buildings_gdf.to_file('[Link]', driver='GeoJSON')
Why: We used [Link] to orient buildings randomly; this mimics real cities
where buildings align to streets. The floors column will be used to estimate population.
We create a population density raster on the same grid as the DEM. The density will be
higher near the town centers, with a kernel density-like pattern. We can use a simple
Gaussian blur of a point density.
python
# Smooth
pop_profile = dem_profile.copy()
pop_profile.update(dtype='float32')
[Link](pop_raster.astype('float32'), 1)
Why: We use a kernel density approach, placing Gaussian blobs at town centers. The
scaling ensures a plausible total population. In a real project, you would use dasymetric
mapping from census data.
We simulate a 5-meter water rise. The river’s normal water surface elevation can be
approximated by the DEM values at points along the river line. We’ll sample the DEM at the
river start and end, take the mean as the base water level, then add 5 m.
python
# Sample DEM at river start and end to get base water elevation
base_elev_start = next([Link]([river_start]))[0]
base_elev_end = next([Link]([river_end]))[0]
base_water_level = (base_elev_start + base_elev_end) / 2.0
river_buffer = river_line.buffer(200)
flood_zone_gdf.to_file('flood_zone.geojson', driver='GeoJSON')
Why:
• Sampling the DEM at river ends gives a realistic water surface, assuming the river
slopes naturally.
• We buffer the river first to limit the analysis extent; this avoids processing the whole
DEM.
• The threshold condition <= flood_level identifies all pixels that would be underwater.
python
if 'index_right' in at_risk.columns:
at_risk = at_risk.drop(columns='index_right')
at_risk.to_file('buildings_at_risk.geojson', driver='GeoJSON')
Why intersects not within: A building partially inside the flood zone is still at
risk. intersects captures any overlap, including edges.
We’ll compute the total population within the flood zone using the population density
raster.
python
pop_array = pop_clip[0]
total_pop_at_risk = [Link](pop_array)
We’ll project all layers to EPSG:4326 for web display. Folium works directly with
GeoDataFrames if they are in lat/lon.
python
# Reproject to WGS84
buildings_gdf_wgs = buildings_gdf.to_crs('EPSG:4326')
at_risk_wgs = at_risk.to_crs('EPSG:4326')
flood_zone_wgs = flood_zone_gdf.to_crs('EPSG:4326')
# Add river
[Link]().add_to(m)
[Link]('flood_map.html')
Why: We separate layers into different GeoJson objects so the user can toggle them. The
tooltip shows building details.
• Raster resolution: 5 m for a 2000 m extent yields 160,000 pixels—tiny. Real cities
might have 0.5 m resolution over tens of kilometers; use windowed processing
and [Link] with crop=True.
• CRS alignment: We ensured all data is in the same UTM CRS before any spatial
operation. Reprojecting rasters on the fly can be slow; do it once at the beginning.
• Zonal stats for many polygons: For thousands of buildings, loop over each building
and mask separately, or use rasterstats library which is optimized. Our single mask
of the entire flood zone is efficient.
2. CRS is the first thing to check – every file has a .crs, and we assert they match.
6. Interactive maps are the final deliverable; always include layer control and
informative tooltips.
I have now walked you through the entire pipeline. The complete script, when assembled,
runs end-to-end and produces all outputs. You can find the consolidated version below (I
will provide it as a single block for your convenience, but you should write it incrementally).
• Include a DEM with levees (raised linear features) that protect some areas, and see
how the flood zone changes.
• Calculate damages by assigning a monetary value per building floor area and
summing over at-risk buildings.
• Optimize: time the script and replace any slow loops with vectorized or parallel
operations.
• Deploy: wrap the core logic in a function and call it from a Flask web API.
This is how you move from a script to a solution. You now have the foundation to tackle any
geospatial problem with confidence.
Imagine a fleet of 10,000 delivery vehicles, each reporting its position every 10 seconds for
a month. The resulting dataset contains over 2.5 billion GPS points. Your task:
• Read this massive dataset (stored as partitioned Parquet files) without loading it all
into memory.
• Compute per-vehicle statistics: total distance traveled, average speed, idle time.
• Snap each cleaned point to the nearest road segment (using a road network
shapefile with 2 million edges) via spatial join.
This problem cannot be solved with ordinary GeoPandas or rasterio because the data
exceeds RAM. We need lazy, out-of-core, and parallel execution.
Dask extends the familiar Pandas/NumPy APIs to parallel, distributed computing. It splits
large datasets into chunks—small Pandas DataFrames or NumPy arrays that fit in
memory—and processes them in parallel, either on a single machine (using multiple cores)
or on a cluster. Key concepts:
• Dask Array: a large parallel NumPy array. Used by rioxarray to handle huge rasters.
• Xarray: a library for labeled multi-dimensional arrays (like NetCDF). With rioxarray, it
wraps rasterio and can use Dask arrays as backends, enabling lazy raster
processing.
Why this matters: You can write code that looks exactly like GeoPandas/rasterio but
processes data that is 100 times larger than your available memory, using all CPU cores.
The transition from single-machine to big-data geospatial is smoother than you think.
1. Generate synthetic partitioned Parquet files with GPS data, using Dask to write them
lazily.
We’ll run everything on a local machine with a multi-core scheduler, but the code can be
deployed to a distributed cluster (like Coiled or a SLURM cluster) with zero changes.
4. Step-by-Step Implementation
We’ll create a new script large_trajectory_analysis.py. I’ll explain each chunk of code and
the design decisions behind it.
python
import numpy as np
import pandas as pd
import [Link] as dd
import dask_geopandas as dask_gpd
import [Link] as da
import [Link]
import rasterio
import rioxarray
python
client = Client(cluster)
print(client)
This starts a local cluster with 4 workers, each with 2 threads and 2 GB RAM limit—
simulating a distributed environment. We’ll use the Client to monitor tasks.
We’ll create 10 million points spread across a city (say, 100 km² area), with
a vehicle_id (1000 vehicles) and timestamp. We’ll write them as a Parquet dataset
partitioned by vehicle_id to allow efficient per-vehicle queries.
python
# Parameters
n_vehicles = 1000
# Generate synthetic data lazily with Dask's bag or delayed, but for simplicity we'll use
Pandas on chunks
# Since we're demonstrating Dask, we'll create a Dask DataFrame from delayed Pandas
DataFrames
import [Link] as db
def create_partition(vehicle_ids_chunk):
dfs = []
df = [Link]({
'vehicle_id': vid,
'timestamp': times,
'lon': lons,
'lat': lats,
'speed': speed
})
[Link](df)
return [Link](dfs)
print([Link]())
Why this pattern: We used [Link] to wrap the Pandas creation function,
then dd.from_delayed to build a Dask DataFrame. The .persist() triggers computation and
keeps the data in distributed memory for subsequent operations.
python
def create_geometry(df):
'vehicle_id': [Link](dtype='int'),
'timestamp': [Link](dtype='datetime64[ns]'),
'lon': [Link](dtype='float'),
'lat': [Link](dtype='float'),
'speed': [Link](dtype='float'),
'geometry': [Link](crs='EPSG:4326')
}))
ddf_gpd = ddf_gpd.persist()
Why: map_partitions applies a function to each Dask partition, returning a Dask
DataFrame. We must provide a meta describing the output schema to avoid inference
costs.
Lazily filter points within a bounding box (city limits) and remove anomalies (speed > 150
km/h).
python
mask = ddf_gpd_proj.within(bbox)
filtered = ddf_gpd_proj[mask]
# Speed filter
filtered = [Link]()
Note: Dask GeoPandas partitions data spatially, so within can push down filters to relevant
partitions (using spatial indexing behind the scenes). This makes it much faster than brute
force.
We’ll use a small synthetic road network for demo, but the same code works for a real
shapefile. We’ll generate it as a single GeoDataFrame and then convert to Dask
GeoDataFrame with spatial partitioning.
python
road_lines = []
for _ in range(2000):
x1, y1 = [Link](550000, 560000), [Link](4180000, 4190000)
roads_gdf['road_id'] = range(len(roads_gdf))
roads_ddf = roads_ddf.persist()
We can’t do a direct sjoin because we want the nearest road, not just intersecting. Dask
GeoPandas doesn’t yet have a distributed nearest-join built in, but we can use a cross join
with a spatial distance filter and then pick the minimum per point. This is feasible
because Dask can handle large cross joins when both sides are partitioned.
python
joined = [Link]()
Why sjoin_nearest: It uses a spatial index on each partition and broadcasts the right
GeoDataFrame to each partition, performing efficient nearest queries. This scales well.
python
raster_profile = {
[Link](traffic_data, 1)
import rioxarray
print(rds)
Now we need to extract raster values at the points in filtered. We’ll collect the point
coordinates (as a Dask Series of x,y) and then use rasterio’s sample or map function.
However, doing this lazily for millions of points can be tricky. We can:
• Simpler: use [Link](x=..., y=...) which works with Dask, but requires x and y as
arrays. We can extract [Link].x and .y as Dask Series and convert to Dask
Array.
python
x_coords = [Link].x.to_dask_array()
y_coords = [Link].y.to_dask_array()
Why interp: Xarray’s interp method is lazy and uses Dask, so it scales. We must ensure
coordinates are chunked similarly to the raster.
python
# Add sampled traffic value back to filtered GeoDataFrame? We had to compute sampled,
so we can add as column
filtered['traffic'] = sampled
# To compute distance, we need consecutive points per vehicle; we'll sort by time and
compute using UTM coordinates
# This is a typical window function. Dask can do groupby-apply with a custom function.
def compute_vehicle_stats(df):
df = df.sort_values('timestamp')
shifted = df[['geometry']].shift()
avg_speed = [Link]()
avg_traffic = [Link]()
return [Link]({'total_dist_km': total_dist_km, 'avg_speed': avg_speed, 'avg_traffic':
avg_traffic})
}).compute()
print(vehicle_stats.head())
Why groupby-apply: Dask executes this in parallel: each partition groups by vehicle, but a
vehicle’s points may span multiple partitions. Dask’s groupby with apply handles the
shuffling automatically. The meta helps define the output schema.
We’ll aggregate the cleaned points to a grid and create a heatmap (using Datashader or just
Dask array). Finally, save the results.
python
x = [Link].x.to_dask_array()
y = [Link].y.to_dask_array()
# Define grid
density = [Link]()
dtype='float32', crs='EPSG:32610',
transform=[Link].from_bounds(x_bins[0], y_bins[0], x_bins[-1],
y_bins[-1],
[Link]([Link]('float32'), 1)
vehicle_stats.to_csv('vehicle_stats.csv')
print("Analysis complete.")
1. Start small, then scale – prototype with Pandas on a sample, then switch to Dask
with minimal changes (many Pandas methods are identical in Dask).
3. Lazy until necessary – call .compute() only at the end, or to inspect intermediate
results.
4. Memory limits – set worker memory limits to avoid out-of-memory crashes. Dask
can spill to disk if needed.
5. Raster chunking – align chunks with the access pattern; for point extraction, chunk
spatially to avoid reading the whole raster.
Now you will modify the above script to handle a fleet of 50,000 vehicles (500 million
points). Don’t generate the full dataset—just modify the code to use more partitions and a
smaller memory limit. Also:
• Replace the synthetic road network with a real OpenStreetMap extract (download a
small .[Link] file, use osmnx to convert to GeoDataFrame, then to Dask
GeoDataFrame).
• Instead of interp, use [Link] in a map_partitions call to extract raster
values more efficiently.
Why this exercise: It forces you to think about memory, partitioning, and lazy operations at
scale—skills that differentiate a journeyman from a master.
1. Real-world road network from OSM, loaded efficiently and converted to a spatially
partitioned Dask GeoDataFrame ready for nearest-neighbor joins.
For a city-scale analysis, downloading the entire planet is wasteful. We use OSMnx, which
can fetch street networks for a defined bounding polygon and return a GeoDataFrame.
Since we need a Dask GeoDataFrame, we'll fetch a moderately large area, convert to a
standard GeoDataFrame, then repartition.
python
import osmnx as ox
# Define the area of interest (same as our GPS data extent, in lat/lon)
# Download street network (all drivable roads) as a graph, then convert to GeoDataFrame
G = ox.graph_from_bbox(north, south, east, west, network_type='drive')
Why OSMnx: It handles the OSM API, parses tags, and simplifies the topology into a clean
edge table. The graph_from_bbox function is memory-efficient and uses a bounding box,
which matches our problem.
For larger areas (e.g., an entire state), you would download pre-extracted OSM extracts
(.[Link]) and use pyrosm or osm2pgsql with a PostGIS database, but OSMnx is ideal for
city scale.
We now have edges_gdf with thousands of road segments. To enable fast spatial joins, we
partition it spatially using the Hilbert curve, which groups nearby geometries into the same
partition.
python
roads_ddf = roads_ddf.persist()
python
joined = [Link]()
If we need to keep all GPS points (even those far from roads), we use how='left' and later
filter dist_to_road to identify outliers. The distance_col will contain the actual distance to
the matched road, invaluable for quality control.
The [Link] method we used earlier builds a large Dask array of coordinates, which
can exhaust scheduler memory for hundreds of millions of points. A more memory-friendly
approach is to use [Link] inside each Dask partition: for a chunk of GPS points,
we open the raster (using a shared, read-only handle or by re-opening inside the function)
and extract values only for that chunk. This limits peak memory to the chunk size.
We'll use the traffic density raster created earlier (traffic_density.tif). For best performance,
we ensure it is a Cloud-Optimized GeoTIFF (COG) or at least a tiled GeoTIFF with internal
overviews, but for our scale we can use a regular GeoTIFF with caching.
We'll use rasterio's block caching: when multiple chunks access the same raster, the
operating system's disk cache helps, but we can also explicitly
use [Link] with GDAL_CACHEMAX.
python
import rasterio
def extract_traffic(df, raster_path):
# df is a GeoDataFrame partition
values = list([Link](coords))
df['traffic_density'] = traffic
return df
Why [Link]: It is written in C and processes coordinates in bulk, much faster than a
Python loop with [Link] and array indexing. For millions of points per partition, this is
critical.
We now apply this function to the Dask GeoDataFrame. We need to pass the raster path as
an additional argument, using map_partitions with args or kwargs. We must also supply the
output meta to avoid inference.
python
meta = filtered_gps._meta.copy()
meta['traffic_density'] = [Link](dtype='float64')
enriched = filtered_gps.map_partitions(
extract_traffic,
meta=meta
)
enriched = [Link]()
Memory management: Each partition loads into memory as a Pandas DataFrame, plus the
extracted column. We should ensure partition sizes are small enough (e.g., 100k–500k
points). We can repartition the Dask GeoDataFrame before
extraction: filtered_gps.repartition(partition_size='100MB').
For multi-band rasters (e.g., satellite imagery), [Link] returns a list of lists, one per
band. We can expand them into separate columns using a loop over band indices. This is
straightforward but must be reflected in the meta.
• File handle contention: If many workers open the same file simultaneously, the OS
can handle it with file caching. For thousands of workers, consider using a
distributed file system or a VRT.
• Windowed raster reading: For very large rasters, we can open the raster with a
specific window that covers the partition's bounding box,
using [Link].from_bounds. This reduces the amount of data read.
However, [Link] does not benefit from windowing; it reads individual pixel
locations. For optimization, we can pre-clip the raster to the partition's extent
using mask and then sample, but that would require writing temporary files—often
not worth the overhead.
We'll now outline the complete extended script, integrating OSM roads and the optimized
raster extraction. We assume the earlier steps (GPS generation, filtering) are already done
and we have filtered as a Dask GeoDataFrame.
python
import osmnx as ox
import dask_geopandas as dask_gpd
edges['road_id'] = range(len(edges))
roads_ddf = roads_ddf.persist()
# Nearest join
joined = dask_gpd.sjoin_nearest(
filtered, roads_ddf,
max_distance=50,
how='left',
distance_col='dist_to_road'
joined = [Link]()
# Section: Raster extraction on the enriched points (which now have road info)
vals = [Link](coords)
df['traffic_density'] = traffic
return df
meta = joined._meta.copy()
meta['traffic_density'] = float
joined = [Link]()
def vehicle_stats(df):
df = df.sort_values('timestamp')
shifted = df[['geometry']].shift()
avg_speed = [Link]()
avg_traffic = df.traffic_density.mean()
return [Link]({
'total_dist_km': total_dist_km,
'avg_speed': avg_speed,
'avg_traffic': avg_traffic
})
}).compute()
stats.to_csv('vehicle_stats_with_roads.csv')
• The OSM road network is partitioned spatially, matching the GPS partitions' spatial
distribution.
• sjoin_nearest broadcasts only the necessary road partitions to each GPS partition.
• The groupby-apply for statistics is a shuffle operation that Dask handles efficiently,
spilling to disk if necessary.
1. Align spatial partitioning: Both left and right GeoDataFrames should use the same
spatial partitioning strategy ('hilbert'). This minimizes data movement during joins.
3. Use map_partitions with care: The function must be serializable and should not rely
on global state. Pass file paths as arguments.
4. Set appropriate partition sizes: Aim for 100–200 MB per partition after filtering.
Use repartition() or partition_size in from_geopandas.
5. Monitor with the dashboard: Keep an eye on the Dask dashboard's "Task Stream"
to identify bottlenecks (e.g., large shuffles, uneven partitions).
7. Lazy writes: Use to_parquet and to_csv on Dask DataFrames to write results lazily;
they'll execute as part of the task graph.
Now that you have seen the pieces, I challenge you to modify the script to:
• Use a pre-downloaded OSM extract for a large state (e.g., California) from
Geofabrik. Load it with pyrosm to get a filtered road network (highways + primary)
and convert to Dask GeoDataFrame.
• Simulate 50 million GPS points from 10,000 vehicles across the state, partitioned by
a spatial grid.
• Use a real traffic density raster (maybe a population density grid) from a public
source.
• Compute per-county aggregation: for each county polygon, compute the total
vehicle distance traveled and average speed. (Hint: use a spatial join to assign
county, then groupby.)
This will test your understanding of distributed spatial joins with large reference layers.
Shall I proceed to the Remote Sensing and Machine Learning chapter, where we classify
land cover from Sentinel-2 imagery using a Random Forest? Or do you wish to continue
refining the big data pipeline? I am, as ever, at your direction.
You are a remote sensing analyst tasked with producing a land cover map for a
20 km × 20 km region containing forests, agriculture, urban areas, and water bodies. You
have:
Your objective: train a Random Forest classifier on the spectral signatures of the training
polygons, then apply it to every pixel in the image to produce a classified GeoTIFF. Finally,
assess accuracy and generate statistics.
The alternative is deep learning (convolutional neural networks), but Random Forest
remains the state-of-the-art for per-pixel classification with moderate training data and is
far easier to implement and interpret.
3. Data Architecture
We need:
• Sentinel-2 image: We will simulate a multi-band GeoTIFF using synthetic data (real
Sentinel-2 would be a 12-band stack, but for education, we generate realistic
spectral signatures per land cover class and plant them in a raster).
All data must be in the same projected CRS (e.g., UTM) with pixel grid perfectly aligned.
4. Pipeline Design
1. Generate synthetic Sentinel-2 imagery with known classes, then write as
multi-band GeoTIFF.
3. Build feature matrix X (spectral bands per pixel) and target vector y (class labels).
5. Predict over the entire image by reading all bands, reshaping to (pixels, bands),
applying the model, and reshaping back to the image grid.
7. Assess accuracy using a confusion matrix and kappa coefficient (optional but
recommended).
5.1 Imports
python
import numpy as np
import rasterio
import pandas as pd
python
# Parameters
# Define spectral signatures: mean and std for each band [Blue, Green, Red, RedEdge, NIR,
SWIR]
signatures = {
'forest': ([500, 800, 400, 2000, 4000, 1500], [50, 80, 40, 200, 300, 150]),
'agriculture':([800, 1000, 1200, 2500, 3500, 1800], [100, 120, 150, 300, 350, 200]),
'urban': ([1200, 1100, 1300, 1500, 1800, 1000], [80, 70, 90, 100, 120, 80]),
'water': ([1500, 1200, 900, 700, 500, 300], [40, 30, 20, 15, 10, 5])
# Create ground truth label map (0: forest, 1: agriculture, 2: urban, 3: water)
# Center of image
bands = []
for b in range(6):
n_pixels = [Link](mask)
[Link](band)
img_array = [Link](bands)
profile = {
'driver': 'GTiff',
'height': height,
'width': width,
'count': 6,
'dtype': 'uint16',
'crs': crs,
'transform': transform_img,
'compress': 'lzw'
for i in range(6):
[Link](img_array[i], i+1)
Why this generation method: It ensures our classifier has meaningful spectral differences
to learn. The concentric pattern is easy to visualize and validate.
We create random points within known class regions, buffer them into small polygons, and
assign the correct class. This simulates field surveys.
python
[Link](42)
train_points = []
train_gdf.to_file('training_polygons.gpkg', driver='GPKG')
Why buffer: A single pixel may be noisy; by buffering 30 m (3 pixels), we average the
spectral signature over a small area, making training more robust.
Now the critical step: for each training polygon, we extract all pixel values from the image.
We’ll use [Link] to crop the image to each polygon, then flatten and collect.
For hundreds of polygons, looping with mask is acceptable. For thousands, we'd
use rasterstats for speed, but manual control teaches the inner workings.
python
samples = []
labels = []
geom = [Link]
[Link](pix)
[Link](row['class_name'])
X = [Link](samples, dtype=np.float32)
y = [Link](labels)
Why: We convert the extracted pixels into a feature matrix X (n_samples, n_bands) and
target vector y. This is the standard input for scikit-learn.
We split the data to evaluate the classifier before full image prediction.
python
[Link](X_train, y_train)
# Validation accuracy
y_pred = [Link](X_val)
print("Classification Report:")
print(classification_report(y_val, y_pred))
Parameter choices:
• max_depth=10: prevents overfitting to noise; the default None can overfit on easy
data.
We must now read the entire image, reshape to (pixels, bands), predict, and reshape back.
python
profile = [Link]
crs = [Link]
transform = [Link]
# Predict
class_ids = [Link](pixels)
Why reshape_as_image: rasterio stores bands as the first axis (BIP), whereas most ML
expects (pixels, bands). The helper aligns the axes.
We write the classification result as a single-band integer raster with a color interpretation.
python
# Create a dictionary mapping class name to integer (alphabetical order may not match our
initial numeric order)
unique_classes = sorted([Link](y))
class_map_int = [Link](class_to_int.get)(class_map)
else:
out_profile = [Link]()
[Link](class_map_int.astype('uint8'), 1)
colors = {
python
total_pixels = class_map_int.size
class_name = int_to_class[cls_id]
python
importances = clf.feature_importances_
band_names = ['Blue','Green','Red','RedEdge','NIR','SWIR']
print(f"{name}: {imp:.4f}")
We can plot the classified map using matplotlib for a quick look.
python
[Link](figsize=(8,6))
[Link](ticks=range(len(unique_classes)), label='Class')
[Link]()
1. Extract training data carefully – use [Link] or rasterstats to get pure pixel
vectors. Avoid including boundary mixed pixels if possible.
4. Predict in chunks for large images – loop over windows rather than loading the
entire raster into memory. Use [Link](window=window) for each tile.
5. Validate with independent polygons – never test on pixels used for training.
6. Use raster block processing for very large rasters; write results incrementally to a
new file.
8. Document the spectral signatures of your classes; they are crucial for
transferability.
• Create training polygons by selecting small areas manually in QGIS and exporting
to GeoPackage. (Or use a provided sample.)
• Add topographic features: compute slope and aspect from a DEM and append as
additional bands to improve classification in mountainous areas.
• Implement a tile-based prediction loop: read the image in 256×256 tiles, classify,
and write to output raster using [Link].
This challenge will test your integration skills and prepare you for operational remote
sensing workflows.
1. Read and resample all bands to a common 10-m grid, stacking them into a
10-band GeoTIFF.
2. Acquire a DEM (e.g., SRTM 1-arc-second) and compute slope and aspect,
resampled to the same grid.
3. Merge spectral bands and topographic bands into a single multi-band image.
4. Load training polygons (manual digitizing) and extract pixel values from this
merged image to build the training set.
6. Predict in tiles (256×256 pixels) to handle large images, writing directly to an output
GeoTIFF.
7. Assess accuracy using a separate validation polygon set and compute the
confusion matrix and kappa.
You will need: rasterio, geopandas, scikit-learn, numpy, scipy, matplotlib. We also
use os and glob for file handling.
python
import numpy as np
import rasterio
We will target a 10-m resolution GeoTIFF containing these bands (order matters):
B02, B03, B04, B05, B06, B07, B08, B8A, B11, B12.
We’ll use B02 as the spatial reference (10 m grid). We find all JP2 files, identify the
resolution, and for 20-m bands we reproject to the 10-m transform.
python
safe_path =
'S2A_MSIL2A_20230101T100031_N0500_R122_T33UUB_20230101T120000.SAFE'
# Inside: GRANULE/.../IMG_DATA/
# Define band mapping to file names (Sentinel-2 naming convention: e.g., ..._B02_10m.jp2)
band_map_10m = {
band_map_20m = {
files = [Link](pattern)
if not files:
return files[0]
# Load the 10m bands metadata from B02 to get reference CRS, transform, and
dimensions
ref_crs = [Link]
ref_transform = [Link]
ref_width = [Link]
ref_height = [Link]
ref_profile = [Link]
bands = []
band_names = []
f = find_file(band_map_10m[bname])
[Link]([Link](1))
band_names.append(bname)
# 2. Resample 20m bands to 10m grid
f = find_file(band_map_20m[bname])
reproject(
source=[Link](1),
destination=dest,
src_transform=[Link],
src_crs=[Link],
dst_transform=ref_transform,
dst_crs=ref_crs,
resampling=[Link]
[Link](dest)
band_names.append(bname)
img_stacked = [Link](bands)
Why bilinear resampling: For continuous spectral data, bilinear preserves radiometry
better than nearest neighbor. For the classification itself we may later read as integers, but
the resampling is a separate step.
We now write this stacked image as a multi-band GeoTIFF for later use.
python
out_meta = ref_profile.copy()
out_meta.update(count=len(bands), dtype=img_stacked.dtype, compress='lzw')
stacked_path = 'sentinel2_stacked.tif'
for i in range(len(bands)):
[Link](img_stacked[i], i+1)
We need a DEM covering the same extent. Use SRTM or ALOS. We'll assume you have a
GeoTIFF [Link]. We must reproject it to match the Sentinel-2 grid and then compute slope
and aspect.
python
dem_path = '[Link]'
ref_transform = src_template.transform
ref_crs = src_template.crs
ref_width = src_template.width
ref_height = src_template.height
reproject(
source=dem_src.read(1),
destination=dem_reproj,
src_transform=dem_src.transform,
src_crs=dem_src.crs,
dst_transform=ref_transform,
dst_crs=ref_crs,
resampling=[Link]
combined_path = 'combined_features.tif'
combined_profile = out_meta.copy()
combined_profile.update(count=combined_bands.shape[0], dtype='float32',
compress='lzw')
for i in range(combined_bands.shape[0]):
[Link](combined_bands[i].astype('float32'), i+1)
python
train_gdf = gpd.read_file('training_polygons.gpkg')
if train_gdf.crs != ref_crs:
train_gdf = train_gdf.to_crs(ref_crs)
samples = []
labels = []
geom = [Link]
[Link](pix)
[Link](row['class_name'])
X = [Link](samples, dtype=np.float32)
y = [Link](labels)
Why use NaN as nodata: We set nodata to NaN when we wrote the float32 raster, so we
can easily filter invalid pixels.
python
[Link](X_train, y_train)
print("Model trained.")
We will read the combined raster in 256×256 pixel windows, predict each window, and
write the classified tile to the output raster. This prevents loading the entire image into
memory.
python
out_cls_profile = combined_profile.copy()
class_names = sorted([Link](y))
tile_size = 256
w = Window(col_start, row_start,
# Predict
pred_str = [Link](pixels)
[Link](pred_img, 1, window=w)
Why tiling: For a 10 m resolution image covering 20×20 km, we have 2000×2000 pixels,
which easily fits in memory. But a full Sentinel-2 tile is 10980×10980 pixels at 10 m,
requiring ~1.2 GB just for 10 bands. Tiling keeps memory usage constant.
8. Step 6: Accuracy Assessment with Independent Validation Polygons
We need a separate set of polygons, ideally not overlapping the training ones. We'll load
them and extract pixel values exactly as we did for training, then predict and compute
metrics.
python
val_gdf = gpd.read_file('validation_polygons.gpkg')
if val_gdf.crs != ref_crs:
val_gdf = val_gdf.to_crs(ref_crs)
X_val_independent = []
y_val_true = []
geom = [Link]
valid = ~[Link](out_img).any(axis=0)
X_val_independent.append(pix)
y_val_true.append(row['class_name'])
y_val_pred = [Link](X_val_independent)
print("Confusion Matrix:")
print("\nClassification Report:")
print(classification_report(y_val_true, y_val_pred, target_names=class_names))
COG enables efficient streaming over HTTP. We can convert the output classification raster
using rasterio's COG driver or use the rio-cogeo plugin. Simple approach: when writing,
set driver='COG' (rasterio>=1.3.0 supports it directly) and specify overviews.
python
cog_profile = out_cls_profile.copy()
cog_profile.update(
driver='COG',
tiled=True,
blockxsize=256,
blockysize=256,
compress='lzw',
overviews='AUTO',
# Write tiles
data = [Link](window=w)
[Link](pred_int, 1, window=w)
Why COG: A COG stores the image in tiles with internal overviews (pyramid). A web map
server can fetch only the required tiles at the appropriate resolution, making it ideal for
interactive maps.
• Use sentinelhub Python library to fetch the image directly by bounding box and time
range, skipping manual download.
Your task now is to adapt this script to your own area of interest. Obtain a Sentinel-2 image
(from Copernicus Open Access Hub or sentinel-hub-py), create training and validation
polygons in QGIS (or use existing land cover maps), and run the pipeline. Then, critically
evaluate the confusion matrix. Which classes are most confused? Consider adding more
features (e.g., NDVI, NDWI, texture) to address those.
• NDVI (Normalized Difference Vegetation Index) = (NIR – Red) / (NIR + Red) — high for
healthy green vegetation, low for soil/water/urban.
By adding these indices as extra bands, the Random Forest can exploit nonlinear
combinations of raw bands more easily.
Texture features quantify spatial patterns within a neighborhood. A single pixel’s spectral
values can be noisy; texture describes the variation around it. The most common method is
the Gray-Level Co-occurrence Matrix (GLCM), from which we derive measures
like contrast, dissimilarity, homogeneity, energy, correlation, and entropy. For a 10 m
image, a 3×3 or 5×5 window captures local texture that distinguishes rough forest canopy
from smooth water or structured urban areas.
In our earlier classification, confusion often arises where spectral signatures overlap. For
instance, wet soil and asphalt may look similar in raw bands, but their texture differs
dramatically—asphalt has low local variance, while soil may show micro-patterns.
Similarly, NDVI can separate vegetation from non-vegetation better than raw NIR and Red
separately because the ratio compensates for shadows and illumination changes. Adding
these derived features gives the model a richer, more discriminative representation.
We already have a stacked image with bands in order: B02, B03, B04, B05, B06, B07, B08,
B8A, B11, B12. Let's recall their meanings: B04=Red, B08=NIR, B03=Green, B11=SWIR1,
B12=SWIR2.
We'll compute NDVI, NDWI, and NDBI as additional layers. We'll do this directly on
the img_stacked array.
python
red = img_stacked[2].astype('float32')
nir = img_stacked[6].astype('float32')
green = img_stacked[1].astype('float32')
swir1 = img_stacked[8].astype('float32')
eps = 1e-6
Why these three: NDVI for vegetation, NDWI for water, NDBI for urban—they target the
main classes we expect.
python
"""
"""
pad = window_size // 2
for r in range(rows):
for c in range(cols):
if win_max - win_min == 0:
entropy_val = 0.0
else:
entropy[r, c] = entropy_val
return entropy
# Compute texture on two representative bands
Performance note: The above loop is extremely slow for large images. In practice, you
would:
• Compute texture only on a subset of pixels (training samples) rather than the full
image, then use those values to train the classifier; at prediction time, you can either
pre-compute the full texture rasters or skip texture if real-time speed is needed.
For our educational purpose, we will pre-compute texture on the training pixels
only—extracting texture windows around each sampled pixel—which is far more
efficient and still enriches the training set.
Better approach: During training data extraction, for each pixel coordinate, extract a small
window from the original bands and compute texture features on the fly. We'll implement
this optimized method.
Instead of creating full-size texture rasters, we modify the training extraction loop. For each
sampled pixel, we extract a 5×5 neighborhood from the stacked image, compute GLCM
features, and append both the spectral value and the texture values to the feature vector.
We'll go back to the step where we loop over training polygons and
use [Link]. The mask function returns a window of the image; we can iterate
over every valid pixel within that window, extract its neighborhood from the original full
image, compute texture, and build the sample. Since we are already looping, this is
manageable.
python
# Re-define window size for texture
texture_window = 5
pad = texture_window // 2
samples = []
labels = []
geom = [Link]
# Mask to get pixels within polygon (but we need to know their row/col)
# Using mask may shift coordinates; better to rasterize polygon to get pixel indices
mask_poly = geometry_mask(
# Spectral: full_img[:, r, c]
# pad with edge reflection (numpy doesn't do that easily), skip edge pixels
continue
red_window = window_bands[2].astype('float32')
nir_window = window_bands[6].astype('float32')
def glcm_homogeneity(win):
if win_max - win_min == 0:
return 0.0
tex_red = glcm_homogeneity(red_window)
tex_nir = glcm_homogeneity(nir_window)
# But we don't have indices here because indices weren't precomputed per pixel;
# better to compute indices on the full image beforehand and stack them,
[Link](feature_vec)
[Link](row['class_name'])
Why this method: It avoids generating huge texture rasters and focuses computation only
on training/validation pixels. For prediction, we either pre-compute full texture rasters
(which can be done tile-by-tile) or we use only spectral+index features for prediction (a
trade-off). For this course, we'll assume we pre-compute texture rasters using a faster
method like [Link].convolve2d with GLCM alternatives (e.g., local variance as a
simple texture proxy). I'll show a faster alternative: local standard deviation as a simple
texture measure.
python
Why local std: High texture areas (urban, forest edges) exhibit high standard deviation;
smooth areas (water, bare soil) have low std. It is a simple, fast, and interpretable texture
measure.
We modify the creation of combined_features.tif to include NDVI, NDWI, NDBI, and local
std for Red, NIR, and SWIR1.
python
# Compute indices
features = [Link]([
img_stacked,
feature_profile = ref_profile.copy()
for i in range([Link][0]):
[Link](features[i], i+1)
We repeat the training extraction using this new feature set, then train the Random Forest
and compare accuracy. We'll also compare feature importance to see which new features
helped.
python
# ... same masking logic but using the new raster ...
1. Start with simple indices (NDVI, NDWI) before complex ones; they often bring the
most gain.
2. Texture window size should match the scale of spatial variation—5×5 works for
10 m resolution.
3. Normalize features if using algorithms like SVM; for Random Forest, it's optional.
4. Avoid data leakage: compute texture only on the training set if you won't have full
texture rasters at inference time. Here we pre-compute the full raster, which is fine
because we predict on the same image.
5. Use feature importance to prune useless features and reduce overfitting.
You now have a state-of-the-art land cover classification pipeline. The challenge is to apply
it to a new area, experimenting with feature combinations and evaluating what improves
accuracy most.
If you wish, I can now guide you through spatial optimization—for example, finding the
best location for a new facility using multi-criteria raster overlay—or geostatistical
interpolation (kriging) to create a continuous surface from point measurements. Which
path shall we explore next?
You are tasked with identifying the most suitable locations for a new wind farm within a
20 km × 20 km study area. The decision must balance:
Each criterion is represented as a raster layer. You must combine them into a
single suitability score per pixel, then identify contiguous patches of high suitability that
meet a minimum area (e.g., 5 hectares).
This is a classic weighted linear combination model: each input raster is normalized to a
common scale (e.g., 0–100), then multiplied by a weight reflecting its importance, and
summed. The result is a final suitability raster. Finally, vectorization extracts candidate
sites.
3. Data Architecture
We will generate synthetic rasters for five criteria, all sharing the same CRS and grid:
We’ll also create a binary constraint mask (e.g., water bodies and steep slopes >15° are
completely excluded).
Finally, we’ll produce a suitability raster, vectorize it into candidate polygons, and filter by
area.
4. Step-by-Step Implementation
python
import numpy as np
import rasterio
# Grid parameters
We model a prevailing wind from the southwest, gradually decreasing towards the
northeast, with some random noise and a ridge where wind accelerates. Wind values
between 3 and 12 m/s.
python
[Link](0)
We create a synthetic DEM with gentle undulations and a couple of steep ridges.
python
We randomly place a few road lines, rasterize them, and compute Euclidean distance.
python
road_geom = [
# Rasterize roads
Protected area polygons (randomly placed circles) are rasterized, and distance computed
similarly.
python
protected_polys = [
We synthesize land cover using noise and thresholds (simulating a classification map).
Classes: 0=water, 1=forest, 2=agriculture, 3=open, 4=urban.
python
[Link](1)
python
[Link]([Link](dtype), 1)
write_raster('[Link]', wind)
write_raster('[Link]', slope)
write_raster('dist_roads.tif', dist_roads)
write_raster('dist_protected.tif', dist_protected)
Each criterion has different units and ranges. We must transform them to a suitability
score where 100 means “most suitable” and 0 means “least suitable”, according to the
goal:
• Slope: threshold-based; slopes < 5° optimal (100), 5–10° moderate (50), >15° poor
(0).
We'll also create a constraint mask (pixels excluded entirely) for water, urban, and slopes
> 20°.
python
slope_suit = np.zeros_like(slope)
# Distance to roads: closer is better; we use a decay function: suitability = 100 * exp(-dist /
5000)
road_suit = 100 * [Link](-dist_roads / 5000.0)
lc_suit[lc == 2] = 80 # agriculture
lc_suit[lc == 1] = 20 # forest
lc_suit[lc == 0] = 0 # water
lc_suit[lc == 4] = 0 # urban
Why exponential decay for roads: Transport costs often follow a distance-decay;
exponential ensures sharp drop-off beyond a certain range.
• Wind: 0.35
• Slope: 0.25
python
weights = {
'wind': 0.35,
'slope': 0.25,
'roads': 0.20,
'protected': 0.10,
'landcover': 0.10
weights['slope'] * slope_suit +
weights['roads'] * road_suit +
weights['protected'] * prot_suit +
weights['landcover'] * lc_suit)
suitability[constraint_mask] = 0.0
write_raster('[Link]', suitability)
python
threshold = 70
results = (
min_area_ha = 5.0
sites_gdf.to_file('candidate_sites.gpkg', driver='GPKG')
Why threshold and area filter: A wind farm needs a contiguous area of adequate size.
Small isolated patches are infeasible.
We can quickly plot the suitability and the final sites using matplotlib.
python
fig, ax = [Link](1, 2, figsize=(12, 5))
ax[0].set_title('Suitability Score')
[Link]()
2. Define constraints early – excluded areas must be set to zero suitability and
removed before vectorization.
3. Validate weights – use sensitivity analysis (vary weights and observe changes) to
ensure robustness.
6. Check for edge effects – distance computations assume the raster extent is the
study area boundary; consider adding a buffer.
7. Area unit consistency – CRS must be projected (meters) for meaningful area
calculations.
Extend this pipeline to use real wind speed data (e.g., Global Wind Atlas GeoTIFF), a real
DEM (SRTM), OpenStreetMap roads, and a protected area layer (WDPA). Implement
the Analytic Hierarchy Process (AHP) to derive weights from pairwise comparison
matrices instead of fixed weights. Finally, use [Link] to attach the nearest road
distance and average wind speed to each candidate polygon and output a report.
This will test your ability to integrate multiple real-world datasets, perform advanced raster
math, and produce a polished deliverable.
You have 200 soil moisture sensors distributed across a 5 km × 5 km agricultural field. Each
sensor reports a moisture percentage. You must:
• Provide an uncertainty map (kriging variance) showing where predictions are less
reliable.
• Ensure the interpolation respects the spatial structure of the data (i.e., nearby
points are more similar).
Ordinary kriging is the optimal linear unbiased predictor under assumptions of stationarity
and a known variogram. It outperforms inverse distance weighting because it adapts to the
spatial correlation pattern.
2. Why Kriging?
Simple interpolation methods (IDW, spline) are deterministic and ignore the spatial
autocorrelation structure. Kriging models that structure via the variogram (or
semivariogram), which describes how data similarity decays with distance. It yields best
linear unbiased estimates and, crucially, kriging variance as a measure of confidence.
This uncertainty information is vital for decision-making—for instance, we may avoid
scheduling irrigation in areas of high uncertainty.
3. Theory in Brief
The core of kriging is solving a system of linear equations based on the variogram model.
The steps:
1. Compute an experimental variogram from the data: for all pairs of points
separated by distance ℎ, compute the average squared difference. This yields a set
of points (ℎ, 𝛾(ℎ)).
4. Solve for weights, compute the predicted value and kriging variance.
We will use the PyKrige library, which implements ordinary kriging efficiently in Python.
We'll also explore scikit-learn's GaussianProcessRegressor as an alternative that unifies
kriging with a broader probabilistic framework.
4. Data Generation
We create a synthetic true moisture field (a continuous function) and sample 200 points
from it with random noise to simulate sensor measurements.
python
import numpy as np
[Link](42)
n_sensors = 200
sensor_x = [Link](xmin, xmax, n_sensors)
5. Variogram Analysis
We first compute the experimental variogram and fit a model. PyKrige can do this internally,
but we'll perform it manually for understanding.
python
# bin by distance
mean_gamma = []
mean_dist = []
for i in range(n_lags):
if [Link](mask):
mean_gamma.append([Link](gamma[mask]))
mean_dist.append([Link](dists[mask]))
lags = 15
# Plot and fit a spherical model (using simple manual fit or [Link])
if h <= range_:
else:
return sill
# Plot variogram
[Link]('Semivariance')
[Link]()
[Link]()
Why manual variogram: To appreciate the model selection. In production, you'd let
PyKrige fit it automatically.
Now we perform kriging on a regular grid using the fitted variogram parameters.
python
grid_res = 5
variogram_model='spherical',
nlags=lags)
Key point: variogram_parameters are [sill, range, nugget] for the chosen model. The order
matches PyKrige's specification.
python
axes[0].set_title('True Moisture')
[Link](im1, ax=axes[0])
axes[1].set_title('Kriging Prediction')
[Link](im2, ax=axes[1])
axes[2].set_title('Kriging Variance')
[Link](im3, ax=axes[2])
[Link]()
python
# Cross validation
variogram_model='spherical',
rmse = [Link]([Link](cv_results['residuals']**2))
python
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True)
[Link](X_train, sensor_values)
# Predict on grid (this can be slow for large grids; use batch prediction)
1. Exploratory data analysis first – check for trends, anisotropy, and outliers.
Consider detrending before kriging if a strong trend exists.
4. Handle anisotropy if the spatial correlation varies with direction (e.g., along wind
direction). PyKrige supports anisotropic models.
6. Be aware of edge effects – predictions near the boundary have higher variance;
convey this with the uncertainty map.
7. Save the variogram model parameters – they are essential for reproducibility and
for applying the model to future data.
Obtain a real dataset of soil moisture samples (e.g., from a public soil database or an IoT
sensor network). Perform:
• Compute and fit an anisotropic variogram if the data show directional dependence.
This will solidify your geostatistics skills and give you a powerful tool for environmental
monitoring.
ADVANCED PROBLEM 6 — Building a Geospatial REST API with FastAPI and DuckDB
Spatial
You have a dataset of points of interest (POIs), road networks, and land parcels. You want
to expose spatial query endpoints:
Users should receive standard GeoJSON. The service must be fast, lightweight, and easily
deployable. We will use FastAPI (a modern async web framework) and DuckDB (an
in-process analytical database with spatial extension) to handle spatial queries without a
separate database server.
FastAPI is built for speed, automatic OpenAPI documentation, and async support. It's the
best choice for building Python APIs in 2025. DuckDB is a single-file database that can
read GeoParquet, CSV, and JSON directly, with spatial functions
(ST_Intersects, ST_Distance, ST_Buffer) available via the spatial extension. Together, they
allow us to embed a fully-featured spatial query engine inside a tiny web service.
3. System Architecture
• Data preparation: Convert POIs and roads to GeoParquet files (or load from existing
shapefiles).
• Database layer: DuckDB with spatial extension, registered in the FastAPI app
lifecycle.
• API layer: FastAPI endpoints that construct SQL queries, execute them, and return
GeoJSON using geopandas (or directly format JSON).
4. Step-by-Step Implementation
python
import duckdb
import json
con = None
FastAPI's startup event is perfect to load DuckDB, install spatial, and register data files.
python
@app.on_event("startup")
global con
con.install_extension('spatial')
con.load_extension('spatial')
# Load data from GeoParquet files (should exist in the same directory)
try:
# POIs
# Roads
except Exception as e:
def shutdown_event():
if con:
[Link]()
Why GeoParquet: It preserves CRS and geometry columns, loads instantly, and is
queryable with DuckDB's native parquet reader.
We'll write a function that takes a DuckDB result (list of tuples) and a geometry column
name, and returns a GeoJSON-like dict. We'll use geopandas for reliable conversion, but to
avoid heavy dependencies we can also build the GeoJSON manually
using [Link].
For simplicity, we'll use geopandas to read the query result and then .to_json().
python
"""Convert DuckDB result to GeoJSON dict. query_result must have a geometry column
as WKB."""
# DuckDB spatial functions return geometry as WKB, we can parse with shapely
features = []
pass
# Better: use geopandas read_postgis? DuckDB can output GeoDataFrame via arrow
DuckDB can return an Arrow table, which geopandas can read. Simpler: execute
with fetchdf() (which returns pandas DataFrame) and then convert. But DuckDB geometry
columns in pandas become bytes (WKB). We'll use [Link].from_wkb.
python
def query_geodataframe(sql: str) -> [Link]:
df = [Link](sql).fetchdf()
geom_col = None
geom_col = col
break
if geom_col is None:
df[geom_col] = df[geom_col].apply([Link])
return gdf
python
@[Link]("/pois")
):
# or project to a local UTM. Easiest: use DuckDB's ST_DWithin with spheroid distance?
# DuckDB spatial currently supports planar CRS; for geographic, we can cast to
EPSG:3857 (web mercator) and use meters (approximate).
# We'll transform the input point to 3857, compute buffer, and intersect.
sql = f"""
SELECT *, geom
FROM pois
WHERE ST_DWithin(
ST_Transform(geom, 'EPSG:3857'),
{radius_m}
"""
try:
gdf = query_geodataframe(sql)
return JSONResponse(content=[Link](gdf.to_json()))
except Exception as e:
Why transform to 3857: ST_DWithin in DuckDB works in the units of the CRS. For
EPSG:4326, radius would be in degrees. To use meters, we reproject to a projected CRS.
EPSG:3857 is good enough for local radius queries (<100 km). For ultimate accuracy, we'd
compute a suitable UTM zone dynamically.
We'll accept a GeoJSON polygon in the request body, use it to filter parcels.
python
from fastapi import Body
class GeoJSONFeature(BaseModel):
geometry: dict
properties: dict = {}
@[Link]("/spatial-join")
geom_dict = [Link]
polygon = shape(geom_dict)
wkt = [Link]
sql = f"""
SELECT *, geom
FROM parcels
"""
try:
gdf = query_geodataframe(sql)
return JSONResponse(content=[Link](gdf.to_json()))
except Exception as e:
@[Link]("/nearest-road")
):
sql = f"""
SELECT *, geom,
ST_Distance(
ST_Transform(geom, 'EPSG:3857'),
) AS distance_m
FROM roads
LIMIT 1
"""
try:
gdf = query_geodataframe(sql)
return JSONResponse(content=[Link](gdf.to_json()))
except Exception as e:
python
@[Link]("/health")
def health():
bash
Then open [Link] to see the auto-generated Swagger UI. You can test
all endpoints directly.
6. Going to Production
• /buffer – POST a GeoJSON point and radius, return the buffered polygon as
GeoJSON.
• /zonal-stats – POST a GeoJSON polygon and a raster path (or raster name), return
mean, min, max of raster values within the polygon (using DuckDB's raster
extension if available, or rasterio).
This will make your API a full-fledged geoprocessing engine accessible over HTTP.
1. Problem Statement
A logistics company wants to monitor its delivery fleet of 50 vehicles. Each vehicle sends
its position (latitude, longitude, speed) every 5 seconds. The system must:
• Ingest live positions via WebSocket (or a simple HTTP POST) without overwhelming
the database.
• Store positions efficiently with spatial indexing for fast historical queries.
• Broadcast the latest positions to all connected web clients (a live map).
• Serve historical trajectories via REST (e.g., last hour for a specific vehicle).
We’ll simulate vehicle movement using Python’s asyncio and a simple moving point
generator, but the ingestion layer works identically with real IoT devices.
2. Architecture
• DuckDB for storing positions (each row: vehicle_id, timestamp, lon, lat, speed).
We’ll use a persistent file to survive restarts.
• Background task that simulates 50 vehicles moving along random routes, sending
positions via HTTP POST to the /ingest endpoint.
• WebSocket manager that keeps track of connected clients and broadcasts new
positions.
• Frontend – a simple HTML page with [Link] that connects to the WebSocket,
updates markers in real time, and offers a slider to fetch historical tracks.
3. Step-by-Step Implementation
We’ll build a single [Link] plus a static/[Link]. I’ll explain every design choice.
text
realtime_tracker/
├── [Link]
├── static/
│ └── [Link]
python
import asyncio
import json
import time
import random
import math
import duckdb
import uvicorn
app = FastAPI()
# Database connection
db = [Link]("[Link]")
db.install_extension('spatial')
db.load_extension('spatial')
[Link]("""
vehicle_id VARCHAR,
ts TIMESTAMP,
lon DOUBLE,
lat DOUBLE,
speed DOUBLE,
geom GEOMETRY
""")
Why DuckDB: We can run SQL queries on the live table while ingesting, and the spatial
index enables fast radius searches. The geom column is automatically maintained via a
trigger? DuckDB’s spatial extension doesn't auto-create geometry from lon/lat, so we’ll
insert ST_MakePoint(lon, lat).
python
class Position(BaseModel):
vehicle_id: str
lon: float
lat: float
speed: float
python
class ConnectionManager:
def __init__(self):
await [Link]()
self.active_connections.add(websocket)
self.active_connections.discard(websocket)
payload = [Link](message)
try:
await connection.send_text(payload)
except Exception:
[Link](connection)
manager = ConnectionManager()
Vehicles (or the simulator) send positions here. We’ll store in DuckDB and broadcast to
WebSocket clients.
python
@[Link]("/ingest")
# Store in DB
ts = [Link]([Link])
[Link]("""
INSERT INTO positions VALUES (?, ?, ?, ?, ?, ST_GeomFromText(?, 'EPSG:4326'))
await [Link]({
"type": "new_position",
"vehicle_id": pos.vehicle_id,
"ts": [Link],
"lon": [Link],
"lat": [Link],
"speed": [Link]
})
python
@[Link]("/ws")
await [Link](websocket)
try:
df = [Link]("""
FROM positions
GROUP BY vehicle_id
""").fetchdf()
initial = []
[Link]({
"vehicle_id": row["vehicle_id"],
"ts": row["last_ts"].timestamp(),
"lon": row["lon"],
"lat": row["lat"],
"speed": row["speed"]
})
while True:
except WebSocketDisconnect:
[Link](websocket)
python
@[Link]("/trajectory/{vehicle_id}")
df = [Link]("""
FROM positions
ORDER BY ts
""", (vehicle_id, since)).fetchdf()
if [Link]:
geojson = {
"type": "Feature",
"geometry": {
"type": "LineString",
"coordinates": coords
},
"properties": {
"vehicle_id": vehicle_id,
"start": [Link]().isoformat(),
"end": [Link]().isoformat()
return JSONResponse(content=geojson)
A simple Leaflet map that connects to the WebSocket, updates markers, and fetches
trajectory on click.
html
<!DOCTYPE html>
<html>
<head>
<script src="[Link]
</head>
<body>
<script>
[Link]('[Link]
[Link] = function(event) {
[Link](function(v) {
addOrUpdateMarker(v);
});
addOrUpdateMarker(data);
};
function addOrUpdateMarker(pos) {
var id = pos.vehicle_id;
if (markers[id]) {
markers[id].setLatLng([[Link], [Link]]);
markers[id].getPopup().setContent('Vehicle ' + id + '<br>Speed: ' +
[Link](1) + ' m/s');
} else {
.addTo(map)
[Link]('click', function() {
loadTrajectory(id);
});
markers[id] = marker;
function loadTrajectory(id) {
fetch('/trajectory/' + id + '?minutes=30')
.then(data => {
if ([Link]) {
});
</script>
</body>
</html>
3.9 Vehicle Simulator (Run Separately)
To test the system, we create a script [Link] that generates moving vehicles.
python
import asyncio
import aiohttp
import random
import math
import time
BASE_URL = "[Link]
# Each vehicle has a current position and a heading that changes occasionally
vehicles = {}
vehicles[vid] = {
v = vehicles[vid]
# Move
# Send to API
data = {
"vehicle_id": vid,
"timestamp": [Link](),
"lon": v["lon"],
"lat": v["lat"],
"speed": v["speed"]
pass
while True:
await [Link](*tasks)
if __name__ == "__main__":
[Link](simulate())
1. WebSocket for push, REST for pull – Use WebSocket for low-latency broadcasts,
REST for historical queries.
3. Spatial indexing – Always create spatial indices for radius and nearest-neighbor
queries.
5. Asynchronous all the way – Use asyncio and async database drivers (DuckDB’s
Python API is synchronous but fast; for massive concurrency, consider a connection
pool).
• Add a /vehicles endpoint that returns the latest position of all vehicles (faster than
full history).
• Containerize the whole system with Docker Compose (add Redis for pub/sub if
needed for scale).
Now you have a complete real-time geospatial system. From here, the step to production
involves adding authentication, persistent storage, and deploying to a cloud VM.
You began with no programming knowledge—not even what an IDE was. Now you can:
• Write clean, professional Python with proper naming, docstrings, error handling,
and testing.
• Apply machine learning to satellite imagery for land cover classification with
feature engineering.
• Build scalable pipelines with Dask and Dask-GeoPandas for datasets too large for
memory.
• Deploy APIs and real-time systems with FastAPI, WebSockets, and DuckDB
spatial.
2. Scale to Production
• Personal Project: Create a web app that solves a problem in your community—air
quality monitoring, deforestation alerts, or flood early warning.
4. Stay Connected
• Join communities: GIS Stack Exchange, Spatial Community Slack, Pangeo for big
geoscience.
"You came asking how to write Python for geospatial. You now know not just the syntax, but
the why behind every line. Remember these principles:
• Start small, then scale. A script that works on a sample can become a production
system.
• Share your knowledge. Teach the next beginner; it will deepen your own
understanding.
You are no longer my student. You are a colleague. I look forward to seeing what you build—
perhaps one day I will use a library you created, or read a paper you authored.
The terminal is now yours. Write code that maps a better world."