Python Geospatial Analysis: 5 Practical Techniques to Master Location Data Today

#programming #devto #python #softwareengineering

As a best-selling author, I invite you to explore my books on Amazon. Don't forget to follow me on Medium and show your support. Thank you! Your support means the world!

I work with maps and location data almost every day. It started when I needed to figure out the best delivery routes for a small business. That led me down a path of discovering how Python can make sense of anything with a latitude and longitude. I want to share the main techniques I use. This isn't about complex theory. It's about practical code you can run today to solve real problems. Let's walk through them one by one.

First, you need to speak the language of maps. Every map uses a coordinate system. Think of it as the rules for a game of Battleship. "A-10" only means something if both players agree on where "A" starts and what "10" represents. The earth is a 3D sphere, but our screens and papers are flat. A Coordinate Reference System (CRS) is the set of rules that flattens the sphere.

The most common system is WGS84, which uses latitude and longitude. Your phone uses this. But for measuring distance or area on a local map, a projected system like UTM is better. It converts curves into a flat grid. The pyproj library is your translator between these systems.

import pyproj
from pyproj import CRS, Transformer

# Define two different "languages" for location.
lat_lon_system = CRS.from_epsg(4326)  # WGS84, what your phone uses.
utm_system = CRS.from_epsg(32610)     # UTM zone 10N, a flat grid for California.

# Create a translator.
translator = Transformer.from_crs(lat_lon_system, utm_system, always_xy=True)

# A point in San Francisco in latitude/longitude.
lat, lon = 37.7749, -122.4194

# Translate it to the flat UTM grid.
easting, northing = translator.transform(lon, lat)
print(f"SF in UTM: {easting:.0f}, {northing:.0f}")

Why does this matter? If you try to calculate the area of a field using latitude and longitude degrees, you'll get a nonsense answer. You must project it first. pyproj also calculates real-world distances on the curved earth, which is essential for logistics.

from pyproj import Geod

# Set up a calculator for the Earth's shape.
geod = Geod(ellps='WGS84')

# Coordinates for San Francisco and Los Angeles.
sf = (-122.4194, 37.7749)
la = (-118.2437, 34.0522)

# Get the true distance.
azimuth1, azimuth2, distance_meters = geod.inv(sf[0], sf[1], la[0], la[1])
distance_km = distance_meters / 1000
print(f"Distance from SF to LA: {distance_km:.1f} km")

Once locations are in a consistent system, you can ask spatial questions. This is where geometry comes in. Is this delivery point inside the service zone? How far is the fire hydrant from the building? If I expand this road by 5 meters, what does it touch? The Shapely library answers these questions.

Shapely deals with points, lines, and polygons. It's like having digital graph paper and a set of rulers.

from shapely.geometry import Point, Polygon, LineString

# Create shapes.
warehouse = Point(10, 20)
delivery_zone = Polygon([(0, 0), (50, 0), (50, 50), (0, 50)])
river = LineString([(60, 10), (40, 30), (20, 5)])

# Ask questions.
is_inside = warehouse.within(delivery_zone)
distance_to_river = warehouse.distance(river)

print(f"Is the warehouse in the zone? {is_inside}")
print(f"Distance to the river: {distance_to_river:.1f} units")

A common task is creating a buffer. Imagine drawing a circle of 1 km around a store. That's a buffer. You can then check which customers live within that circle.

# Create a 1.5 unit buffer around the warehouse.
service_area = warehouse.buffer(1.5)
print(f"Service area size: {service_area.area:.2f} square units")

# Check if a customer point is in the service area.
customer = Point(11, 20.5)
can_we_deliver = customer.within(service_area)
print(f"Can we deliver to the customer? {can_we_deliver}")

Real-world data is messy. You have city boundaries, road lines, and points of interest. You need a way to store them, mix them with non-spatial data (like population), and perform bulk operations. This is where GeoPandas shines. It's Pandas for maps.

A GeoDataFrame is a table where one column is a Shapely geometry. You can filter, merge, and group just like a normal DataFrame, but with spatial superpowers.

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

# Create some city data with locations.
data = {
    'city': ['SF', 'LA', 'NYC'],
    'population': [883000, 3979000, 8336000],
    'geometry': [
        Point(-122.4194, 37.7749),
        Point(-118.2437, 34.0522),
        Point(-74.0060, 40.7128)
    ]
}

# Make a GeoDataFrame.
gdf = gpd.GeoDataFrame(data, crs='EPSG:4326')
print(gdf.head())

One of the most powerful features is the spatial join. Let's say you have a file of coffee shops (points) and a file of neighborhood boundaries (polygons). A spatial join can tell you which neighborhood each coffee shop is in, all in one step.

# Imagine 'neighborhoods_gdf' is loaded from a file.
# This line joins each coffee shop point to the polygon that contains it.
# shops_with_hood_info = gpd.sjoin(coffee_shops_gdf, neighborhoods_gdf, how='left', predicate='within')

Not all data is made of shapes. Some, like satellite images or elevation models, are grids of pixels or cells. This is raster data. Each cell has a value, like temperature or height. The Rasterio library is the standard tool for this.

While a full raster analysis is deep, the basic idea is simple: read a grid, inspect it, and process it.

import rasterio
import numpy as np

# Typical workflow to read a GeoTIFF (e.g., a satellite image).
# with rasterio.open('elevation.tif') as src:
#     elevation_grid = src.read(1)  # Read the first band.
#     profile = src.profile  # Save the map's location info.

# Let's create a simple synthetic grid as an example.
height, width = 5, 5
# Make a small hill.
simple_hill = np.array([
    [100, 101, 102, 101, 100],
    [101, 103, 105, 103, 101],
    [102, 105, 110, 105, 102],
    [101, 103, 105, 103, 101],
    [100, 101, 102, 101, 100]
])

print("Our tiny digital elevation model:")
print(simple_hill)
print(f"Highest point: {simple_hill.max()} meters")

You can calculate slope, find valleys, or combine rasters. For example, you could overlay a rainfall raster on a soil-type raster to estimate erosion.

Maps are for finding your way. Network analysis is about modeling connections—roads, pipes, or social links. The core idea is a graph: intersections are nodes, and roads are edges with a cost (like travel time).

Python's NetworkX library is perfect for this. You can find the shortest path, the best location for a new store, or all areas reachable within a 10-minute drive.

import networkx as nx

# Create a simple road network.
G = nx.Graph()

# Add intersections (nodes).
G.add_node('A', pos=(0, 0))
G.add_node('B', pos=(2, 0))
G.add_node('C', pos=(2, 2))
G.add_node('D', pos=(0, 2))

# Add road segments (edges) with travel time.
G.add_edge('A', 'B', weight=4)  # 4 minutes
G.add_edge('B', 'C', weight=3)
G.add_edge('C', 'D', weight=3)
G.add_edge('D', 'A', weight=5)
G.add_edge('A', 'C', weight=7)  # A diagonal road.

# Find the fastest path from A to C.
fastest_path = nx.shortest_path(G, 'A', 'C', weight='weight')
travel_time = nx.shortest_path_length(G, 'A', 'C', weight='weight')
print(f"Fastest path from A to C: {fastest_path}")
print(f"Travel time: {travel_time} minutes")

A service area analysis is very useful. For an ambulance station, what addresses can be reached in 8 minutes?

# Calculate travel time from node 'A' to all others.
times = nx.single_source_dijkstra_path_length(G, 'A', weight='weight')
print("Travel times from station A:")
for node, time in times.items():
    print(f"  To {node}: {time} min")

Static maps have their place, but interactive maps let you explore. You can zoom, click, and toggle layers. Folium makes this easy by generating Leaflet.js maps with Python code.

You can start with a base map of the world and add your own data on top.

import folium

# Center the map on San Francisco.
sf_map = folium.Map(location=[37.7749, -122.4194], zoom_start=13)

# Add a marker for City Hall.
city_hall = [37.7793, -122.4192]
folium.Marker(
    location=city_hall,
    popup="San Francisco City Hall",
    tooltip="Click for info",
    icon=folium.Icon(color="green", icon="flag")
).add_to(sf_map)

# Draw a circle around it (e.g., a 1km radius).
folium.Circle(
    radius=1000,
    location=city_hall,
    color="blue",
    fill=True,
    fill_opacity=0.2
).add_to(sf_map)

# To see the map, save it to an HTML file.
# sf_map.save("my_first_map.html")
print("Map created! Open 'my_first_map.html' in your browser.")

Folium can color regions based on data (choropleth maps), plot lines from GPX files, or cluster thousands of points. It turns your analysis into something anyone can use and understand.

Finally, the true power comes from mixing these techniques. You might use GeoPandas to filter census tracts, Shapely to calculate their centroids, and NetworkX to find the shortest path between them. Or use Rasterio to load a satellite image, classify it with Scikit-learn, and publish the results with Folium.

Here's a small, combined example. Let's find public parks in an area and calculate a walkable area around them.

# This is a conceptual outline of the steps.
# 1. Use GeoPandas to load park polygons and street lines.
# parks_gdf = gpd.read_file('parks.geojson')
# streets_gdf = gpd.read_file('streets.geojson')

# 2. Find the main entrance (a point on the park boundary).
# park_entrances = parks_gdf.geometry.representative_point()

# 3. Build a walking network from the streets.
# (Convert street lines to a NetworkX graph)

# 4. For each entrance, calculate the area reachable by a 10-minute walk.
# service_areas = []
# for entrance in park_entrances:
#     area = calculate_walkable_area(graph, entrance, time_limit=600) # 600 seconds
#     service_areas.append(area)

# 5. Plot the parks and their walkable areas on a Folium map.

This workflow solves a real urban planning question: "Where should we focus park improvements to serve the most residents?"

The goal isn't to memorize every function. It's to know which tool to reach for. Need to change map projections? Use pyproj. Asking "inside or outside"? Use shapely. Working with a table of places? Use geopandas. Modeling roads or pipes? Use networkx. Making a map to share? Use folium.

Start with a small project. Plot your favorite coffee shops on a map. Calculate how far you ran last week. Find the optimal route for your weekend errands. Each problem you solve will teach you which of these eight techniques you need, and how they fit together. The code is here. The maps are waiting.

📘 Checkout my latest ebook for free on my channel!

Be sure to like, share, comment, and subscribe to the channel!

101 Books

101 Books is an AI-driven publishing company co-founded by author Aarav Joshi. By leveraging advanced AI technology, we keep our publishing costs incredibly low—some books are priced as low as $4—making quality knowledge accessible to everyone.

Check out our book Golang Clean Code available on Amazon.

Stay tuned for updates and exciting news. When shopping for books, search for Aarav Joshi to find more of our titles. Use the provided link to enjoy special discounts!