{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Shapely and geometric objects\n", "\n", "In this lesson, you will learn how to create and manipulate geometries in Python using the [Shapely Python Package](https://shapely.readthedocs.io/en/stable/manual.html).\n", "\n", "**Sources:**\n", "\n", "These materials are partly based on [Shapely-documentation](https://shapely.readthedocs.io/en/stable/manual.html) and [Westra\n", "E. (2013), Chapter 3](https://www.packtpub.com/application-development/python-geospatial-development-second-edition).\n", "\n", "## Spatial data model\n", "\n", "![Spatial data model](img/SpatialDataModel.PNG)\n", "\n", "*Fundamental geometric objects that can be used in Python with* [Shapely](https://shapely.readthedocs.io/en/stable/manual.html).\n", "\n", "The most fundamental geometric objects are `Points`, `Lines` and `Polygons` which are the basic ingredients when working with spatial data in vector format. \n", "Python has a specific module called [Shapely](https://shapely.readthedocs.io/en/stable/manual.html) for doing various geometric operations. Basic knowledge of using Shapely is fundamental for understanding how geometries are stored and handled in GeoPandas.\n", "\n", "**Geometric objects consist of coordinate tuples where:**\n", "\n", "- `Point` -object represents a single point in space. Points can be either two-dimensional (x, y) or three dimensional (x, y, z).\n", "- `LineString` -object (i.e. a line) represents a sequence of points joined together to form a line. Hence, a line consist of a list of at least two coordinate tuples\n", "- `Polygon` -object represents a filled area that consists of a list of at least three coordinate tuples that forms the outerior ring and a (possible) list of hole polygons.\n", "\n", "**It is also possible to have a collection of geometric objects (e.g. Polygons with multiple parts):**\n", "\n", "- `MultiPoint` -object represents a collection of points and consists of a list of coordinate-tuples\n", "- `MultiLineString` -object represents a collection of lines and consists of a list of line-like sequences\n", "- `MultiPolygon` -object represents a collection of polygons that consists of a list of polygon-like sequences that construct from exterior ring and (possible) hole list tuples\n", "\n", "**Useful attributes and methods in Shapely include:**\n", "\n", "- Creating lines and polygons based on a collection of point objects.\n", "- Calculating areas/length/bounds etc. of input geometries\n", "- Conducting geometric operations based on the input geometries such as `union`, `difference`, `distance` etc.\n", "- Conducting spatial queries between geometries such as `intersects`, `touches`, `crosses`, `within` etc.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Tuple**\n", "\n", "[Tuple](https://docs.python.org/3/tutorial/datastructures.html#tuples-and-sequences) is a Python data structure that consists of a number of values separated by commas. Coordinate pairs are often represented as a tuple. For example:\n", "\n", "```\n", "(60.192059, 24.945831)\n", "``` \n", "\n", "Tuples belong to [sequence data types](https://docs.python.org/3/library/stdtypes.html#typesseq) in Python. Other sequence data types are lists and ranges. Tuples have many similarities with lists and ranges, but they are often used for different purposes. The main difference between tuples and lists is that tuples are [immutable](https://docs.python.org/3/glossary.html#term-immutable), which means that the contents of a tuple cannot be altered (while lists are mutable; you can, for example, add and remove values from lists).\n", "
\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Point\n", "\n", "Creating point is easy, you pass x and y coordinates into `Point()` -object (+ possibly also z -coordinate):" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Import necessary geometric objects from shapely module\n", "from shapely.geometry import Point, LineString, Polygon\n", "\n", "# Create Point geometric object(s) with coordinates\n", "point1 = Point(2.2, 4.2)\n", "point2 = Point(7.2, -25.1)\n", "point3 = Point(9.26, -2.456)\n", "point3D = Point(9.26, -2.456, 0.57)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see what these variables now contain: " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "point1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we see here, Jupyter notebook is able to display the shape directly on the screen.\n", "\n", "We can use the print statement to get information about the actual definition of these objects:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (2.2 4.2)\n", "POINT Z (9.26 -2.456 0.57)\n" ] } ], "source": [ "print(point1)\n", "print(point3D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3D-point can be recognized from the capital Z -letter in front of the coordinates.\n", "\n", "Let's also check the data type of a point:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "shapely.geometry.point.Point" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(point1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the type of the point is shapely's Point. The point object is represented in a specific format based on\n", "[GEOS](https://trac.osgeo.org/geos) C++ library that is one of the standard libraries behind various Geographic Information Systems. It runs under the hood e.g. in [QGIS](http://www.qgis.org/en/site/). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Point attributes and functions\n", "\n", "Points and other shapely objects have useful built-in [attributes and methods](https://shapely.readthedocs.io/en/stable/manual.html#general-attributes-and-methods). Using the available attributes, we can for example extract the coordinate values of a Point and calculate the Euclidian distance between points." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`geom_type` attribute contains information about the geometry type of the Shapely object:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'Point'" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "point1.geom_type" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extracting the coordinates of a Point can be done in a couple of different ways:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`coords` attribute contains the coordinate information as a `CoordinateSequence` which is another data type related to Shapely." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(2.2, 4.2)]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get xy coordinate tuple\n", "list(point1.coords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we have a coordinate tuple inside a list. Using the attributes `x` and `y` it is possible to get the coordinates directly as plain decimal numbers." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Read x and y coordinates separately\n", "x = point1.x\n", "y = point1.y" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.2 4.2\n" ] } ], "source": [ "print( x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to calculate the distance between two objects using the [distance](https://shapely.readthedocs.io/en/stable/manual.html#object.distance) method. In our example the distance is calculated in a cartesian coordinate system. When working with real GIS data the distance is based on the used coordinate reference system. always check what is the unit of measurement (for example, meters) in the coordinate reference system you are using.\n", "\n", "Let's calculate the distance between `point1` and `point2`:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (2.2 4.2)\n", "POINT (7.2 -25.1)\n" ] } ], "source": [ "# Check input data\n", "print(point1)\n", "print(point2)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Distance between the points is 29.72 units\n" ] } ], "source": [ "# Calculate the distance between point1 and point2\n", "dist = point1.distance(point2)\n", "\n", "# Print out a nicely formatted info message\n", "print(\"Distance between the points is {0:.2f} units\".format(dist))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## LineString\n", "\n", "\n", "Creating LineString -objects is fairly similar to creating Shapely Points. \n", "\n", "Now instead using a single coordinate-tuple we can construct the line using either a list of shapely Point -objects or pass the points as coordinate-tuples:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Create a LineString from our Point objects\n", "line = LineString([point1, point2, point3])" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# It is also possible to produce the same outcome using coordinate tuples\n", "line2 = LineString([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check if lines are identical\n", "line == line2 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see how our line looks like: " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "line" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456)\n" ] } ], "source": [ "print(line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see from above, the `line` -variable constitutes of multiple coordinate-pairs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check also the data type:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "shapely.geometry.linestring.LineString" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check data type of the line object\n", "type(line)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'LineString'" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check geometry type of the line object\n", "line.geom_type" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### LineString attributes and functions\n", "\n", "\n", "`LineString` -object has many useful built-in attributes and functionalities. It is for instance possible to extract the coordinates or the length of a LineString (line), calculate the centroid of the line, create points along the line at specific distance, calculate the closest distance from a line to specified Point and simplify the geometry. See full list of functionalities from [Shapely documentation](http://toblerity.org/shapely/manual.html). Here, we go through a few of them.\n", "\n", "We can extract the coordinates of a LineString similarly as with `Point`" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "[(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get xy coordinate tuples\n", "list(line.coords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we have a list of coordinate tuples (x,y) inside a list.\n", "\n", "If you would need to access all x-coordinates or all y-coordinates of the line, you can do it directly using the `xy` attribute: " ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Extract x and y coordinates separately\n", "xcoords = list(line.xy[0])\n", "ycoords = list(line.xy[1])" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2.2, 7.2, 9.26]\n", "[4.2, -25.1, -2.456]\n" ] } ], "source": [ "print(xcoords)\n", "print(ycoords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is possible to retrieve specific attributes such as lenght of the line and center of the line (centroid) straight from the LineString object itself:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length of our line: 52.46 units\n" ] } ], "source": [ "# Get the lenght of the line\n", "l_length = line.length\n", "print(\"Length of our line: {0:.2f} units\".format(l_length))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (6.229961354035622 -11.89241115757239)\n" ] } ], "source": [ "# Get the centroid of the line\n", "print(line.centroid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the centroid of the line is again a Shapely Point object. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Polygon\n", "\n", "\n", "Creating a `Polygon` -object continues the same logic of how `Point` and `LineString` were created but Polygon object only accepts a sequence of coordinates as input. \n", "\n", "Polygon needs **at least three coordinate-tuples** (three points are reguired to form a surface):" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Create a Polygon from the coordinates\n", "poly = Polygon([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also use information from the Shapely Point objects created earlier, but we can't use the point objects directly. Instead, we need to get information of the x,y coordinate pairs as a sequence. We can achieve this by using a list comprehension." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# Create a Polygon based on information from the Shapely points\n", "poly2 = Polygon([[p.x, p.y] for p in [point1, point2, point3]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to understand what just happened, let's check what the list comprehension produces:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[2.2, 4.2], [7.2, -25.1], [9.26, -2.456]]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[[p.x, p.y] for p in [point1, point2, point3]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This list of lists was passed as input for creating the Polygon." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check that polygon objects created using two different approaches are identical\n", "poly == poly2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see how our Polygon looks like" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "poly" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((2.2 4.2, 7.2 -25.1, 9.26 -2.456, 2.2 4.2))\n" ] } ], "source": [ "print(poly)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that `Polygon` representation has double parentheses around the coordinates (i.e. `POLYGON (())` ). This is because Polygon can also have holes inside of it. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check also the data type:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "shapely.geometry.polygon.Polygon" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Data type\n", "type(poly)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'Polygon'" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Geometry type\n", "poly.geom_type" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "# Check the help for Polygon objects:\n", "#help(Polygon)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "As the help of [Polygon](https://shapely.readthedocs.io/en/stable/manual.html#polygons) -object tells, a Polygon can be constructed using exterior coordinates and interior coordinates (optional) where the interior coordinates creates a hole inside the Polygon:\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "Help on Polygon in module shapely.geometry.polygon object:\n", " class Polygon(shapely.geometry.base.BaseGeometry)\n", " | A two-dimensional figure bounded by a linear ring\n", " |\n", " | A polygon has a non-zero area. It may have one or more negative-space\n", " | \"holes\" which are also bounded by linear rings. If any rings cross each\n", " | other, the feature is invalid and operations on it may fail.\n", " |\n", " | Attributes\n", " | ----------\n", " | exterior : LinearRing\n", " | The ring which bounds the positive space of the polygon.\n", " | interiors : sequence\n", " | A sequence of rings which bound all existing holes.\n", " \n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see how we can create a `Polygon` with a hole:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "# Define the outer border\n", "border = [(-180, 90), (-180, -90), (180, -90), (180, 90)]" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90))\n" ] } ], "source": [ "# Outer polygon\n", "world = Polygon(shell=border)\n", "print(world)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "world" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "# Let's create a single big hole where we leave ten units at the boundaries\n", "# Note: there could be multiple holes, so we need to provide list of coordinates for the hole inside a list\n", "hole = [[(-170, 80), (-170, -80), (170, -80), (170, 80)]]" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90), (-170 80, -170 -80, 170 -80, 170 80, -170 80))\n" ] } ], "source": [ "# Now we can construct our Polygon with the hole inside\n", "frame = Polygon(shell=border, holes=hole)\n", "print(frame)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see what we have now:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "frame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see the `Polygon` has now two different tuples of coordinates. The first one represents the **outerior** and the second one represents the **hole** inside of the Polygon." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Polygon attributes and functions\n", "\n", "\n", "We can again access different attributes directly from the `Polygon` object itself that can be really useful for many analyses, such as `area`, `centroid`, `bounding box`, `exterior`, and `exterior-length`. See a full list of methods in the [Shapely User Manual](https://shapely.readthedocs.io/en/stable/manual.html#the-shapely-user-manual).\n", "\n", "Here, we can see a few of the available attributes and how to access them:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Polygon centroid: POINT (-0 -0)\n", "Polygon Area: 64800.0\n", "Polygon Bounding Box: (-180.0, -90.0, 180.0, 90.0)\n", "Polygon Exterior: LINEARRING (-180 90, -180 -90, 180 -90, 180 90, -180 90)\n", "Polygon Exterior Length: 1080.0\n" ] } ], "source": [ "# Print the outputs\n", "print(\"Polygon centroid: \", world.centroid)\n", "print(\"Polygon Area: \", world.area)\n", "print(\"Polygon Bounding Box: \", world.bounds)\n", "print(\"Polygon Exterior: \", world.exterior)\n", "print(\"Polygon Exterior Length: \", world.exterior.length)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see above, it is again fairly straightforward to access different attributes from the `Polygon` -object. Note that distance metrics will make more sense when we start working with data in a projected coordinate system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Check your understanding\n", "\n", "Plot these shapes using Shapely!\n", "\n", "- **Pentagon**, example coords: `(30, 2.01), (31.91, 0.62), (31.18, -1.63), (28.82, -1.63), (28.09, 0.62)` \n", "- **Triangle** \n", "- **Square** \n", "- **Cicrle** \n" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Pentagon - Coordinates borrowed from this thread: https://tex.stackexchange.com/questions/179843/make-a-polygon-with-automatically-labelled-nodes-according-to-their-coordinates\n", "Polygon([(30, 2.01), (31.91, 0.62), (31.18, -1.63), (28.82, -1.63), (28.09, 0.62)])" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Triangle\n", "Polygon([(0,0), (2,4), (4,0)])" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Square\n", "Polygon([(0,0), (0,4), (4,4), (4,0)])" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Circle (using a buffer around a point)\n", "point = Point((0,0))\n", "point.buffer(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geometry collections (optional)\n", "\n", "\n", "In some occassions it is useful to store multiple geometries (for example, several points or several polygons) in a single feature. A practical example would be a country that is composed of several islands. In such case, all these polygons share the same attributes on the country-level and it might be reasonable to store that country as geometry collection that contains all the polygons. The attribute table would then contain one row of information with country-level attributes, and the geometry related to those attributes would represent several polygon. \n", "\n", "In Shapely, collections of points are implemented by using a MultiPoint -object, collections of curves by using a MultiLineString -object, and collections of surfaces by a MultiPolygon -object. " ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "# Import constructors for creating geometry collections\n", "from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start by creating MultiPoint and MultilineString objects:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MULTIPOINT (2.2 4.2, 7.2 -25.1, 9.26 -2.456)\n", "MULTILINESTRING ((2.2 4.2, 7.2 -25.1), (7.2 -25.1, 9.26 -2.456))\n" ] } ], "source": [ "# Create a MultiPoint object of our points 1,2 and 3\n", "multi_point = MultiPoint([point1, point2, point3])\n", "\n", "# It is also possible to pass coordinate tuples inside\n", "multi_point2 = MultiPoint([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])\n", "\n", "# We can also create a MultiLineString with two lines\n", "line1 = LineString([point1, point2])\n", "line2 = LineString([point2, point3])\n", "multi_line = MultiLineString([line1, line2])\n", "\n", "# Print object definitions\n", "print(multi_point)\n", "print(multi_line)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "multi_point" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "multi_line" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MultiPolygons are constructed in a similar manner. Let's create a bounding box for \"the world\" by combinin two separate polygons that represent the western and eastern hemispheres. " ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((-180 90, -180 -90, 0 -90, 0 90, -180 90), (-170 80, -170 -80, -10 -80, -10 80, -170 80))\n" ] } ], "source": [ "# Let's create the exterior of the western part of the world\n", "west_exterior = [(-180, 90), (-180, -90), (0, -90), (0, 90)]\n", "\n", "# Let's create a hole --> remember there can be multiple holes, thus we need to have a list of hole(s). \n", "# Here we have just one.\n", "west_hole = [[(-170, 80), (-170, -80), (-10, -80), (-10, 80)]]\n", "\n", "# Create the Polygon\n", "west_poly = Polygon(shell=west_exterior, holes=west_hole)\n", "\n", "# Print object definition\n", "print(west_poly)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "west_poly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Shapely also has a tool for creating [a bounding box](https://en.wikipedia.org/wiki/Minimum_bounding_box) based on minimum and maximum x and y coordinates. Instead of using the Polygon constructor, let's use the [box](https://shapely.readthedocs.io/en/stable/manual.html#shapely.geometry.box) constructor for creating the polygon: " ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import box" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((180 -90, 180 90, 0 90, 0 -90, 180 -90))\n" ] } ], "source": [ "# Specify the bbox extent (lower-left corner coordinates and upper-right corner coordinates)\n", "min_x, min_y = 0, -90\n", "max_x, max_y = 180, 90\n", "\n", "# Create the polygon using Shapely\n", "east_poly = box(minx=min_x, miny=min_y, maxx=max_x, maxy=max_y)\n", "\n", "# Print object definition\n", "print(east_poly)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "east_poly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can combine the two polygons into a MultiPolygon:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MULTIPOLYGON (((-180 90, -180 -90, 0 -90, 0 90, -180 90), (-170 80, -170 -80, -10 -80, -10 80, -170 80)), ((180 -90, 180 90, 0 90, 0 -90, 180 -90)))\n" ] } ], "source": [ "# Let's create our MultiPolygon. We can pass multiple Polygon -objects into our MultiPolygon as a list\n", "multi_poly = MultiPolygon([west_poly, east_poly])\n", "\n", "# Print object definition\n", "print(multi_poly)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "multi_poly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the outputs are similar to the basic geometric objects that we created previously but now these objects contain multiple features of those points, lines or polygons.\n", "\n", "### Convex hull and envelope\n", "\n", "Convex hull refers to the smalles possible polygon that contains all objects in a collection. Alongside with the minimum bounding box, convex hull is a useful shape when aiming to describe the extent of your data. \n", "\n", "Let's create a convex hull around our multi_point object:" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check input geometry\n", "multi_point" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Convex Hull (smallest polygon around the geometry collection)\n", "multi_point.convex_hull" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Envelope (smalles rectangular polygon around the geometry collection): \n", "multi_point.envelope" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other useful attributes \n", "lenght of the geometry collection:" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of objects in our MultiLine: 2\n", "Number of objects in our MultiPolygon: 2\n" ] } ], "source": [ "print(\"Number of objects in our MultiLine:\", len(multi_line))\n", "print(\"Number of objects in our MultiPolygon:\", len(multi_poly))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Area:" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Area of our MultiPolygon: 39200.0\n", "Area of our Western Hemisphere polygon: 6800.0\n" ] } ], "source": [ "# Print outputs:\n", "print(\"Area of our MultiPolygon:\", multi_poly.area)\n", "print(\"Area of our Western Hemisphere polygon:\", multi_poly[0].area)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the above we can see that MultiPolygons have exactly the same attributes available as single geometric objects but now the information such as area calculates the area of **ALL** of the individual -objects combined. We can also access individual objects inside the geometry collections using indices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can check if we have a \"valid\" MultiPolygon. MultiPolygon is thought as valid if the individual polygons does notintersect with each other. \n", "Here, because the polygons have a common 0-meridian, we should NOT have a valid polygon. We can check the validity of an object from the **is_valid** -attribute that tells if the polygons or lines intersect with each other. This can be really useful information when trying to find topological errors from your data:" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Is polygon valid?: False\n" ] } ], "source": [ "print(\"Is polygon valid?: \", multi_poly.is_valid)" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" } }, "nbformat": 4, "nbformat_minor": 4 }