{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Geometric objects - Spatial data model\n", "\n", "**Sources:**\n", "\n", "These materials are partly based on [Shapely-documentation](http://toblerity.org/shapely/manual.html) and [Westra\n", "E. (2013), Chapter 3](https://www.packtpub.com/application-development/python-geospatial-development-second-edition).\n", "\n", "## Overview of geometric objects and Shapely -module\n", "\n", "![Spatial data model](img/SpatialDataModel.PNG)\n", "*Fundamental geometric objects that can be used in Python with* [Shapely](http://toblerity.org/shapely/manual.html) *module*\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. Python has a specific module called [Shapely](http://toblerity.org/shapely/manual.html) that can be used to create and work with `Geometric Objects`. There are many useful functionalities that you can do with Shapely such as:\n", "\n", "- Create a `Line` or `Polygon` from a `Collection` of `Point` -geometries\n", "- Calculate areas/length/bounds etc. of input geometries\n", "- Conduct geometric operations based on the input geometries such as `Union`, `Difference`, `Distance` etc.\n", "- Conduct spatial queries between geometries such `Intersects`, `Touches`, `Crosses`, `Within` etc.\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" ] }, { "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 }, "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)\n", "\n", "# What is the type of the point?\n", "point_type = type(point1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Let's see what the variables look like" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (2.2 4.2)\n", "POINT Z (9.26 -2.456 0.57)\n", "\n" ] } ], "source": [ "print(point1)\n", "print(point3D)\n", "print(type(point1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the type of the point is shapely's Point which is represented in a specific format that is based on\n", "[GEOS](https://trac.osgeo.org/geos) C++ library that is one of the standard libraries in GIS. It runs under the hood e.g. in [QGIS](http://www.qgis.org/en/site/). 3D-point can be recognized from the capital Z -letter in front of the coordinates." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Point attributes and functions\n", "\n", "Point -object has some built-in attributes that can be accessed and also some useful functionalities. One of the most useful ones are the ability to extract the coordinates of a Point and calculate the Euclidian distance between points.\n", "\n", "- Extracting the coordinates of a Point can be done in a couple of different ways:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "shapely.coords.CoordinateSequence" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get the coordinates\n", "point_coords = point1.coords\n", "\n", "# What is the type of this?\n", "type(point_coords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see, the data type of our `point_coords` variable is a Shapely CoordinateSequence.\n", "\n", "- Let's see how we can get out the actual coordinates from this object:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "xy variable:\n", " (array('d', [2.2]), array('d', [4.2])) \n", "\n", "x variable:\n", " 2.2 \n", "\n", "y variable:\n", " 4.2\n" ] } ], "source": [ "# Get x and y coordinates\n", "xy = point_coords.xy\n", "\n", "# Get only x coordinates of Point1\n", "x = point1.x\n", "\n", "# Whatabout y coordinate?\n", "y = point1.y\n", "\n", "# Print out\n", "print(\"xy variable:\\n\", xy, \"\\n\")\n", "print(\"x variable:\\n\", x, \"\\n\")\n", "print(\"y variable:\\n\", y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see from above the `xy` -variable contains a tuple where x and y coordinates are stored inside numpy arrays.\n", "Using the attributes `point1.x` and `point1.y` it is possible to get the coordinates directly as plain decimal numbers.\n", "\n", "- It is also possible to calculate the distance between points which can be useful in many applications. The returned distance is based on the projection of the points (e.g. degrees in WGS84, meters in UTM):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Distance between the points is 29.72 decimal degrees\n" ] } ], "source": [ "# Calculate the distance between point1 and point2\n", "point_dist = point1.distance(point2)\n", "\n", "print(\"Distance between the points is {0:.2f} decimal degrees\".format(point_dist))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## LineString\n", "\n", "\n", "Creating a LineString -object is fairly similar to how Point is created. \n", "\n", "- Now instead using a single coordinate-tuple we can construct the line using either a list of shapely Point -objects or pass coordinate-tuples:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "line variable: \n", " LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456) \n", "\n", "line2 variable: \n", " LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456) \n", "\n", "type of the line: \n", " \n" ] } ], "source": [ "# Create a LineString from our Point objects\n", "line = LineString([point1, point2, point3])\n", "\n", "# It is also possible to use coordinate tuples having the same outcome\n", "line2 = LineString([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])\n", "\n", "# Print the results\n", "print(\"line variable: \\n\", line, \"\\n\")\n", "print(\"line2 variable: \\n\", line2, \"\\n\")\n", "print(\"type of the line: \\n\", type(line))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see from above, the `line` -variable constitutes of multiple coordinate-pairs and the type of the data is shapely LineString." ] }, { "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`\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array('d', [2.2, 7.2, 9.26]), array('d', [4.2, -25.1, -2.456]))\n" ] } ], "source": [ "# Get x and y coordinates of the line\n", "lxy = line.xy\n", "\n", "print(lxy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see, the coordinates are again stored as a numpy arrays where first array includes all x-coordinates and the second one all the y-coordinates respectively.\n", "\n", "- We can extract only x or y coordinates by referring to those arrays as follows\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "line_x:\n", " array('d', [2.2, 7.2, 9.26]) \n", "\n", "line_y:\n", " array('d', [4.2, -25.1, -2.456])\n" ] } ], "source": [ "# Extract x coordinates\n", "line_x = lxy[0]\n", "\n", "# Extract y coordinates straight from the LineObject by referring to a array at index 1\n", "line_y = line.xy[1]\n", "\n", "print('line_x:\\n', line_x, '\\n')\n", "\n", "print('line_y:\\n', line_y)" ] }, { "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": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length of our line: 52.46\n", "Centroid of our line: POINT (6.229961354035622 -11.89241115757239)\n", "Type of the centroid: \n" ] } ], "source": [ "# Get the lenght of the line\n", "l_length = line.length\n", "\n", "# Get the centroid of the line\n", "l_centroid = line.centroid\n", "\n", "# What type is the centroid?\n", "centroid_type = type(l_centroid)\n", "\n", "# Print the outputs\n", "print(\"Length of our line: {0:.2f}\".format(l_length))\n", "print(\"Centroid of our line: \", l_centroid)\n", "print(\"Type of the centroid:\", centroid_type)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nice! These are already fairly useful information for many different GIS tasks, and we didn't even calculate anything yet! These attributes are built-in in every LineString object that is created. \n", "\n", "Notice that the centroid that is returned is `Point` -object that has its own functions as was described earlier." ] }, { "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 coordinate-tuples as input. \n", "\n", "- Polygon needs at least three coordinate-tuples (that basically forms a triangle):\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "poly: POLYGON ((2.2 4.2, 7.2 -25.1, 9.26 -2.456, 2.2 4.2))\n", "poly2: POLYGON ((2.2 4.2, 7.2 -25.1, 9.26 -2.456, 2.2 4.2))\n", "Geometry type as text: Polygon\n", "Geometry how Python shows it: \n" ] } ], "source": [ "# Create a Polygon from the coordinates\n", "poly = Polygon([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])\n", "\n", "# We can also use our previously created Point objects (same outcome)\n", "# --> notice that Polygon object requires x,y coordinates as input\n", "poly2 = Polygon([[p.x, p.y] for p in [point1, point2, point3]])\n", "\n", "# Geometry type can be accessed as a String\n", "poly_type = poly.geom_type\n", "\n", "# Using the Python's type function gives the type in a different format\n", "poly_type2 = type(poly)\n", "\n", "# Let's see how our Polygon looks like\n", "print('poly:', poly)\n", "print('poly2:', poly2)\n", "print(\"Geometry type as text:\", poly_type)\n", "print(\"Geometry how Python shows it:\", poly_type2)" ] }, { "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. \n", "\n", "As the help of Polygon -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 inside" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Let's create a bounding box of the world and make a whole in it\n", "\n", "# First we define our exterior\n", "world_exterior = [(-180, 90), (-180, -90), (180, -90), (180, 90)]\n", "\n", "# Let's create a single big hole where we leave ten decimal degrees at the boundaries of the world\n", "# Notice: there could be multiple holes, thus we need to provide a list of holes\n", "hole = [[(-170, 80), (-170, -80), (170, -80), (170, 80)]]\n", "\n", "# World without a hole\n", "world = Polygon(shell=world_exterior)\n", "\n", "# Now we can construct our Polygon with the hole inside\n", "world_has_a_hole = Polygon(shell=world_exterior, holes=hole)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Let's see what we have now:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "world: POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90))\n", "world_has_a_hole: POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90), (-170 80, -170 -80, 170 -80, 170 80, -170 80))\n", "type: \n" ] } ], "source": [ "print('world:', world)\n", "print('world_has_a_hole:', world_has_a_hole)\n", "print('type:', type(world_has_a_hole))" ] }, { "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`.\n", "\n", "- Here, we can see a few of the available attributes and how to access them:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Poly centroid: POINT (-0 -0)\n", "Poly Area: 64800.0\n", "Poly Bounding Box: (-180.0, -90.0, 180.0, 90.0)\n", "Poly Exterior: LINEARRING (-180 90, -180 -90, 180 -90, 180 90, -180 90)\n", "Poly Exterior Length: 1080.0\n" ] } ], "source": [ "# Get the centroid of the Polygon\n", "world_centroid = world.centroid\n", "\n", "# Get the area of the Polygon\n", "world_area = world.area\n", "\n", "# Get the bounds of the Polygon (i.e. bounding box)\n", "world_bbox = world.bounds\n", "\n", "# Get the exterior of the Polygon\n", "world_ext = world.exterior\n", "\n", "# Get the length of the exterior\n", "world_ext_length = world_ext.length\n", "\n", "# Print the outputs\n", "print(\"Poly centroid: \", world_centroid)\n", "print(\"Poly Area: \", world_area)\n", "print(\"Poly Bounding Box: \", world_bbox)\n", "print(\"Poly Exterior: \", world_ext)\n", "print(\"Poly Exterior Length: \", world_ext_length)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see above, it is again fairly straightforward to access different attributes from the `Polygon` -object. Notice, that the extrerior lenght is given here with decimal degrees because we passed latitude and longitude coordinates into our Polygon. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geometry collections (optional)\n", "\n", "\n", "In some occassions it is useful to store e.g. multiple lines or polygons under a single feature (i.e. a single row in a Shapefile represents more than one line or polygon object). 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. These collections are not computationally significant, but are useful for modeling certain kinds of features. A Y-shaped line feature\n", "(such as road), or multiple polygons (e.g. islands on a like), can be presented nicely as a whole by a using MultiLineString or MultiPolygon accordingly. Creating and visualizing a minimum [bounding box](https://en.wikipedia.org/wiki/Minimum_bounding_box) e.g. around your data points is a really useful function for many purposes (e.g. trying to understand the extent of your data), here we demonstrate how to create one using Shapely.\n", "\n", "- Geometry collections can be constructed in a following manner:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MultiPoint: MULTIPOINT (2.2 4.2, 7.2 -25.1, 9.26 -2.456)\n", "MultiLine: MULTILINESTRING ((2.2 4.2, 7.2 -25.1), (7.2 -25.1, 9.26 -2.456))\n", "Bounding box: POLYGON ((180 -90, 180 90, 0 90, 0 -90, 180 -90))\n", "MultiPoly: 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": [ "# Import collections of geometric objects + bounding box\n", "from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon, box\n", "\n", "# 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", "# MultiPolygon can be done in a similar manner\n", "# Let's divide our world into western and eastern hemispheres with a hole on the western hemisphere\n", "# --------------------------------------------------------------------------------------------------\n", "\n", "# 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", "# Let's create the Polygon of our Eastern hemisphere polygon using bounding box\n", "# For bounding box we need to specify the lower-left corner coordinates and upper-right coordinates\n", "min_x, min_y = 0, -90\n", "max_x, max_y = 180, 90\n", "\n", "# Create the polygon using box() function\n", "east_poly_box = box(minx=min_x, miny=min_y, maxx=max_x, maxy=max_y)\n", "\n", "# 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_box])\n", "\n", "# Print outputs\n", "print(\"MultiPoint:\", multi_point)\n", "print(\"MultiLine: \", multi_line)\n", "print(\"Bounding box: \", east_poly_box)\n", "print(\"MultiPoly: \", multi_poly)\n" ] }, { "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", "### Geometry collection -objects' attributes and functions\n", "\n", "- We can also get many useful attributes from those objects such as `Convex Hull`:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Convex hull of the points: POLYGON ((7.2 -25.1, 2.2 4.2, 9.26 -2.456, 7.2 -25.1))\n", "Number of lines in MultiLineString: 2\n", "Area of our MultiPolygon: 39200.0\n", "Area of our Western Hemisphere polygon: 6800.0\n", "Is polygon valid?: False\n" ] } ], "source": [ "# Convex Hull of our MultiPoint --> https://en.wikipedia.org/wiki/Convex_hull\n", "convex = multi_point.convex_hull\n", "\n", "# How many lines do we have inside our MultiLineString?\n", "lines_count = len(multi_line)\n", "\n", "# Let's calculate the area of our MultiPolygon\n", "multi_poly_area = multi_poly.area\n", "\n", "# We can also access different items inside our geometry collections. We can e.g. access a single polygon from\n", "# our MultiPolygon -object by referring to the index\n", "\n", "# Let's calculate the area of our Western hemisphere (with a hole) which is at index 0\n", "west_area = multi_poly[0].area\n", "\n", "# We can check if we have a \"valid\" MultiPolygon. MultiPolygon is thought as valid if the individual polygons \n", "# does notintersect with each other. Here, because the polygons have a common 0-meridian, we should NOT have \n", "# a valid polygon. This can be really useful information when trying to find topological errors from your data\n", "valid = multi_poly.is_valid\n", "\n", "# Print outputs\n", "print(\"Convex hull of the points: \", convex)\n", "print(\"Number of lines in MultiLineString:\", lines_count)\n", "print(\"Area of our MultiPolygon:\", multi_poly_area)\n", "print(\"Area of our Western Hemisphere polygon:\", west_area)\n", "print(\"Is polygon valid?: \", valid)" ] }, { "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. There are also some extra features available such as **is_valid** -attribute that tells if the polygons or lines intersect with each other." ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [default]", "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.5.2" } }, "nbformat": 4, "nbformat_minor": 2 }