{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lesson1: Geometric Objects - Spatial Data Model\n", "\n", "- https://kodu.ut.ee/~kmoch/geopython2019/L1/Geometric-Objects.html" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import necessary geometric objects from shapely module\n", "from shapely.geometry import Point, LineString, Polygon" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Create Point geometric object(s) with coordinates\n", "point1 = Point(2.2, 4.2)\n", "\n", "point2 = Point(7.2, -25.1)\n", "\n", "point3 = Point(9.26, -2.456)\n", "\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": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POINT (2.2 4.2)\n", "POINT Z (9.26 -2.456 0.57)\n", "\n" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(point1)\n", "\n", "print(point3D)\n", "\n", "print(type(point1))\n", "\n", "display(point1)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "shapely.coords.CoordinateSequence" ] }, "execution_count": 4, "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": "code", "execution_count": 5, "metadata": {}, "outputs": [], "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" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array('d', [2.2]), array('d', [4.2]))\n", "2.2\n", "4.2\n" ] } ], "source": [ "print(xy)\n", "\n", "print(x)\n", "\n", "print(y)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "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": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3306.1044 for equatorial radius in km\n", "3294.7004 for polar radius in km\n" ] } ], "source": [ "# law of cosines - determines the great-circle distance between two points on a \n", "# sphere given their longitudes and latitudes based on \"basic math\"\n", "import math\n", "\n", "distance = math.acos(math.sin(math.radians(point1.y))*math.sin(math.radians(point2.y))+\n", " math.cos(math.radians(point1.y))*math.cos(math.radians(point2.y))*\n", " math.cos(math.radians(point2.x)-math.radians(point1.x)))*6378\n", "\n", "print( \"{0:8.4f} for equatorial radius in km\".format(distance))\n", "\n", "distance = math.acos(math.sin(math.radians(point1.y))*math.sin(math.radians(point2.y))+\n", " math.cos(math.radians(point1.y))*math.cos(math.radians(point2.y))*\n", " math.cos(math.radians(point2.x)-math.radians(point1.x)))*6356\n", "\n", "print( \"{0:8.4f} for polar radius in km\".format(distance))\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3286.3538 for ellipsoid WGS84 in km\n" ] } ], "source": [ "# with pyproj\n", "import pyproj\n", "\n", "geod = pyproj.Geod(ellps='WGS84')\n", "\n", "angle1,angle2,distance = geod.inv(point1.x, point1.y, point2.x, point2.y)\n", "\n", "print (\"{0:8.4f} for ellipsoid WGS84 in km\".format(distance/1000))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "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)])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456)\n", "LINESTRING (2.2 4.2, 7.2 -25.1, 9.26 -2.456)\n" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(line)\n", "print(line2)\n", "type(line)\n", "display(line)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "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": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array('d', [2.2, 7.2, 9.26])\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", "\n", "print(line_y)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "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": "code", "execution_count": 15, "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", "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" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "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)\n", "print(poly2)\n", "print(\"Geometry type as text:\", poly_type)\n", "print(\"Geometry how Python shows it:\", poly_type2)\n", "\n", "display(poly)\n", "display(poly2)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "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": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90))\n", "POLYGON ((-180 90, -180 -90, 180 -90, 180 90, -180 90), (-170 80, -170 -80, 170 -80, 170 80, -170 80))\n" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(world)\n", "print(world_has_a_hole)\n", "type(world_has_a_hole)\n", "\n", "display(world)\n", "display(world_has_a_hole)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "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" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "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": [ "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": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Country or areaUrban AgglomerationLatitudeLongitudePopulation_2015Unnamed: 5
0JapanTokyo35.689500139.69171038001018NaN
1IndiaDelhi28.66667077.21667025703168NaN
2ChinaShanghai31.220000121.46000023740778NaN
3BrazilS?o Paulo-23.550000-46.64000021066245NaN
4IndiaMumbai (Bombay)19.07397572.88083821042538NaN
\n", "
" ], "text/plain": [ " Country or area Urban Agglomeration Latitude Longitude Population_2015 \\\n", "0 Japan Tokyo 35.689500 139.691710 38001018 \n", "1 India Delhi 28.666670 77.216670 25703168 \n", "2 China Shanghai 31.220000 121.460000 23740778 \n", "3 Brazil S?o Paulo -23.550000 -46.640000 21066245 \n", "4 India Mumbai (Bombay) 19.073975 72.880838 21042538 \n", "\n", " Unnamed: 5 \n", "0 NaN \n", "1 NaN \n", "2 NaN \n", "3 NaN \n", "4 NaN " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import pandas as pd\n", "\n", "# make sure you have the correct path to your working file\n", "# e.g. 'L1/global-city-population-estimates.csv' if you saved the file in your working directory\n", "\n", "df = pd.read_csv('global-city-population-estimates.csv', \";\", encoding='latin1')\n", "\n", "# this option tells pandas to print up to 20 columns, typically a the print function will cut the output for better visibility\n", "# (depending on the size and dimension of the dataframe)\n", "pd.set_option('max_columns',20)\n", "\n", "# print(df.head(5))\n", "display(df.head(5))" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# we make a function, that takes a row object coming from Pandas. The single fields per row are addressed by their column name.\n", "def make_point(row):\n", " return Point(row['Longitude'], row['Latitude'])\n", "\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Country or areaUrban AgglomerationLatitudeLongitudePopulation_2015Unnamed: 5points
0JapanTokyo35.689500139.69171038001018NaNPOINT (139.69171 35.6895)
1IndiaDelhi28.66667077.21667025703168NaNPOINT (77.21666999999999 28.66667)
2ChinaShanghai31.220000121.46000023740778NaNPOINT (121.46 31.22)
3BrazilS?o Paulo-23.550000-46.64000021066245NaNPOINT (-46.64 -23.55)
4IndiaMumbai (Bombay)19.07397572.88083821042538NaNPOINT (72.880838 19.073975)
\n", "
" ], "text/plain": [ " Country or area Urban Agglomeration Latitude Longitude Population_2015 \\\n", "0 Japan Tokyo 35.689500 139.691710 38001018 \n", "1 India Delhi 28.666670 77.216670 25703168 \n", "2 China Shanghai 31.220000 121.460000 23740778 \n", "3 Brazil S?o Paulo -23.550000 -46.640000 21066245 \n", "4 India Mumbai (Bombay) 19.073975 72.880838 21042538 \n", "\n", " Unnamed: 5 points \n", "0 NaN POINT (139.69171 35.6895) \n", "1 NaN POINT (77.21666999999999 28.66667) \n", "2 NaN POINT (121.46 31.22) \n", "3 NaN POINT (-46.64 -23.55) \n", "4 NaN POINT (72.880838 19.073975) " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Go through every row, and make a point out of its lat and lon, by **apply**ing the function from above (downwards row by row -> axis=1)\n", "df['points'] = df.apply(make_point, axis=1)\n", "\n", "display(df.head(5))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "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" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MultiPoint: MULTIPOINT (2.2 4.2, 7.2 -25.1, 9.26 -2.456)\n" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "MultiLine: MULTILINESTRING ((2.2 4.2, 7.2 -25.1), (7.2 -25.1, 9.26 -2.456))\n" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(\"MultiPoint:\", multi_point)\n", "display(multi_point)\n", "\n", "print(\"MultiLine: \", multi_line)\n", "display(multi_line)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "# 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])" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bounding box: POLYGON ((180 -90, 180 90, 0 90, 0 -90, 180 -90))\n" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "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" ] }, { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(\"Bounding box: \", east_poly_box)\n", "display(east_poly_box)\n", "\n", "print(\"MultiPoly: \", multi_poly)\n", "display(multi_poly)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "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" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "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": [ "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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "geopy2019", "language": "python", "name": "geopy2019" }, "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.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }