Shapely

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 1 of 8. You might be used to treating map geometries as basic lists of coordinates, used primarily to draw shapes on a screen. But when you need to calculate if a GPS track crosses a specific border, simple drawing logic falls apart. To do real mathematical spatial analysis, you need The Spatial Data Model. First, we must clear up a common source of confusion regarding this model. Shapely operates strictly on a flat, two-dimensional Cartesian plane. You can provide X, Y, and Z coordinates when creating your geometries. Shapely will parse them and even store that Z value in memory. But it completely ignores the Z-coordinate during any spatial analysis. If you ask for the distance between two points, Shapely calculates the flat 2D distance. It does not calculate 3D volumes or 3D distances. It is a purely planar library. Here is the key insight. Shapely is built on mathematical point-set theory. It does not view a shape as just the handful of vertices you typed into your code. Instead, it treats every geometry as a rigorous set containing an infinite number of points. Under this model, every geometry divides the entire 2D plane into three distinct parts. The interior is the actual substance of the shape. The boundary is the edge that frames the shape. The exterior is the rest of the infinite plane that is not part of the shape. Consider modeling a city park to see how these three concepts apply. You start with a water fountain, represented as a Point geometry. A Point is zero-dimensional. It represents a single, exact location on the plane. Because it has no length or area, its interior consists solely of that one specific coordinate. A Point has no boundary at all. Every other coordinate on the entire Cartesian plane makes up its exterior. Next, you add a walking path leading to the fountain. This is modeled as a LineString. A LineString is a one-dimensional geometry defined by an ordered sequence of coordinate vertices. Even if you only pass it a start and an end coordinate, Shapely understands the LineString as the infinite set of points forming the continuous path between them. The interior is the entire continuous length of the path. The boundary consists of exactly two points: the first starting coordinate and the last ending coordinate. Finally, you model the footprint of the park itself using a Polygon. A Polygon is a two-dimensional geometry with a measurable area. You define it using a closed ring of coordinates, meaning the first and last points are identical. This outer ring forms the boundary of the Polygon. The interior is the infinite set of points enclosed within that boundary. Polygons can also contain holes, such as a restricted conservation zone inside the park. The edge of that hole simply acts as an inner boundary, cleanly separating the valid interior of the park from the empty exterior space inside the hole. Shapely does not see your geometry as a hollow wireframe of coordinates; it sees a continuous region of space mathematically divided into interior, boundary, and exterior. If you find these episodes helpful and want to support the show, search for DevStoriesEU on Patreon. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 2 of 8. You are processing a batch of spatial coordinates. You write a loop to turn each coordinate pair into a shape, one at a time. It works, but as your dataset grows to millions of points, your application grinds to an absolute halt. The problem is Python overhead, and the fix requires bypassing standard object creation entirely. This episode is about Geometry Creation and Serialization. Creating shapes is the logical starting point. If you are dealing with a single location, instantiating a geometry is straightforward. You import the singular Point class and pass it an X and Y value. It is easy to read, but it is slow. Here is the key insight. Shapely provides vectorized universal functions that operate directly on NumPy arrays. Instead of using the singular Point class in a Python loop, you use the plural shapely dot points function. You pass it a two-dimensional array of coordinates, and it pushes the entire iteration down to the highly optimized C level. The result is an array of geometry objects returned in a fraction of the time. This applies across the board, with plural functions available for line strings, polygons, and linear rings. Say you are building a routing tool and need to load a dataset of a thousand delivery drop-off points. The data arrives as an array of GeoJSON payloads. This brings us to serialization. You need to parse those strings into actual geometry objects in memory. You might be tempted to use standard JSON parsing and extract the coordinates manually to build your shapes. Do not do that. Shapely has built-in input and output functions for the three major spatial formats: Well-Known Text, Well-Known Binary, and GeoJSON. Like geometry creation, these parsing functions are vectorized. You take your entire array of a thousand GeoJSON strings and pass it directly into the shapely dot from geojson function. It parses the whole batch at once and returns a high-performance array of geometries. If your data was in text format, you would use shapely dot from wkt. If it was in binary format, you would use shapely dot from wkb. That covers inputs. What about outputs? Once you have finished processing your drop-off points, you need to save them or send them to a database. You reverse the process using the export functions. If you are debugging and need to read the coordinates with your own eyes, you use shapely dot to wkt. Well-Known Text gives you a human-readable string, spelling out the geometry type followed by the coordinates. However, text parsing is expensive and takes up space. If you are sending this data to a spatial database, or storing it for later, you should use shapely dot to wkb. Well-Known Binary is the raw byte representation of the geometry. It is significantly smaller, it avoids precision loss from floating-point string conversion, and it is much faster for machines to read and write. You pass your geometry array into shapely dot to wkb, and you get an array of byte strings ready for storage. Whenever you are moving spatial data in or out of your application, remember that Python loops are the enemy. The single most important habit for performance is feeding full arrays into plural Shapely functions, letting the underlying C library do the heavy lifting in a single pass. Thanks for listening, happy coding everyone!

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 3 of 8. Regular distance tells you the closest two lines get to each other. But if you want to know the absolute furthest they drift apart, regular distance is completely useless. To answer that, you need to understand Measurement and Properties. Every geometry object in Shapely comes with built-in attributes that describe its physical dimensions. The most basic ones are area, length, and bounds. The area property gives you the two-dimensional space covered by a shape. For a polygon, this is a positive number. For points and lines, the area is exactly zero. The length property returns the one-dimensional span of a geometry. For a line, it is the distance from start to finish along the path. The bounds property provides the exact coordinate limits of your shape. It returns a four-number tuple representing the minimum X, minimum Y, maximum X, and maximum Y. This is your bounding box. Consider a scenario where you have two GPS tracks from runners, both represented as LineStrings. You read the length property to see the total distance each person ran. You read the bounds property to define the exact rectangular viewport needed to display both of their routes on a screen. The math gets heavier when you need to measure the space between those two tracks. The standard distance measurement returns the absolute shortest gap between two geometries. It scans track A, scans track B, and finds the single closest pair of points. If the two runners crossed paths at an intersection halfway through their run, the standard distance between their tracks is zero. Here is the key insight. People frequently confuse regular distance with Hausdorff distance. Regular distance finds the nearest points. Hausdorff distance finds the furthest distance between nearest points. Think of Hausdorff distance as measuring maximum drift. It represents the longest gap you would be forced to cross if someone placed you on the worst possible spot on one track and told you to walk to the nearest point on the other track. If our two runners started at the same place, split up by two miles in the middle of their workout, and finished at the same place, the regular distance is zero. The Hausdorff distance is two miles. It defines the worst-case separation across the entire length of the two geometries. There is one more measurement function called minimum clearance. This one does not compare two different objects. It measures the structural stability of a single geometry. Minimum clearance calculates the smallest distance any node can move before the geometry collapses into an invalid state, like a line overlapping itself. If your geometry has a tiny minimum clearance, the coordinates are squeezed tightly together. A slight rounding error during a database export could break the shape entirely. If you want to simplify the runner data to save file size, checking the minimum clearance tells you exactly how far you can shift the points before the geometry ruins itself. Knowing how close two shapes get handles the basic collisions. Knowing the absolute maximum limit of their separation dictates how well you truly understand the spatial relationship between them. Thanks for spending a few minutes with me. Until next time, take it easy.

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 4 of 8. You are routing a new subway line near a protected wetland. You need to know if the track cuts straight through the water or just skirts the edge. How does the computer actually figure that out mathematically? It uses a clever three-by-three matrix to answer the question. This brings us to spatial predicates and the Dimensionally Extended nine-Intersection Model, or DE-9IM. Spatial predicates are binary functions. They compare two geometries and return exactly one thing: True or False. They do not create new shapes. They do not calculate overlapping areas. They simply answer yes or no questions about topological relationships. Under the hood, Shapely evaluates these relationships using the DE-9IM. Every geometry has three parts: an interior, a boundary, and an exterior. The DE-9IM is a matrix that tests the intersections between these three parts of the first geometry against the three parts of the second geometry. Three parts times three parts gives you nine possible intersections. Depending on which of those intersections result in an empty set, a point, a line, or an area, Shapely determines the exact relationship. Let us apply this to the subway line and the wetland. The broadest predicate is intersects. If the track shares any point, line, or area with the wetland, intersects returns True. It is the baseline catch-all. If intersects is False, the two geometries are completely separate in space. But you often need more precision. You want to know if the subway line actively cuts through the wetland. For this, you use crosses. A line crosses a polygon if its interior intersects the polygon's interior, but the line also extends outside the polygon. The DE-9IM matrix checks that the intersection of their interiors is a line, and that the line's interior also intersects the polygon's exterior. If both are true, crosses returns True. What if the track runs exactly along the border of the wetland without entering the water? That is where touches comes in. Two geometries touch if they share at least one point, but their interiors do not intersect at all. Only their boundaries interact. If the track veers even one millimeter into the wetland's interior, touches becomes False. Now consider a subway station built completely inside the wetland. Here, you check containment. People often mix up contains and within, but they are just inverse operations of each other. If the wetland contains the station, then the station is within the wetland. Geometrically, it means the station's interior and boundary are entirely enclosed by the wetland's interior. No part of the station intersects the wetland's exterior. You can write your code checking if the wetland contains the station, or checking if the station is within the wetland. The result is exactly the same. Here is the key insight. You do not need to memorize the nine-intersection matrix to route your subway, but knowing it exists explains why these binary checks are so reliable. Every named spatial predicate is just a specific pattern of true and false flags in that underlying mathematical matrix. Thanks for hanging out. Hope you picked up something new.

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 5 of 8. Venn diagrams are not just for logical sets or database queries. You can use the exact same mathematical principles to physically carve up real-world geography. To do that, you need Set-Theoretic Operations. In Shapely, every geometry is essentially an infinite set of points on a two-dimensional plane. Because they are sets, you can manipulate them using classic set theory. First, we must clear up a common source of confusion. Spatial predicates, like checking if two shapes touch or cross, evaluate a condition and return a simple True or False. Set-theoretic operations do something completely different. They return entirely new geometry objects. They take two input shapes, analyze their overlapping points, and construct a brand new third shape as the output. Consider a concrete scenario. You have two polygon datasets. One represents an old flood zone map. The other is a newly updated flood risk model. You need to calculate the exact area where these two models agree. For this, you use the intersection method. You call intersection on the first geometry and pass the second geometry as the argument. Shapely looks at the spatial overlap and returns a new polygon containing only the points that exist in both the old map and the new map. If the two shapes do not overlap at all, it returns an empty geometry. Now, what if you need to find the strictly new risk areas? These are the neighborhoods that were safe under the old map but are now in the updated flood zone. Here, you use the difference method. You take your new flood model geometry and call difference, passing in the old map geometry. This returns a shape containing all the points in the new model, explicitly minus any points that were already in the old one. Order is critical here. If you reverse it and call difference on the old map, passing the new model, you get an entirely different result. You get the areas that are no longer considered at risk. It is just spatial subtraction. Sometimes you do not want to compare shapes, you just want to combine them. If emergency services need a master map of anywhere that has ever been flagged as a flood risk, you use the union method. You call union on one shape, pass the other, and Shapely merges them together. If the shapes overlap, the internal boundaries dissolve. The output is a single continuous geometry representing all points from both original shapes. Finally, there is the symmetric difference method. Think of this as the spatial equivalent of an exclusive OR operation. It returns a new geometry containing points that are in either the old map or the new map, but absolutely not in both. In our flood scenario, this single operation gives you exactly the areas of disagreement between the two models. It outputs the newly added risk zones alongside the newly removed risk zones, completely hollowing out the areas where the models agree. Here is the key insight. When you stop thinking about spatial boundaries as complex geographic concepts and start treating them as basic mathematical sets, calculating complex overlaps and exclusions becomes a predictable, completely standard process. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 6 of 8. A simple mathematical point has exactly zero area. It is just a location. But what if you need to represent a blast radius around that location? You can instantly transform that zero-dimensional point into a massive, precise polygon. This is the power of constructive operations, specifically buffers and hulls. Constructive operations take an existing geometry and build a new, different geometry out of it based on its spatial properties. You can think of them as safety wrappers around your spatial data. The most widely used constructive operation is the buffer. You pass a distance value to the buffer function, and Shapely returns a polygon representing all points within that exact distance from your original geometry. If you have a line representing the path of a chemical spill along a highway, and you need to establish a 50-meter safety zone, you call buffer on that line geometry with a distance of 50. The output is a new polygon that traces the entire route, expanding outward exactly 50 meters in every direction. Here is the key insight. Buffers do not just grow. They can also shrink. If you pass a positive distance, the geometry inflates. If you apply a negative distance to a polygon, the geometry deflates. A negative buffer pulls the boundaries inward. This is incredibly useful for finding the safe inner core of a zone, isolating the area that is a strict distance away from the dangerous edges. Sometimes you do not need an exact buffer around complex shapes. You just need the absolute limits of your data. Calling the envelope property gives you exactly that. It returns the bounding box, which is the smallest, perfectly upright rectangle that completely encloses your geometry. It is mathematically cheap to calculate and perfect for fast spatial indexing. But what happens when your data is not a solid shape, but a scattered collection of points? If you have dozens of contaminated soil samples and you need to draw a single continuous boundary around the affected area, you need a hull. Shapely gives you two main ways to do this. The convex hull calculates the smallest convex polygon that contains all your geometries. You can picture this as stretching a rubber band around a board full of pegs. The rubber band snaps tight across the outermost pegs. Because it is convex, the boundary never dents inward. There are no interior angles greater than 180 degrees. While the convex hull is reliable, it often includes a large amount of empty space if your points form a crescent or an irregular cluster. If you need a tighter fit, you use a concave hull. Instead of a rigid rubber band, a concave hull acts like shrink-wrap. It allows the boundary to cave inward to follow the actual footprint of your data. You control how tightly it wraps by adjusting a parameter that limits the maximum allowed length of the boundary edges. Tighter limits force the boundary to fold into the gaps between points, giving you a highly accurate, realistic map of that contaminated soil. The difference between these tools comes down to how closely they hug your data. The envelope gives you a quick loose box, the convex hull stretches tightly across the extreme outer points, the concave hull shrink-wraps the exact shape, and the buffer strictly pads the edges by a defined distance. If you want to help keep these episodes coming, you can support the show by searching for DevStoriesEU on Patreon. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 7 of 8. Try rendering a one-million-point coastline on a web map, and you will quickly crash the user's browser. You need a way to reduce that data footprint to one percent of its original size while making it look completely identical. This requires Advanced Constructive Operations. The core tool for this data cleanup is the simplify function. Listeners sometimes think simplifying a shape just means skipping every second or third point in an array. Doing that destroys the geometry. It creates self-intersecting lines, collapses small islands, and ruins the topological structure. Shapely's simplify function evaluates the actual geometry. You pass it your complex coastline and a tolerance value. The function removes vertices that do not deviate from the core shape by more than that tolerance. Crucially, Shapely gives you a flag to preserve topology. When enabled, the algorithm actively checks if removing a point will cause a polygon to cross over itself or split into invalid pieces. It guarantees your highly compressed coastline remains mathematically valid. Once you simplify a heavy base layer, you often need other layers to align with it perfectly. Say you have land parcels or property lines that end exactly at the water. Because you just simplified the coastline, its edges have shifted slightly. This creates microscopic overlaps or empty spaces, known as slivers. You fix this using the snap function. You provide your property line geometry, your coastline geometry, and a strict tolerance distance. Snap looks at the vertices of the first shape. If a vertex sits within the tolerance distance of a vertex or segment on the second shape, it gets pulled directly onto it. Points outside the tolerance stay exactly where they are. This creates a watertight, flush fit between adjacent geometries without distorting the rest of the shape. That covers cleaning up existing shapes. The other side of constructive operations is dividing empty space using points. If you take the vertices of a geometry, you can pass them into the delaunay triangles function. This connects those points to form a continuous mesh of non-overlapping triangles. The logic specifically draws these triangles to maximize the smallest internal angles. This prevents long, needle-like triangles, creating a balanced mesh highly useful for rendering terrain. If you only want the outline of these connections rather than filled shapes, you can set a flag to return just the line edges. Closely related is the voronoi polygons function. While Delaunay connects points, Voronoi draws boundaries between them. You provide a geometry, and Shapely uses its vertices to divide the surrounding area into completely distinct polygonal regions. Any location inside a specific Voronoi polygon is closer to its central generating point than to any other point in the entire set. If you are analyzing a scatter of delivery hubs, Voronoi polygons instantly map the exact coverage zone for each hub based on pure proximity. Just like Delaunay, you can instruct the function to return only the boundary edges if you do not need solid polygons. Here is the key insight. Real-world spatial data is dense, disconnected, and rarely aligns perfectly. Constructive operations like simplify, snap, and Voronoi give you the tools to enforce mathematical order, making the data light enough to render and precise enough to analyze. Thanks for spending a few minutes with me. Until next time, take it easy.

Hi, this is Alex from DEV STORIES DOT EU. Shapely, episode 8 of 8. Comparing a single point against a million polygons takes massive compute time, yet some applications do this in milliseconds. They do not calculate every distance. They filter out ninety-nine percent of the data before doing a single math equation. This is high-performance spatial indexing with STRtree. Say you want to find the nearest coffee shop out of ten thousand locations to a user's current GPS coordinate. If you calculate the exact distance to every single shop, your system will crawl to a halt. Instead, spatial indexing uses an R-tree. Think of an R-tree as a hierarchical folder system for two-dimensional space. It draws a simple rectangle, called a bounding box or envelope, around clusters of nearby geometries. Then it draws bigger boxes around those boxes. Checking if two simple rectangles overlap is a computationally cheap operation. In Shapely, the specific implementation is the STRtree, which stands for Sort-Tile-Recursive. You create one by passing an array of geometries into the constructor. Under the hood, this bulk-loads the index. The tree is built once and is immutable. You cannot add or remove items after it is created. You must prepare your entire dataset up front. Once the tree exists, you search it using the query method. You pass an input geometry, and the method returns a NumPy array of integer indices. These indices map directly to the original sequence of geometries you used to build the tree. If you input multiple geometries at once, query returns a two-dimensional array, pairing the indices of your inputs with the matching indices in the tree. Here is the key insight. The base index of an STRtree only checks bounding boxes, not exact geometries. If you query a complex polygon, the tree uses the square envelope of that polygon to find matches. This gives you a shortlist of candidates containing false positives. To get an exact answer, Shapely allows you to pass a spatial predicate, like intersects or contains, directly into the query method as an argument. The tree first uses the cheap bounding boxes to find candidates, then automatically evaluates your exact predicate on that shortlist. This keeps the heavy mathematical lifting down in the compiled C code. Beyond overlaps, the STRtree excels at proximity. Returning to our coffee shop scenario, you do not need to check envelopes manually. You use the query nearest method. You pass your user's GPS point, and the tree traverses its hierarchical boxes to find the absolute closest geometry without measuring against all ten thousand. It immediately returns the index of the nearest shop. You can also specify a parameter to return an array of the top N nearest neighbors, or set a maximum distance threshold to limit the search radius entirely. There is also a distinct nearest method, which compares two separate trees against each other to find the closest matching pairs of geometries between them. The actual geometry operations are always the most expensive part of spatial analysis. The STRtree acts as a ruthless filter. It prevents you from wasting CPU cycles on geometries that are nowhere near each other, providing a tiny, manageable shortlist so your exact mathematical checks only run when they absolutely have to. This wraps up our series on Shapely. I highly encourage you to explore the official documentation and try these indexing concepts hands-on in your own datasets. If you have workflows or tools you want covered in future series, visit devstories dot eu and drop us a suggestion. Thanks for tuning in. Until next time!