{ "cells": [ { "cell_type": "markdown", "id": "9e87802f-2461-4371-b489-b9bbe58fc5cb", "metadata": {}, "source": [ "# Geometric data manipulations" ] }, { "cell_type": "markdown", "id": "1eb5a35f", "metadata": {}, "source": [ "In this section we will use the Helsinki Region Travel Time Matrix data that constist of 13231 statistical grid squares (250m x 250m) to demonstrate some of the most common geometry manipulation functions available in geopandas. \n", "\n", "As the geometries in GeoDataFrames are eventually Shapely objects, we can use all of Shapely's tools for geometry manipulation directly via geopandas.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "62278989-accb-4a59-b6ee-7d31f98f6e75", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 POLYGON ((382000.000 6697750.000, 381750.000 6...\n", "1 POLYGON ((382250.000 6697750.000, 382000.000 6...\n", "2 POLYGON ((382500.000 6697750.000, 382250.000 6...\n", "3 POLYGON ((382750.000 6697750.000, 382500.000 6...\n", "4 POLYGON ((381250.000 6697500.000, 381000.000 6...\n", "Name: geometry, dtype: geometry" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import geopandas as gpd\n", "from pathlib import Path\n", "\n", "input_folder = Path(\"../data/Helsinki\")\n", "fp = input_folder / \"TravelTimes_to_5975375_RailwayStation.shp\"\n", "\n", "data = gpd.read_file(fp)\n", "\n", "data.geometry.head()" ] }, { "cell_type": "code", "execution_count": 11, "id": "f32e4fb2-9e2c-44c1-bcef-5191b415aa28", "metadata": {}, "outputs": [], "source": [ "# Plot grid geometry" ] }, { "cell_type": "markdown", "id": "e5fbcdae-5a13-4239-9988-c00d6ef7d96e", "metadata": { "tags": [] }, "source": [ "## Centroid\n", "\n", "Extracting the centroid of geometric features is useful in a multitude of use cases. For example, the values in the Travel Time Matrix data set have originally been calculated from the center point of each grid square. We can find out the geometric centroid of each grid square in geopandas via the centroid-attribute:" ] }, { "cell_type": "code", "execution_count": 2, "id": "13e2a48e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 POINT (381875.000 6697875.000)\n", "1 POINT (382125.000 6697875.000)\n", "2 POINT (382375.000 6697875.000)\n", "3 POINT (382625.000 6697875.000)\n", "4 POINT (381125.000 6697625.000)\n", " ... \n", "13226 POINT (372875.000 6665625.000)\n", "13227 POINT (373125.000 6665625.000)\n", "13228 POINT (372375.000 6665375.000)\n", "13229 POINT (372625.000 6665375.000)\n", "13230 POINT (372875.000 6665375.000)\n", "Length: 13231, dtype: geometry" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Polygon centroids\n", "data.centroid" ] }, { "cell_type": "code", "execution_count": 3, "id": "ea10f0dd-f15b-4624-a93c-9ea30e37130e", "metadata": {}, "outputs": [], "source": [ "# plot centroid geometries" ] }, { "cell_type": "markdown", "id": "efef7987", "metadata": {}, "source": [ "## Unary union\n", "\n", "Extracting only the outlines of the set of grid squares is possible through creating a geometric union among all geometries in the GeoDataFrame/GeoSeries. This could be useful, for example, when delineating the outlines of a study area." ] }, { "cell_type": "code", "execution_count": 4, "id": "993def06", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Polygon centroids\n", "data.unary_union" ] }, { "cell_type": "markdown", "id": "7584c825-8be2-486b-87f8-13add35eee89", "metadata": {}, "source": [ "Note that the previous operation is identical to calling `data['geometry'].unary_union`." ] }, { "cell_type": "markdown", "id": "b36ec0ac", "metadata": {}, "source": [ "## Data extent\n", "\n", "Sometimes it is enough to describe the approximate extent of the data using a bounding polygon." ] }, { "cell_type": "markdown", "id": "83e4d86c-c9a0-4ee2-bdc0-72845faa4158", "metadata": {}, "source": [ "### Bounding box\n", "\n", "A minimum bounding rectangle, i.e. a bounding box or an envelope is the smallest rectangular polygon surrounding a geometric object. In a GeoDataFrame, the `envelope` attribute returns the bounding rectangle for each geometry:" ] }, { "cell_type": "code", "execution_count": 5, "id": "8c0513d4-03c8-42d1-917d-aaf43c5dab38", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 POLYGON ((381750.000 6697750.000, 382000.000 6...\n", "1 POLYGON ((382000.000 6697750.000, 382250.000 6...\n", "2 POLYGON ((382250.000 6697750.000, 382500.000 6...\n", "3 POLYGON ((382500.000 6697750.000, 382750.000 6...\n", "4 POLYGON ((381000.000 6697500.000, 381250.000 6...\n", " ... \n", "13226 POLYGON ((372750.000 6665500.000, 373000.000 6...\n", "13227 POLYGON ((373000.000 6665500.000, 373250.000 6...\n", "13228 POLYGON ((372250.000 6665250.000, 372500.000 6...\n", "13229 POLYGON ((372500.000 6665250.000, 372750.000 6...\n", "13230 POLYGON ((372750.000 6665250.000, 373000.000 6...\n", "Length: 13231, dtype: geometry" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.envelope" ] }, { "cell_type": "markdown", "id": "4edd8c53-c564-491b-b21a-bbc98b070a58", "metadata": {}, "source": [ "In order to get the bounding rectangle for the whole layer, we can first create an union of all geometries, and then create the bounding rectangle for those:" ] }, { "cell_type": "code", "execution_count": 6, "id": "dd2eb312", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.unary_union.envelope" ] }, { "cell_type": "markdown", "id": "acd68be9-1145-46ca-94cb-935eb00c59b0", "metadata": {}, "source": [ "Corner coordinates of a GeoDataFrame can be directly fetched via the `total_bounds`attribute, while the `bounds` attribute returns the bounding coordinates of each feature:" ] }, { "cell_type": "code", "execution_count": 7, "id": "f19b82ec-97bf-4264-950e-b91fa7d22ca7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 361500.00014042, 6665250.00004393, 403750.00013197,\n", " 6698000.00003802])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.total_bounds" ] }, { "cell_type": "code", "execution_count": 8, "id": "c73aee97-3941-4faa-954d-850f9c34ae83", "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", "
minxminymaxxmaxy
0381750.0001366.697750e+06382000.0001366.698000e+06
1382000.0001366.697750e+06382250.0001366.698000e+06
2382250.0001366.697750e+06382500.0001366.698000e+06
3382500.0001366.697750e+06382750.0001366.698000e+06
4381000.0001366.697500e+06381250.0001366.697750e+06
\n", "
" ], "text/plain": [ " minx miny maxx maxy\n", "0 381750.000136 6.697750e+06 382000.000136 6.698000e+06\n", "1 382000.000136 6.697750e+06 382250.000136 6.698000e+06\n", "2 382250.000136 6.697750e+06 382500.000136 6.698000e+06\n", "3 382500.000136 6.697750e+06 382750.000136 6.698000e+06\n", "4 381000.000136 6.697500e+06 381250.000136 6.697750e+06" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.bounds.head()" ] }, { "cell_type": "markdown", "id": "5c3bda05-d74d-4358-88f4-0beb7a71f4a7", "metadata": {}, "source": [ "### Convex hull\n", "\n", "A bit more detailed delineation of the data extent can be extracted using a convex hull which represents the smalles possible polygon that contains all points in an object. \n", "\n", "In order to create a covex hull for all grid squares, in stead of individual grid squares, we need to first create an union of all polygons: " ] }, { "cell_type": "code", "execution_count": 9, "id": "8a81a3c2-f4f2-423e-9737-d5030d1cbe96", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.unary_union.convex_hull" ] }, { "cell_type": "markdown", "id": "a6df7c4e-2df9-4d76-aa9f-0b6927697d3f", "metadata": {}, "source": [ "## Buffer" ] }, { "cell_type": "code", "execution_count": 10, "id": "045169e9-b23e-4165-b5be-f27dd167aadf", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 5 km buffer for the travel time matrix extent\n", "data.unary_union.buffer(5000)" ] }, { "cell_type": "markdown", "id": "05e5f7e9", "metadata": {}, "source": [ "\n", "\n", "## Dissolving and merging geometries\n", "\n", "Data aggregation refers to a process where we combine data into groups. When doing spatial data aggregation, we merge the geometries together into coarser units (based on some attribute), and can also calculate summary statistics for these combined geometries from the original, more detailed values. For example, suppose that we are interested in studying continents, but we only have country-level data like the country dataset. If we aggregate the data by continent, we would convert the country-level data into a continent-level dataset.\n", "\n", "In this section, we will aggregate our travel time data by car travel times (column `car_r_t`), i.e. the grid cells that have the same travel time to Railway Station will be merged together.\n", "\n", "- For doing the aggregation we will use a function called `dissolve()` that takes as input the column that will be used for conducting the aggregation:\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f1ddc1d7", "metadata": {}, "outputs": [], "source": [ "# Conduct the aggregation\n", "dissolved = intersection.dissolve(by=\"car_r_t\")\n", "\n", "# What did we get\n", "dissolved.head()" ] }, { "cell_type": "markdown", "id": "152bde89", "metadata": {}, "source": [ "- Let's compare the number of cells in the layers before and after the aggregation:" ] }, { "cell_type": "code", "execution_count": null, "id": "0976aa20", "metadata": {}, "outputs": [], "source": [ "print(\"Rows in original intersection GeoDataFrame:\", len(intersection))\n", "print(\"Rows in dissolved layer:\", len(dissolved))" ] }, { "cell_type": "markdown", "id": "194b621b", "metadata": {}, "source": [ "Indeed the number of rows in our data has decreased and the Polygons were merged together.\n", "\n", "What actually happened here? Let's take a closer look. \n", "\n", "- Let's see what columns we have now in our GeoDataFrame:" ] }, { "cell_type": "code", "execution_count": null, "id": "517500ed", "metadata": {}, "outputs": [], "source": [ "dissolved.columns" ] }, { "cell_type": "markdown", "id": "c68b18fd", "metadata": {}, "source": [ "As we can see, the column that we used for conducting the aggregation (`car_r_t`) can not be found from the columns list anymore. What happened to it?\n", "\n", "- Let's take a look at the indices of our GeoDataFrame:" ] }, { "cell_type": "code", "execution_count": null, "id": "32fa4a1c", "metadata": {}, "outputs": [], "source": [ "dissolved.index" ] }, { "cell_type": "markdown", "id": "97133b03", "metadata": {}, "source": [ "Aha! Well now we understand where our column went. It is now used as index in our `dissolved` GeoDataFrame. \n", "\n", "- Now, we can for example select only such geometries from the layer that are for example exactly 15 minutes away from the Helsinki Railway Station:" ] }, { "cell_type": "code", "execution_count": null, "id": "6aa6cf43", "metadata": {}, "outputs": [], "source": [ "# Select only geometries that are within 15 minutes away\n", "dissolved.loc[15]" ] }, { "cell_type": "code", "execution_count": null, "id": "c8b12863", "metadata": {}, "outputs": [], "source": [ "# See the data type\n", "type(dissolved.loc[15])" ] }, { "cell_type": "code", "execution_count": null, "id": "7c66bda8", "metadata": {}, "outputs": [], "source": [ "# See the data\n", "dissolved.loc[15].head()" ] }, { "cell_type": "markdown", "id": "68918dff", "metadata": {}, "source": [ "As we can see, as a result, we have now a Pandas `Series` object containing basically one row from our original aggregated GeoDataFrame.\n", "\n", "Let's also visualize those 15 minute grid cells.\n", "\n", "- First, we need to convert the selected row back to a GeoDataFrame:" ] }, { "cell_type": "code", "execution_count": null, "id": "d72a88c9", "metadata": {}, "outputs": [], "source": [ "# Create a GeoDataFrame\n", "selection = gpd.GeoDataFrame([dissolved.loc[15]], crs=dissolved.crs)" ] }, { "cell_type": "markdown", "id": "772b207f", "metadata": {}, "source": [ "- Plot the selection on top of the entire grid:" ] }, { "cell_type": "code", "execution_count": null, "id": "40045aa0", "metadata": {}, "outputs": [], "source": [ "# Plot all the grid cells, and the grid cells that are 15 minutes a way from the Railway Station\n", "ax = dissolved.plot(facecolor=\"gray\")\n", "selection.plot(ax=ax, facecolor=\"red\")" ] }, { "cell_type": "code", "execution_count": null, "id": "d07187be-cb7b-44c0-8a7f-519fb4abe0bb", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "a6694453-7486-4da7-bc72-a5a618b235aa", "metadata": {}, "source": [ "## Simplifying geometries" ] }, { "cell_type": "code", "execution_count": null, "id": "0df8ef9f-1c05-4618-be22-c006f5fc914b", "metadata": {}, "outputs": [], "source": [ "data.unary_union.simplify(tolerance=500)" ] }, { "cell_type": "markdown", "id": "43ae3595", "metadata": {}, "source": [ "## Simplifying geometries\n", "\n", "Sometimes it might be useful to be able to simplify geometries. This could be something to consider for example when you have very detailed spatial features that cover the whole world. If you make a map that covers the whole world, it is unnecessary to have really detailed geometries because it is simply impossible to see those small details from your map. Furthermore, it takes a long time to actually render a large quantity of features into a map. Here, we will see how it is possible to simplify geometric features in Python.\n", "\n", "As an example we will use data representing the Amazon river in South America, and simplify it's geometries.\n", "\n", "- Let's first read the data and see how the river looks like:" ] }, { "cell_type": "code", "execution_count": null, "id": "aa9964f4", "metadata": {}, "outputs": [], "source": [ "import geopandas as gpd\n", "\n", "# File path\n", "fp = \"data/Amazon_river.shp\"\n", "data = gpd.read_file(fp)\n", "\n", "# Print crs\n", "print(data.crs)\n", "\n", "# Plot the river\n", "data.plot()" ] }, { "cell_type": "markdown", "id": "e26624fc", "metadata": {}, "source": [ "The LineString that is presented here is quite detailed, so let's see how we can generalize them a bit. As we can see from the coordinate reference system, the data is projected in a metric system using [Mercator projection based on SIRGAS datum](http://spatialreference.org/ref/sr-org/7868/). \n", "\n", "- Generalization can be done easily by using a Shapely function called `.simplify()`. The `tolerance` parameter can be used to adjusts how much geometries should be generalized. **The tolerance value is tied to the coordinate system of the geometries**. Hence, the value we pass here is 20 000 **meters** (20 kilometers)." ] }, { "cell_type": "code", "execution_count": null, "id": "d102a662", "metadata": {}, "outputs": [], "source": [ "# Generalize geometry\n", "data[\"geom_gen\"] = data.simplify(tolerance=20000)\n", "\n", "# Set geometry to be our new simlified geometry\n", "data = data.set_geometry(\"geom_gen\")\n", "\n", "# Plot\n", "data.plot()" ] }, { "cell_type": "markdown", "id": "8a5a7dad", "metadata": {}, "source": [ "Nice! As a result, now we have simplified our LineString quite significantly as we can see from the map.\n", "\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }