{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 3.5.1 - Advanced - Geospatial Analysis\n", "\n", "## Outline\n", "\n", "### Prerequisites\n", "\n", "- Intermediate R skills\n", "- Theoretical understanding of multiple regression\n", "- Basic geometry\n", "\n", "### Outcomes\n", "\n", "After completing this notebook, you will be able to:\n", "\n", "- Manipulate geospatial objects in R using the `sf` package\n", "- Perform geospatial operations on real world data\n", "- Understand applications of geospatial analysis in the context of\n", " economic research\n", "\n", "### References\n", "\n", "- [Geocomputation with R](https://r.geocompx.org/intro)\n", "\n", "## Introduction\n", "\n", "In this notebook, we’ll introduce the basics of geospatial analysis with\n", "vector data. We’ll use the `sf` package and data from `spData` for our\n", "examples." ], "id": "e6633415-e2ea-432f-92fb-d011319fc404" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# load packages\n", "library(tidyverse)\n", "library(ggplot2)\n", "library(sandwich)\n", "library(lmtest)\n", "library(sf)\n", "library(spData)\n", "\n", "# load self-tests\n", "source('advanced_geospatial_tests.r')" ], "id": "09f193d9-c3c1-4f9d-ba13-b0fbdf46729f" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1: Geospatial Data\n", "\n", "Geospatial data is data coded to represent physical objects. For\n", "example, let’s say we’re interested in analyzing the accessibility of\n", "healthcare in BC. A dataset containing locations of emergency\n", "departments (coded in longitude and latitude coordinates) and population\n", "density per municipality would consitute a geospatial dataset. We could\n", "use the geospatial data to estimate the share of BC residents within\n", "10km of an emergency department. That analysis would involve calculating\n", "relationships between physical objects (i.e., location of emergency\n", "departments) and overlaying numeric data (i.e., population density per\n", "municipality) on the physical objects.\n", "\n", "The first step to learning how to conduct geospatial analyses like this\n", "in R is to understand how geospatial data is stored - that’s what we’ll\n", "cover in this section.\n", "\n", "### Vector vs Raster Data\n", "\n", "In R, geospatial data can be stored as either **vector data** or\n", "**raster data**.\n", "\n", "- **Vector data** is coded in a collection of mathematical objects,\n", " such as points, lines, and polygons. Geospatial objects with vector\n", " data have discrete boundaries and are associated with specific\n", " locations through a **coordinate reference systems (CRS)**.\n", "- **Raster data** is coded as 2-D cells with constant size, called\n", " pixels. These pixels are accompanied with information that links\n", " them with a specific location.\n", "\n", "Vector data is most widely used in the social sciences because\n", "applications in politics or economics typically require discrete\n", "boundaries of administrative regions (e.g., country or state borders).\n", "For this reason, we’ll focus on conducting geospatial analysis with\n", "vector data.\n", "\n", "Vector data codes geospatial objects with the following elements:\n", "\n", "- **Points**, coded as `c(x, y)`: the most basic element, used when\n", " the area of the objects are not meaningful (e.g., locations of\n", " emergency departments in BC).\n", "- **Linestrings**, coded as `rbind(c(x1, y1), c(x2, y2), c(x3, y3))`:\n", " a series of points, used when the length of an object is meaningful\n", " but the area is not (e.g., rivers, roads).\n", "- **Polygons**, coded as\n", " `list(rbind(c(x1, y1), c(x2, y2), c(x3, y3), c(x1, y1)))`: a series\n", " of *closed* points, used when the area of an object is meaningful\n", " (e.g., Canadian provinces, metropolitan areas, protected areas).\n", "\n", "> **Note**: to be closed objects, polygons need to start and end at the\n", "> same points; notice above that the polygon starts at `c(x1, y1)` and\n", "> also ends at `c(x1, y1)`.\n", "\n", "What exactly are these (x, y) coordinates? How does R (and the user)\n", "know which locations these coordinates represent?\n", "\n", "### Coordinate Reference Systems\n", "\n", "A coordinate reference system can be thought of as the base map in which\n", "your geospatial objects will be overlayed. There are two main types of\n", "CRS for vector data: **geographic CRS** and **projected CRS**.\n", "\n", "- **Geographic CRS** map locations with longitude and latitude\n", " coordinates. The x’s and y’s for the points, linestrings, and\n", " polygons introduced above will simply be longitude and latitude\n", " coordinates on the base map.\n", "\n", "> **Note**: Geographic coordinates are spherical! This means that you\n", "> cannot use the distance formula you learned in high-school to\n", "> calculate the distance between two points coded in\n", "> `c(longitude, latitude)` format. More on this later.\n", "\n", "- **Projected CRS** map locations with Cartesian coordinates on a flat\n", " surface. The x’s and y’s for the points, linestrings, and polygons\n", " introduced above will simply be x’s and y’s of a regular xy plane\n", " grid. There are different ways to project the earth on a flat\n", " surface, and that implies different ways to store objects and the\n", " relationships between them (e.g., conic, cylindrical, equidistant,\n", " equal-area…). More on this later.\n", "\n", "### The `sf` package\n", "\n", "The `sf` package is currently the most widely used package for\n", "manipulating geospatial data in vector format in R. The package supports\n", "all of the elements described above (i.e., points, linestrings,\n", "polygons), as well as any combination of those objects (i.e.,\n", "multipoints, multilinestrings, multipolygons, and geometry collections).\n", "We’ll introduce them as needed throughout this notebook.\n", "\n", "The beauty of the `sf` package is that it is compatible with\n", "`tidyverse`: geometric objects are stored in dataframes, and we can\n", "manipulate those objects with the typical `tidyverse` functions we use\n", "with non-spatial datasets.\n", "\n", "To illustrate, let’s create some geospatial data from scratch." ], "id": "2e92625c-2945-4580-8c8d-909b75740af0" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# creating a point\n", "a_point <- c(0,1)\n", "class(a_point)" ], "id": "f30b284a-2b8a-47e7-a49c-205378be1aa7" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let’s transform our point into a geospatial object using the\n", "`st_point()` function." ], "id": "a8a4d99e-83f6-43e8-b25e-64ef9b011265" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# transforming point into geospatial object\n", "geo_point <- st_point(a_point)\n", "geo_point\n", "class(geo_point)" ], "id": "e991fe3f-1b2b-46ad-ba4d-b01b53b7ffd8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the data is coded as `POINT (0 1)` and the type of the\n", "object is `sfg`, which stands for **simple feature geometry**. The\n", "functions to transform R numeric vectors into `sfg` objects are:\n", "\n", "- `st_point()`\n", "- `st_linestring()`\n", "- `st_polygon()`\n", "\n", "We can bind these `sfg` objects into a simple feature column with the\n", "function `st_sfc()`." ], "id": "6d9acc06-a5c1-495a-a84e-b26845539bcc" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# creating two more points\n", "another_point <- c(1,2)\n", "yet_another_point <- c(3,3)\n", "\n", "# transforming points into `sfg` objects\n", "geo_point_2 <- st_point(another_point)\n", "geo_point_3 <- st_point(yet_another_point)\n", "\n", "# combining `sfg` objects into a simple feature column \n", "sfc_obj <- st_sfc(geo_point, geo_point_2, geo_point_3)\n", "class(sfc_obj)" ], "id": "25e5ee43-ea92-4817-8ccc-22bd95bfd5d9" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the type of the object is `sfc`. R also recognizes that it is a\n", "`sfc_POINT`, because we only passed points to the simple features\n", "column. Simple feature columns support different types of simple feature\n", "objects in a same column (e.g., a column with a point, a linestring, and\n", "a geometry collection).\n", "\n", "Since we have coded our elements as geospatial data, we can now plot the\n", "points in space." ], "id": "ca0b7f67-a606-4cbc-82de-d73815258254" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot(sfc_obj)" ], "id": "42a1892b-464d-4b70-929d-5f2e38a063a6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let’s take a look at the ouput of `sfc_obj`." ], "id": "33016c10-e4c4-4a41-91a3-8bcdbb388a15" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sfc_obj" ], "id": "4a00f5a0-a5b1-4588-a462-5b6b92f106ad" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output of this object gives us the following information:\n", "\n", "- We only have points as geometric objects in the column\n", "- The dimension of our objects is 2-D (i.e., we only passed 2\n", " coordinates for each point `c(x, y)`)\n", "- The bounds of our plot are `[c(0,1), c(3,1), c(3,3), c(0,3)]`\n", "- We have not specified a CRS (i.e., we don’t have a base map)\n", "\n", "We can specify the CRS for our geometric data when creating the `sfc`\n", "object with the parameter `crs`. Here, we have chosen the `crs`\n", "“EPSG:4326”, which is a basic map of the world." ], "id": "e7da0d28-cc06-4694-9b2a-27918f0bc115" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sfc_obj_with_crs <- st_sfc(geo_point, geo_point_2, geo_point_3, crs = \"EPSG:4326\")\n", "sfc_obj_with_crs" ], "id": "6e619860-d651-4285-8277-25ecf6a59fba" }, { "cell_type": "markdown", "metadata": {}, "source": [ "See now that R knows that our data refers to the Geodetic CRS WGS 84. A\n", "simple search of the term will tell you that this is one of the most\n", "widely used geographic CRS’s, in which the `c(x, y)` represents latitude\n", "and longitude coordinates.\n", "\n", "Once we specify the CRS, R knows which locations on Earth our geometric\n", "objects refer to. This allows us to overlay objects we create on\n", "existing geospatial objects that share the same CRS.\n", "\n", "Let’s see this in practice. Let’s load the `world` dataset from the\n", "`spData` package. We’ll use this dataset, which contains country\n", "borders, to overlay two lines connecting Salvador, Brazil, to (1)\n", "Luanda, Angola, and to (2) Maputo, Mozambique." ], "id": "0438fe7e-dba8-466e-8335-01534ea037a4" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# check the data and plot the `world` geometry\n", "head(world)\n", "plot(world[\"name_long\"]) # specify that we want to plot the attribute `name_long`" ], "id": "5baf2ffa-03a2-47ba-9bf2-0100a425e7c8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see from the output above that the country boundaries are stored\n", "as multipolygons, and that the CRS for the `world` geometry is WGS 84.\n", "\n", "As shown earlier, the argument used to transform numeric vectors to this\n", "CRS is `crs = \"EPSG:4326\"`. Let’s use that information to draw the lines\n", "from Salvador to Luanda and Maputo and overlay them on the plot. We’ll\n", "need the coordinates of these cities, which we can find by searching for\n", "them online.\n", "\n", "> **Note**: for this CRS, we need to specify the coordinates as\n", "> `c(longitude, latitude)`. The longitude and latitude of Salvador are\n", "> approximately `c(-38.5, -13)` and those of Luanda and Maputo are\n", "> `c(13.2, -8.8)` and `c(32.6, -26)` respectively." ], "id": "5b6da087-d603-4227-b029-6129ffc9c266" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create lines with `st_linestring()`\n", "line_luanda <- st_linestring(rbind(c(-38.5, -13), c(13.2, -8.8)))\n", "line_maputo <- st_linestring(rbind(c(-38.5, -13), c(32.6, -26)))\n", "\n", "# create simple feature column\n", "lines_sfc <- st_sfc(line_luanda, line_maputo, crs = \"EPSG:4326\")\n", "\n", "plot(world[\"name_long\"], reset = FALSE)\n", "plot(lines_sfc, add = TRUE, col =\"red\") # overlay the lines" ], "id": "2f58b70d-3f14-413e-93fb-2a6746f11dd1" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `plot()` function is the easiest way to plot geospatial objects in\n", "R. The default of this function plots one map for each attribute (try\n", "running `plot(world)` to see for yourself), so we specify `name_long` to\n", "plot a single map for the output. The parameters `reset = FALSE` and\n", "`add = TRUE` are needed to overlay plots.\n", "\n", "### Geospatial Operations\n", "\n", "We now turn our attention to objects that contain both simple features\n", "and other data types as columns. The `world` dataset, which we used in\n", "the previous section, has those features. Let’s examine this dataset\n", "further." ], "id": "e1dc90f9-a454-4888-a275-176787202b97" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "str(world)\n", "class(world)" ], "id": "8606232b-6be7-4748-8851-e310e60c0227" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that `world` contains the boundaries of countries as\n", "multipolygons, but also contains attributes of those countries as other\n", "types of data. Recognizing this, R categorizes this dataset as both `sf`\n", "and `data.frame`.\n", "\n", "We can create `sf` objects with categorical and numeric attributes using\n", "the function `st_sf()`. Let’s show this by creating a custom dataset for\n", "three cities in Portugal." ], "id": "2bc7e011-309a-46e0-966f-8767992348c1" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create a data.frame object containing the cities and their populations\n", "city_attr <- data.frame( \n", " name = c(\"Lisboa\", \"Coimbra\", \"Porto\"),\n", " population = c(545796, 106655, 231800) \n", ")\n", "\n", "# specify points as c(lon, lat) coordinates of cities\n", "l_coord <- st_point(c(-9.1, 38.7)) \n", "c_coord <- st_point(c(-8.4, 40.2))\n", "p_coord <- st_point(c(-8.6, 41.1))\n", "\n", "# create a column with the specified CRS\n", "coordinates <- st_sfc(l_coord, c_coord, p_coord, crs = \"EPSG:4326\") \n", "\n", "# create `sf` object with the data.frame and the sfc\n", "city_data <- st_sf(city_attr, geometry = coordinates) \n", "city_data" ], "id": "70daf619-8af0-4fc2-b192-dc32d14f4fcb" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class(city_data)" ], "id": "c7cd86a4-d866-4fdf-affb-4889eebdd14b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The image below illustrates the process of creating `sf` objects from\n", "`sfc` and `data.frame` objects, which we did above.\n", "\n", "
\n", "\n", "
Lovelace, Robin and Nowosad, Jakub and\n", "Muenchow, Jannes, (2019). Geocomputation with R.
\n", "
\n", "\n", "Now let’s suppose we want to merge `city_data` with data from the major\n", "rivers in Portugal that pass through those cities.\n", "\n", "> **Note**: These are not the actual coordinates of the rivers. The\n", "> points are simulated and made so that they pass through the major\n", "> cities." ], "attachments": { "media/sf_diagram.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8cAAACtCAAAAAC441InAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAAmJLR0QA/4ePzL8AAAAHdElN\nRQfmARkNNSsa7cs/AAAn/ElEQVR42u2deUATRxfAX7jFE1RExIqgcnriTbzxJtaKWK0l2kMQr1Ss\nYr1I1bZ41JqKBx5fJdTaindQrFJBBU8QVBCwAqJcIjdyCCHz/UG4NMcmu8kmdH7/QHZn376Z5O3M\nvjfzhoEAg8FoOTp0K4DBYEiD7RiD0X6wHWMw2o8e3QpgWgWxgaqVb7md7hpqNrg/xlDBiz9VKj7l\nPN0V1HBwf4yhhHbHVSk98BDd9dNwcH+MwWg/2I4xGO0H2zEGo/1gO8ZgtB9sxxiM9oPtGIPRfrAd\nYzDaD7ZjDEb7wXaMwWg/2I4xGO0H2zEGo/1gO8ZgtB9sxxiM9oPtGIPRfuiz4+Ldb+iuPEZlhDPN\nZJ2u5vzWezLdOrYm6LJj9LvdZrtjOFlna2W6Q/PvtuD90z+Kvghzp1vH1gRNdvx82hKfN1tWD4+l\nu/4YFdGx2f+Z694/e64LOC6lW8XWBC12XMV1qo3ntuMk9x7FKae7BTDUUxGdBQAAtQ9iygFyZmZn\nlgGU3kqoFZ/OzCnJLAB4Dql5ANmRzwEAirMhNRFA+CADAADSbpXSXQttgg47jhxy+PB1OwDocers\nRTs+3U2AoZr/Lak5dQEAYkdUvOqbAI/e5p5/BscWGJ4cVQYAAC/Pv006/yCk34Ifh3igpTvbrlwC\n5Wu6BXmvHbw26vOVfcMByr65fsPqLN310CaQusnx1PEqbfxU4W/glqF2HTBUE9q56f8bH9Ui5NEF\noTnTEBrujdDI2Qgh0wD0lnFSXKTdJoTQ+PZPsqMfw110GXJRKcyvRZt0jyI03g2hL28jNK9jWaPI\nfY50V1DDUXeePdHRtTZ3hjd9NuYuWOa4doMB3Y8zDHVsn6EH0AsAtjKg8F1DVCJ0IGQxWoQoutg4\ngYXwwlDIgTfmHWCQHgyoGwnQ+wGU/tb2FBS1SRxFd1W0BjXbcfzS1K0rdFscso0I+fbkwYl0NwSG\nMh6OEP/jmLytb3uR+MPEq7HDGKIWBfUAQG/m/8rbgvi4AQCAbjUkwR6cyFUh1Pp+XMoZZpbI0X3v\nKIOdOm0KGweTWw1v88T/BH/97XyjhqM+FzZMZnxYuITZ0dfp/YNClAQAUEt3TbQHddqxwOnSJYGl\nhBMmvMiHtjyRwgIxGon930IAJIIanyltQIgAGNUAtw7NAyRqCCqjxujyvlh3EMJ7EwlsdX4AgOdn\n6K6J9qA+O06b5r4wcaqUk2Pi/TeNS6S7MTCUsO7V0jfpEUV/VdWee8jPfnEfuj3OjXkHRxN+MUx5\nDBMWAdRUlQNAcTkAvBMG3j4N8f8WQRVAFZQDVFVCN3bo3N83bcQzRQijy1XPfWp3zzMPWyD9pUd3\nJDvat2AM9ndpKU/DmuZ69Le5vL98qAOzj/n1h3Od/+4yumf4jQVD8qKyVtf9M8bmbJeZr47UVun2\nPhdvXG1l7Pjout730Y9cf0c1FtWn2ubZPrrRXtTfrTbmVr/dbRpF3o9dRncNNRuGeqZGRi0r3uHJ\nkFNIsKKOhx/B2snppQXkhUgn8BAerMlELePqPPakoYlseWYMrOSvP2Nl0t0iGIz2oQY7FvGdnsTw\nOxMoacx9UOjAraG7TTAYbUP1dpzgsmxt7EiChQfEHNw/7Da9TYLBaB2qtuOK9UO7JPvpEi7PYKeO\nHMNW6bsWRgvJW3s2m24dNBkVT5sRLNcXTFfsEtOghT62Py2R+zaN+S+hm3C4zHLEyBHOxiq+Ufm5\nEzWfu3dSU7X276JAyO9MFdtx+ooInx/bKnzZ2IQDvicOOqhSM4yW0fUapEfHnd9YZ+vMdLFX1Siy\n9sofF4zmGqxfMePzGYbqqFaJ/kbSMr6uApXace2Bjc4JSlmjPoe1YpAyTwBMK8bamg0V8XFxP2Z2\n6M90GdWF8hvE8f8snXxorjHwIvls3VkeM4i/DiqN2WLSIpYAqNKOb/rk7VF6dGx9WbA8LFDBETmm\n1dOWyQTIiYuJ/rWqO9PFeRh1nWbyXyfSR21Y2AUAQNfVtSqM727m7sGku8JEUZUdF313dOEeMo9M\n1sRtrOkHetLRJhjNxsKCBcLUmOjDyXoDXJydHclLzAkNjXFYuKh305E2Hh5Fp/n77D08beiuLiFU\n86aB+LZ3o/nkRj5tA2IL7HfU0dEoGI1Hz9GLn1R82S3d18mCxY2oJCGqjM/qtdv5YRK3d8vjpl7R\nGeyTfYbytGEpnkrs+LGLz/IH5NeAD4o5sMv5rtqbBKMtdHTlCgrSAqwjZnZwZB9OkrRgDsmx8LoI\ntgXHJPwlb7Cks738UhNdA3qy+BV0V1UeKrDjSu6wzk+5VKx40GEnDnBhF6q9UTBahDWbF10U5QU/\nOXVirhe833c+sxfIuDaO051V/Ntrvqt0P45jQFaExYoebIGQ7orKhHo7Ftgf/UPQiyJh5vzrsU58\nnOYaI5O2TA4/I/uEa9ynZhbzeNHvms7cy5rNein5oqfcPsPjtr8WeMjpcnSZQa+PFLv34kTTXU0Z\nUG3HWXPdZ6dQuWhp3KNvvCckq7NJMNqJBYt7rSyRa3R4bPuhHH5S/cE7uqK/+0qYsp/NYzqGfp4W\n7dWBiOg2HoI8/7ixDtzndFdSGtTacS3PPv8hrx2lMvX9Eo0Grq9WZ6NgtJUP3F83aqG25kfbf1qU\nKuWzrH52jk/iWhGXbOoV/WLRn32H8vLprqNEKF1/fMsnx3+lKjxngmWG+6eSF4NRGac/n69K8Slv\nFVx/LEq5d/duEsPxiQgAQAct2NtVfObd1dAzhm7sSUpNbUgKCX4zwdOdshlKP1yOIS1D//JkSvvj\nIs6EIakclTjAWYkzZ7KyVCEZQw1WKjVjsJut4AU6Dl8EPSq5Nq6+mxKh09Y8EQCIojk95xUfz5Pl\n2ZKFY0B2lPWKrvM0zutFWX+MQr41UWX22vilqd+vUMNEOUxrYteWxhcyPYdjbUJDXoxizyf0SiyD\n6mshF9q7e1Ix10vj+uMnY3yWPVFlEurBd3ZuGXZfhTfAtEJuN7mu6x4Pd/rLK4ugZ0sWRqxTeXvS\nx9pz/6W7ek1QY8eV3KEmSZSEjGVo6pXiNMq7TB2NgmktxCDQMTQA0DM1hY4zT8T4dadGrgn7WubS\nC/2G8l7TXUMxlIyrBStrf2KrQ9vIZSU71HIjTKsgu4+5dV8ri9f3LxvN9XRR8JVYGJ02ygEAcjoI\nO0kukRTCz5/gOYdEfIbYuLr0Rtkk6Q8gqsbV2fPmfJyiHuua8NDba2KKWm6FaQWYV2Vc+zx3zZbq\n43lBTAXNuObTTun9swF4KUaBCZKLOAZkRVmvNJsnUO3GFy+/cgyeIKcMaTsW8uwz7vLaq7QiTbTh\nJuoP5r4jLwjzX0A3idtnfNwP+QIPfZnl8nI/PHYqa9D26z1gXb+JehsD06RcqMMMev0XzDX3jiY8\nrs113qNghqKtNjan/5BThqwdxwz5fts9Z5JCFKHP3/xDTlfVeEOMtpLFc3YK9c6K9pLXy9R5Spjc\nccMYGOPgdvp0AMa6hdJ7XCPWqbxf0sf2Xv+MoFoP/Xq6HClSoB43jKHjEDllyNlxMWfcoBTVhIyl\n45EyY8Y8TfEvYDSUEv7kj36ZnJrkZy6rVPXRQ+ee1XwacfH3ZgcFv554UBEcnxd8CsB3IQBAP6Ng\nWVJM2NcyV1+1ddyRB0QQojsruo7nl8stmMs7cqkwLTg3IfimvKJkbBDx7cKv8M1ISFCOTrz7L+x4\neGUyRhrVgnndVlvfzAjoJ6cgy3rpi/O1E8F2QNOxH2JX9t1uMMjYaJAjxNyrj6WO3S9HUE/Ow0SP\nIEvmYfnWCQCoRhT9VZfpobJTtRdOWvjFb2mmg0RdBsldd0QiH8iz5THrvlNLNrIPGHJ7/xb+oWG0\n3Buj4dRFnjgrmn3BlcBPuzTiAHCOtLVpYcenlzGGT9Uf2BENBEhs1xEAAKwfyRfm6Lgl8sQ639nD\nDWR2jqViLQGuXTO2rRZJLxxd1UFv6xsTE13zgXJvrrQdV+0IcE3qrezVZNHjeKwfvWwb6ZA+pvUR\ntjtG13dJHyJF29tMCpzl/d7BwatKOQ27wiWLt0ExIeTG0pnUx3LniegqmYGoumb/lcd1+Ft6Ejqn\n7AmH+hOtNVKOsN4WwUpeShVhvbvTrQJGI3m5dzA4BOQSKJniAKxidAUSmh0rnwX2iQi5MRFCi6zq\nj10iYifFwa6MXn6pckrlNFieEfTxXzNaVtHzJnobRQi12yRToN5VhJBy78c57NmsZLonZMxM8vKa\nlEqzEhgNhPDLal2v+B3X3N4/euHs29Fif7KZeBRcLNdOCL+S12MA5pyUf7kmssq8/TiF/cNmYlVW\nxo6FPLuk2zz6x7RtuE90cDAZIwlH7vMo543ypmgUnDBY99vtal0oaXaQB5/crnwAIgQAo0urAADg\njbXMu4mivbt9CqeJTjbRZ3RdGZcbYCun2KU0s2Pe/0C9JnJQwo4fjuZuu68ZPqa+Vw8d7B9BtxYY\nTUSHycs6ZTRf9hSNX94CY7JRH8ahP141HjudDAzz4fAqRwQwpXN9nseYz2TcKYlrMz59X/4pluzJ\nJo20XxSZt1teQBgAYDsCxiworHwlv6jC86tL/PfP2ddNqZZVCcXcQPdA9ce+MNpBycWQ65YLvpDS\n9ZUfMdAv92kLR+MWNq1CPFLZMXOR2dk0sJzSEwKf/QoAFcPvS0sd8Ors8QRnz/nELeJt1JSmBUUy\n51cn3DIq7+KZJXgD9m4ytrXSvzwZFPZznTLr87fijgcZb/1eWRKORomuEBcR3b/T3joqdcK0FpIO\nXE6qfLV3CAGvV3WJGGGLwyKPFwih70IlX1QU7Mqw8numvILbRyt/bQNK+LmeT1nkkziFokdlvUS+\nhNRlP9Yxum8mPFBwecjdPDyWSqUwrQTRspmOxv2PWi3utbfHKNlerz/Gi3nS4jDj+N4C4PebK+GK\nasE8c1/rm+kBfemuJoBC/XGlv+GEZIUfF7/LPi38sDA3BCEUvlmBe2R76q0qI/9kw7Qy6sTDUYM2\nugAMRq/j7xSXURtR+FCC5FteHdp4XKwlqSAd/XHkkMOHr9sp+py47Sv7fItUPbtPA0Ba+OcAMO2m\n3DmlTVjwz16046vqUYfRVnTE3qSaqjrQgSkjl5mzIxR1COlNMv1gr4mk9Zb1ni0Vbx+ugJZEC+b6\nnfh6l0KxpsAK3eIf4ma9XW3XOGMm888uvY5WL50GUHoAZQxdogPJJ12mVoaf3n/02tCt+ke/67Xa\nbdK6+hkuszcokvebNWnnktB9VhS3zisKJnEbWFCslHzKjQh6TlVFwVsKhFhSYCRjHjSEJQ1sjo2C\n0guh07t/8sUgUjJfnvwt1dlvgYb5VgkOUII6DL6nWHcfvBiJpiK0vlPToeo9MJwb6qLzD3ptHYdK\nBs1EGWvgOMr0BO+rPDiE0MjZCFUwriGEELoP6YrdMGWisb8S4yZZmJBvX+hPrUpyyT8yw7DT56cr\n1HzbFiyi4peZQVKJnIv+bm3Fwz0DwwDxEDhrrws4BOQoK7QwyIVh7/8vVQ1F2bia2DMvfmnqVkWT\nVT78953h6paHDFf7j/KHKT13TNzYbwh0DJj216erfwb4aFrIWpvJu+7Ud9upqP5J1x1iFJu+bRsR\n8u3JA5Oo+Ak1coh01uzjZylVSA6Z587FdP34fMnZxXVTP3EzVeetW+C+m6yErDEkLq56eO/u3Vft\nh46axwYA0EETDzfswNuDw3l66sgGpdJmVl8LOd91boAmbopMxI5Lt+yffsZSUcnu+wfv+8AKdNsA\ndJh6C65PBIDJ+pGf6gEA6IIegLl4T7sUqO8FTUDRNcYMNos79bOfuyp4mSy6WpGVoEZbSheE3v7o\n483j9QDmV0cL/L4c6TG3h/pu35x2VvTcFwByYqLjYt9Zu6xxHqEPsD4HwNAsqMViBAcuN46/iePK\n/liRzJB1d0JO6sw6M11jXolbQMDPJXC6dEmgsBnDmMsVrqskvl/2MIL8fADQ6V4q4awIasV/Oyl8\nTxNeZLwtT6Twda2AJK6DzWHmzQxe/Xo9I1deVhQz0NKR+x9KZ1YevYNl1mPJU9fQwjQ+h6kPAEwd\nfb11zz5YU+TMyzpl8oUCXq+k9Zbj0wOz+Zrj2WqJXLXSll/35RopITltcvL3Owd9CR+2VEF/sH+E\nGABvP1iVhQAsxZNdi8BeibuOeXhg0+mDTipvN42i7k7o2SyHefNbBhN0mcyApNDQ763dPBRNFal9\n1KXExcXE6/Rj7nZ2aF5Zl1NjDktcuWDAYu27EDq9+yeLB8sVnvnn/545r/+MyoEe1cjpj2t3OFU/\nClDGjOFItfGO8RlgWF7RvH+sBaj9xxdWZN4ASBZ9BXVQByACAKhDAIZ5IBrRJgMAALJ1lLFj0Oek\ndBvMocJfqiVUR3B6jI/7NiuJKyEm6MhNSlsVN7Y3J0LTNjKhkFwBl9XFaX26x42ypCC2Y4tn1sQ/\noqQuQOrIFmSuiRviyM2QJb3oMLM3f8HzWI4mm7Ecf3WkvXmwSEk32tKlVc+cMlA4zPqhaa5Hp34H\n7i7YiRDys4l+Mvsqqj0BXpWVy+A3YWKHPtnIT8/rBFq2CiGE0NYvlfbgXezV4zR5PyBCyOQMaRE8\nlfqrKy56djByC3otu9TLIDf9zp4Xq1WpSQsWLSItIoOQv7omdq+nA7RzWXUqT9kbPfXvo+Oyt0Dy\nycpTbvo9Vt1STSshCv3VstZJ5K07sfCXzso+IFKrntdOMwG4KpzUlPzHxHta9ujeAADp9zqN7giv\n0gG6dnwOMCClBkwG1J7uPRLypt4xBqhzOad87v/KnT9NCaRgL3XTo3PIivj16GPyekim8FLoNWNX\nt0+IJB0uCgu7jCZ6zFbPatPFcJysiBe9M6xkl2hwaDkzB5Hc9yuO/2fpZI+5769FqIvkn9Od5aFK\nzxZl+ztJt2PR7749D46kRFuBeGTddbTJ0p/kF4+4vwFgj+liMrd8vPTR2g2kN6rRYDvOPB8W1Wm6\nx1Tidaz8J/TiO6abAktzlEbVdlz+KCb63puOw1ycXaiJB9RF8s/psdgTm71mxvH/KpFg3NRCmR1L\nfdYk+DzZ/C1F+xveFL+c9RtVQ+Q1zbUyxPP37gtI3XJATMiacwdHU6O/xpEuCL3da5bfeIV6CmMW\n690tQcCakR7uikcfNAZpDi1y6Lq6VoXxp5vPWVQ/k/PFX8fSRn23sAvdtSWO5DH3Wz9dt5eUvw1U\n/AwO52oIFMwUUXDzQi8dzzfkRGjk+3Givz04+N1S0m9Rd8uvLzj4P6Vaqxao6P0456K/Wyfo7hZw\nq0olahcEuTAc/NPEf1TaQmJUPJ9LsFxfMF3BBwIB2nh7AxCZ+vsR9CRQSg6mQQt9bH9a0qpCLnV3\nQs/kDp63QF5KGOnoMJkBSaFh2haNqn0cHRf3tN1AZ/ZY1b0XdPbySjkRso3RY8FBwnkqSRJnRVqE\nEEDyuDp9RYTPj20VlEYEhiqEymBswgHfEwcc1XtT1VF9LexCwci1HqQXXjg6ctMFoft6TnPT0NlJ\nLWlwaHmRdmjJx27b1jtCptp2SGFRMd/OCSTZce2Bjc4JDuqqiGrR57BWDFbNM0ndVFwPPS+ctPUT\niqKY1hzOq3CBe4cZHlPo2UqAGA0OLVc/ihxa8mGo06kyYAB5GQAgwY5v+uTtaUVjUevLguVh+2bQ\nrQZJCi6HXm3rGkhmp90P6enlVRQWtkB90SgFeXYzJjqFaodWa+U9Oy767ujCPVrkpSMAa+K2WdP3\nf0S3GsqTeT4symTaaQUiTIQxZbOrIkJXeTPdPjUnL41ipnZ3ZrsMVWou4X+PFm8CiG97N5pPpRnX\nxmx6AXA4iFBhqcVuoL9J6NA2ILbAYYdqtnWrOcr5vkQlksUk7WD23usQmcdnKWrGZ7m+T5s+/ZsT\nJTnRdxsWv+Caww5L5o7nBAVTQvbBrw7JWaMQmyPwY2IzJkZzO37s4rP8wShKxcfv+6EU4OwZQoWl\nFVMs7Z4EBt0+sMv5LrGyxLa+bMB9oG8guc2Yyyqln0vi2jvxmTczeEo4Xn5O23A/sPHTuasWA9eV\nSimqx+RlRTGP9XXkPiUmWxL5Cq0yK5z15WxOruwySs8k/E/SGIaq8Ddwy6QiJtYisV4YJCBUJyPW\neb8pa6iUYoqn3ZNArqeOZwGRgiPtdjS1grz48QMQIjkpSOTFj6PafCqQFFIX3lrVQ8fZX95uQQgh\nVBL24bF3JhFI2JgN+Jo3Qug1SyhTTKK/M1ivUjYuze3KaUoYIzd+vIuJ5DQcsfnVGDGND3qB/dE/\nBFS8RbZMrKcHAKAj3U/x5tOm/khyMSXS7knAnH891pFPoFOvSdnca9jBAmJSE0EXyL62Vp2e3dnr\nRsverFrgbTE+bu2rWC6B3YLqFkkYEecU64Fuw5cr9N0OAGaDZM+JdeTGpq+KG2vlLXsvFWkUHxzR\ncwvRxc6JekC64TDNEPu5sr656PMDWW/o39EWGcsKWyTWAwCAqrAbgaJbZzzyQjrsNge4d+5fhw1t\n4NlvvbImDJmasav7rtw/2/cMXDkm7EYgPD3TeeQvhVtGAkSEv+vkYurYTqm0exIY92iP9/8OElgL\nWQNxj1cM/Upm3peyPSalll9ui4BvGN+ZNWiqpGJ1UH78mOlnHg3pYkquCc4LJ+2a1UnulfVafHWh\n5MU3TYtCAit0i384GQb7znuMjrr+VmedGRzs3AUAYM7YZXKCN705nDfhoe4dZrBmKByr06uGrJ3b\n+ixky9wKqf5b33rr3er2W1HwyyxbjjZEsLUAhBCq2dtuTCLprv25rQj532mRWA+hK5CAzvSxQqV7\nYCb/fNf5CIVvQuX9JiNk8wJF/oSuQAKq3AsTf5m1/0wfK4TCjZwDIkf0EKK/LWurB5nvzlE67Z4k\nDafq+8mb0SfOlMow0J8WXCF1XP3FcVSwAKHdIESNmkouKXdcLf4eDMDcLwWhN8FuBiYeweWE6lOv\nRTXsbXZMnN4wFiIQOsgRiexGIjRiI0IIobp2/yMktjDYo52xW3CpQm3bkGvCgOG897X0cXX9t16/\nKen8MyhLZ72UgnhcrRB6AHDTJ3/fIvIhuif5rz5aKslNNCc8Ajp86vuJJ5y+BaIlf6WCy2+p3dKe\n9BpfDQAAbThc02++AQiPAJjW094Plnz9yurgED29r9YsN4pXNu2eBGyunPQ96yg7YWmh+PFWA9eu\nLRM++VjyDKKHptC5Yd/YBk2lCCz2Blk0OHtqIG/vju6GL7vPvjyOaBfVXIvGY03pDTM3v2DAzjsA\nKfW7jOl0i/mCiFhTNvtt+LlV3mNqehF3Fzesxa+Bh0/WdLN5K3l0UtL0rcOp3DnQYxvpLxUDoNy+\nqVIYZzj4kNkgSWf0oCmXXnrO46ioEYeMOk39xLdsWkMBq4ZioKsHYA4VoPsSoDfSIZF2T5XM/3nm\n8wblGzSlV4sG3O8N/kec3jCqe1sA1o+QWyrO32tCexM2/9bDBwDABnLL2jBi9ABgbMKBlcfI57Qy\niV7sc/KMvOhzuujT+l/VGc7e0AvSd4/84uO/p/7jY0Aq7d77EMk15pwJAAAMfTRpwVzL/lIm9K43\n3OR0xLOlptKaRXbw/Iag/q9BjfmiL2yh4HLojLaubsRmbjXXooExl792XfmLLgDAsxIQN53YcSUi\n1oRFYWHhoom/KjbL6/tb4nrUDvFcsA6k6d/sW39mo4h8jEx0AOpzWg0hndMq3fomP84XQLZT2ALC\nAABiqsqO3jH8BKQWn7n91PaRPLJp95pDPNcYw0Bn2L6CcLb0NeRpq1MnsNNbakoCfZ0uS2/lBtgC\ndGELXgcZrTBj8UvkX9agRfMmTJucvG5fMAAA9MhKBIDst10NxLKKCDThGz7L3NfojyIBW4nJmobQ\n57vnsRzpby6N3zoA9Ih8BwCppBoOI0Y8HOxx+sx5O2KzNaTy4C7Dc+37ifXqAAHUNculZ9v7m8dQ\nscG07DiMOJZTbQh5IAJR/dk6qC9XBwjCao9tmgcAJNPuNRE1cG9QJBEZBuC89/U9L5k/41/B8rxR\nBgihDho1VRZd3faLr7/mNWY37+TBz//LYm0XJi9HzpX1Wujr5EFTi4vTGwqhDmYYLM6GzH3tDFzq\nm1D4Wt6YK4PH7LbZ4kwun6X44gmhAViuS/6XK9NdLf7WQVgHMCdraTXciiTRcpgmGj1eFf4Gbi/I\nuMz+HJlbMON8y8R6lSsgqLZgqMFTYQh4v80bBffQtY5g0+U0yuscW7V+HSo0HrApKwZG5yNUMNTg\nqSiprXNe5Qr4WbiMYe4wnPMUUZB2DyHK54HYBddeZNaUzIIL1ahRU4mQnAeSIuvaei3QiK5r7jQe\nq09vWP0j+JWjI/o6lgNzETo1ACGE0E1bmTkc1DEPRPytv+zZ+ZlQ9CUYmi2QNjcF+6sVonl+LpI5\nrXIKXxSP7tMysV5qHoDVSxGAw1OAj6rzgTEWXt9EE7rCu6Si7H6jAJ48dTW+DwBjGbdEAMPihdDB\n/BnA4ISb3QreFt5NoCLtniK5xvKaLxiQmp/r0dsXHabr3akBMLO/2aCpROTl5yrTkz56Twr9K8WB\n5SZ1H5J6LeD11dFNr5r16Q2T8wEMRsGL26auegBo4s9DAGCpm5vUBroddvZfB495yq9Xze/SzNMn\nNT9X/bdedR8A+nWHh08chkmTJz/PHqYZLfLsoZA1FpTktLqTL5Y+S1kJyYvvAQDs8tWlIO2e0rnG\niOTZa6apJEjm2UsKE9zuNctjtBx3eHHDbDd7iRPActYf14H7R45Ivlh4NzQ0n9KsXWrJl4lpRotQ\nJYPt9t0YKtYtPhbPz9NR2o4jUy5P0Su9aKFLPu1exbbd05MpyBMkX1MV4Ojol3k+bILJNNmZMcui\nxP8YS7Rji43bv4u/KtFxXp9F008D1y1iFOH9gfYNB9MgZVPPU0jFake93rNj6z+QSrt3saf1ZWWv\nJZJnr4WmH0JJnj2FZnlJpKxQUhMWBnu0VXjmFgHUloceI+bD/NW1BzY6t56cVuRyjWlQ/ur6vD6s\n2dRtn/0qXPC3ivL64HG1uvnwvUuf87jtYE4F3YpRQi1vQEUCrzWk54K2LH7+XxZbLORHowiRzmP2\n+sH6cp4yESaM5iHJf2J9+cy5AeF0a0YBNwdt3RPVSlIGAoARKyg7ynl3z6FckpMnkrhDbX51vvlC\nvMcqRvuR7AdlJXuwWC/p1o0kRd4TnFO9WleKNl0m79Vjt1N2juujlUyQIope39cp1O1pGo/Zuprm\nv42UeIYqc1qpB+pzjWkKjtynaV7RY62V2Av1XQSnx/joFa+SuGQnuWI0C6lxyUExxHNaaSAqyDWm\nQVhzol9883R6d7aghvhFlQK2GSv9++xojhZv74SRjPT5BTrsxAEu7EK6FVSKSu6wzsncVp045iPO\ntbyfiz26zeOXEyleyJ9nNr94X77ASw3bLWLUjqx5Qub867FORHJaaRqU5RrTbDqzBUV8o5VmrMP5\nsgu+PMzq7mt0skjAJrJVMkYLkT3fb9yjb7wnJNOto4JkzXWfneJOtxbqwZjFfy2w3mLB5GVLK5LO\nY1r9iCNMrRw583b1/RKNBq6vpltLBajl2ec/5FG6gYpmY+TKy45y3m3pyJWQrDKJ62hzmHkzA0eY\nWjly09HYXDlzwukK3WoS5tbg77dHkU5tomXoMnmvEj1C7VtGo+qi1/d1CvVITgrAEaZWD4G0UqzE\nmW6sLLoVJUQRZ8KQVA4dqbJox5GblOYVPba3OBpVHcGxHB+9IiuJa0e3Zhg1QGS41ZG3eKnT9ytU\nvvcsWVDIt6ZXJ1Iq8s0LshKK1Fd/aw4n89y5/V0/nl1yNrxu6g43dW01+iFvX5CVoB0dh8bAIOaO\nFh1da3NwBN3KyuaJTzyZNAiSMC0mL6M/NeskCPPmwrl/2rjNnm5MXpTSLA6mQAheJ6EABO0YINfv\nxNe7NHKf3Hoqd/40ZZ8VxUJfUTCjzcBC7W1RbqSv9nu2oIBszkYAAEvsmiMOYTsGiFxWsoNNuLRM\nnibOg6KzC9vAlU7Esu3IRbBS+CNFymEwWocCPqEJ8d5eE4luxCWbzfOL4cSSqwALfSiRlz1vzsfJ\n2Iwx/1kU6I8B4Pnym37fUTCdID11OpRd+dgQotsPJC9NuH/zgIP9VdVCGIzmo5gdA4SubL9/Ct1K\ntyTGJ8t/5X8y1oTBiFH09++RMmPGPNq3CWpGMWfcoJT/ZsgYg2lAYQPoxLv/wo6nKSuTEd8u/Aqf\nupxVGIxWoui4GgBAuH9Ln0PDFL+Oep4tj1lHxfs6BqPdKDMg1eMkO47mlNGtOlRxBxgmcbEZYzDK\n9McAAJdWVgfQHOi5tPLdTzjWhMGA8vuYz0zy8ppE556XOezZLBwyxmAAQHk7hjbcJzqDue9oUlvI\ns0u6zdPgaaIYjDpRdlwNAIBC1nY84EqH1g+X/svV/PVXGIy6IBN4ZbBTpk+dl09CgnKUcIZbpXCw\nGWMwDZCbQGHCu5duyxOpV+VQ28uXT+GsjxhME2TG1QAAINy/ud+hoepT+PmyaBwyxmBaQnpCox4n\nxWEUh1AOZQqo4joJH+KQMQbTEtL9MQCAYFWNeiK5kT5ldEetMRgNhJIFBqykJUtYL1Suay7bdVwK\nNmMM5gOoWShkzH1S6chVYK8hJRAdtku8E4RDxhjMh1AyrgYAQCHfmhyYpDpF45emakHKTgyGFihb\nuMtgp06byn6jIjVLOcPMEnHIGIORDIUL8E14UfEqCiYLnC5dEuDdPjEYKVCaSIP50H/TuETKdUyb\n5r4wcaramgSD0TqoTYijz0npNphDRfLiJmp3OFU/CjBSZ6NgMFoGZX6uRgQrhTwKdy2NWla8wxPv\nM4bByIL6BHWsp19/xsqkSFgee9LQRDY2YwxGJipINGnMfVDoQEkwWcR3ehLD76z2RsFgtAzqx9UA\nAChkjcXB0WSlJPg82fwtjjVhMHJRTeJnBjt15Bh2ASkZFeuHdkn2w2aMwchHVQncTYMi42wPk+js\nBfahAkFPWtoEg9E2VLcRw9iELb7jk5S8OH2G+yePp9PTJBiM1qHCDVX0OY/bDuZUKHFlLW9ARQKv\nLW2NgsFoGarxczUiWK6/b4aiF930yftpCY41YTCEUfEGZ6xkj1mslwpdUuQ9wTnVC5sxBkMcVW9U\n2DYgtsBhB/Ft3RDf9m40vwu9jYLBaBkqHlcDAKAQX8tDI4mVfeTzaO0GA5rbBIPRNtSwcTCDnTjA\nhV1IoGQld1jnZC42YwxGQdTQHwMA3PAp2il3tYNghWgvhSssMJj/DGrojwEAxj1a7T0hWWaRrLnu\ns5OxGWMwSqAmOwZ9v0SjgeurpZ6v5dnnx/Pa0d0cGIxWoqZxNQAACJYZBk6TfOqWT+6Wlep6pmAw\nrQ112g4rcaYbK0vCiSLv8UNSONiMMRglUavxdOQ9yHfivR9MRny7G9f4XeluCQxGe1HnuBoAQHR0\nrc3BEc2PPFmagEPGGAwp1D2Y1fFKcRrtXdb4uZI71DQJh4wxGFKouz8GAIhcVrJDvE2TYKXwFw+6\n2wCD0XbocC5NiPf2mpgCANkecz5OxmaMwZCFFiexETdRfzD3Lc8+9yGvPd0tgMFoP3SMqwEA0Ik1\nZcY7v8TLEzEYCqDLjgGKjy3CsSYMhhLos2MMBkMVeBIVBqP9YDvGYLQfbMcYjPbzf7VOsJw7Kts8\nAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDIyLTAxLTI1VDEzOjUyOjM4KzAwOjAws0zauAAAACV0RVh0\nZGF0ZTptb2RpZnkAMjAyMi0wMS0yNVQxMzo1MjozOCswMDowMMIRYgQAAAAASUVORK5CYII=\n" } }, "id": "bc7e45da-6dda-413a-adf7-c324c0f68c86" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# creating river data\n", "# data.frame object\n", "river_names = data.frame( \n", " river_name = c(\"Tejo\", \"Mondego\", \"Douro\")\n", ")\n", "\n", "# specify points as c(lon, lat) coordinates of rivers\n", "t_coord <- st_linestring(rbind(c(-10, 37), c(-9.1, 38.7), c(-9.1, 39.3))) \n", "m_coord <- st_linestring(rbind(c(-8, 38.5), c(-8.4, 40.2), c(-7.8, 41.9)))\n", "d_coord <- st_linestring(rbind(c(-8, 42),c(-8.6, 41.1), c(-8.8, 40)))\n", "\n", "# create a column with the specified CRS\n", "rivers <- st_sfc(t_coord, m_coord, d_coord, crs = \"EPSG:4326\") \n", "\n", "# create `sf` object with the data.frame and the sfc\n", "river_data <- st_sf(river_names, river_geometry = rivers) \n", "plot(river_data)" ], "id": "298e942c-8fc2-4ddb-b941-353622fb84ad" }, { "cell_type": "markdown", "metadata": {}, "source": [ "One way of joining the data is to simply attach columns with `cbind()`." ], "id": "4e98705e-8b94-4767-8013-3a3abcb2a003" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "portugal_data <- cbind(city_data, river_data)\n", "portugal_data" ], "id": "067bbe5c-cda4-444e-8030-f4c40fc838f0" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way to join datasets is through geospatial operations. These\n", "operations allow us to filter, merge, or subset dataframes based on\n", "whether points, lines, and polygons touch, intersect, or contain other\n", "geospatial objects.\n", "\n", "For example, let’s say we don’t know which river borders which city and\n", "we want to use geospatial operations to find out. We can use `st_join()`\n", "to merge both datasets based on intersections. `st_join()` is basically\n", "a left join in which the matches are determined by the relationships\n", "between the spatial objects." ], "id": "3c1cbced-80b0-4b55-92a8-864b5a817c8d" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# merge city_data with river_data \n", "pt_data <- st_join(city_data, river_data[\"river_name\"]) # specify \"river_name\" to keep it in the merged dataset\n", "pt_data" ], "id": "df8b3972-2450-4099-9239-db5363275e0b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "`st_join()` assigns the names of the rivers to the observations in which\n", "the point (i.e., city) intersects with the line (i.e., river).\n", "\n", "> **Note**: `st_join()` drops the geometry of the second object.\n", "\n", "Similar to `st_join()`, we can filter, merge, or subset `sf` objects\n", "with other geospatial functions. The most common types of geospatial\n", "operations are:\n", "\n", "- `st_intersects()`\n", "- `st_disjoin()`\n", "- `st_contains()`\n", "- `st_overlaps()`\n", "- `st_within()`\n", "\n", "You can learn more about the different geospatial operations in the\n", "[documentation of the\n", "package](https://r-spatial.github.io/sf/reference/st_join.html).\n", "\n", "### Test your Knowledge\n", "\n", "Create a polygon with vertices (0,0), (0,3), (3,6), (6,3), (6,0). Your\n", "answer must be a `sfg` object." ], "id": "e471de7c-3907-40f1-b328-8828dd5959c7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "polygon_1 <- st_...(...(...(...)))\n", "\n", "answer_1 <- polygon_1\n", "\n", "test_1()" ], "id": "e38601d1-4c11-44a0-be3a-e1848dec0206" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Determine whether the following points are contained in the polygon with\n", "a geospatial operation." ], "id": "17d4d70c-bf04-4ff0-be1e-c47ae4c27a9e" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "point_1 <- st_point(c(0,0))\n", "point_2 <- st_point(c(3,4))\n", "point_3 <- st_point(c(3,3))\n", "point_4 <- st_point(c(4,5))\n", "point_5 <- st_point(c(6,6))\n", "points <- st_sfc(point_1, point_2, point_3, point_4, point_5)" ], "id": "8a302fd1-3b67-4070-b4dc-aa3c848a3760" }, { "cell_type": "markdown", "metadata": {}, "source": [ "> **Hint**: you only need to run one function." ], "id": "d46872b8-4745-4202-856e-3ea06b569bb7" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "answer_2 <- st_...(...)\n", "\n", "test_2()" ], "id": "d531534b-74b4-4d3a-9832-5ed314ce7ab5" }, { "cell_type": "markdown", "metadata": {}, "source": [ "> **Think Deeper**: once you get to the correct answer of question 2,\n", "> inspect the output of `answer_2`. What does it mean? How is it\n", "> different than the intersection of two objects?\n", "\n", "Which of the following functions should be used to create an object with\n", "a simple feature column and other attributes?\n", "\n", "1. `as_data_frame()`\n", "2. `st_sfg()`\n", "3. `st_sfc()`\n", "4. `st_sf()`" ], "id": "d2ead71a-002d-4260-bb3c-a264bdf5aaf5" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Replace ... by your answer of \"A\", \"B\", \"C\" or \"D\"\n", "answer_3 <- \"...\"\n", "\n", "test_3()" ], "id": "2c91f3f6-8dc2-4cf7-9eb9-3def02bdc1da" }, { "cell_type": "markdown", "metadata": {}, "source": [ "What are the common coordinate inputs for geographic CRS?\n", "\n", "1. `c(lat, lon)`\n", "2. `c(lon, lat)`\n", "3. `c(x, y)`\n", "4. `c(y, x)`" ], "id": "cd451a17-472e-427d-84d0-cc2692d3fcd5" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Replace ... by your answer of \"A\", \"B\", \"C\" or \"D\"\n", "answer_4 <- \"...\"\n", "\n", "test_4()" ], "id": "109c07d7-0162-4d5c-bed4-1a0780a1e9ee" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2: Hedonic Pricing with Data from Athenian Properties\n", "\n", "In this section, we’re going to estimate the value of apartment\n", "characteristics with a hedonic pricing function for apartments located\n", "in Athens, Greece. We’ll use the geospatial datasets `depmunic` and\n", "`properties` from the package `spData`.\n", "\n", "### Theory\n", "\n", "Typically, economists rely on market prices to estimate consumers’\n", "marginal willingness to pay for a good. When the good in question is not\n", "traded in markets, economists must resort to other strategies. One of\n", "those strategies is the **hedonic pricing method**.\n", "\n", "The main assumption of hedonic pricing is that *the total price of a\n", "good equals the sum of the prices of its components*. Using this\n", "assumption, we can deconstruct the price of a good traded in markets and\n", "attribute part of the total price to each of its non-traded component\n", "parts. For example, we can model the price of an apartment as a function\n", "of it’s size, view, age, location, number of bathrooms, proximity to\n", "public service, etc. If we take two apartments that are exactly equal\n", "but only differ on their view (e.g., one might be in the 2nd floor and\n", "the other on the 30th floor of the same building), their difference in\n", "price must be the added value of the view. If we have data on the prices\n", "and attributes of several apartments, we could leverage hedonic pricing\n", "to estimate how much consumers value a nice view.\n", "\n", "> **Think Deeper**: what are the other underlying assumptions of hedonic\n", "> pricing? When might it not be appropriate to use this method to value\n", "> non-traded goods?\n", "\n", "A model for the hedonic price of an apartment could be estimated with\n", "the following regression:\n", "\n", "$$\n", "P_{i} = \\beta_{0} + \\beta_{1} Size_{i} + \\beta_{2} Age_{i} + \\beta_{3} Bath_{i} + \\beta_{4} View_{i} + \\beta_{5} Room_{i} + \\beta_{6} Park_{i} + ... + \\beta_{k} X_{ki}\n", "$$\n", "\n", "where\n", "\n", "- $P_{i}$ is the price per square meter of apartment $i$\n", "- The variables $Size$, $Age$, $Bath$, $View$, $Room$, $Park$, $...$,\n", " $X_{k}$ are the components of the apartment\n", "- The $\\beta_{k}$s are the effects of those components on the price of\n", " the apartment\n", "\n", "With that in mind, let’s try to estimate the hedonic pricing function of\n", "apartments in Athens.\n", "\n", "### Exploratory Analysis\n", "\n", "Let’s load two datasets from the `spData` package.\n", "\n", "- `depmunic` has information on the 7 districts of the municipality of\n", " Athens. The district boundaries are stored in the data as polygons.\n", "- `properties` has data on properties for sale in the city of Athens.\n", " Property locations are stored in the data as points.\n", "\n", "Let’s start with `depmunic`." ], "id": "980af79e-ae7d-40e6-a79c-5b79cc1652b6" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "head(depmunic)\n", "str(depmunic)" ], "id": "b15fd0ce-8d89-42b9-ab16-a1a84fa8e978" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the dataset has attributes of the districts: `airbnb`,\n", "`museums`, `population`, …\n", "\n", "We can also see that the CRS used for this dataset is a Greek CRS:\n", "Projected CRS: GGRS87 / Greek Grid. This tells us that we’re dealing\n", "with a projected CRS tailored to Greece. A quick search will find that\n", "the units of reference are in meters.\n", "\n", "> **Note**: a crucial part of geospatial data exploration is learning\n", "> about the CRS of the dataset. This not only affects the relationships\n", "> within your own dataset but also determines whether you can merge or\n", "> filter geospatial objects based on other datasets.\n", "\n", "Let’s plot the `depmunic` data. We specify the attribute greenspace\n", "`[\"greensp\"]` to plot a chart of the green space exclusively." ], "id": "bb87cc2a-a58e-4c2f-b9dd-20ea4b6cc786" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot(depmunic[\"greensp\"])" ], "id": "89d27b83-49f3-44e2-8930-8fbd3591addc" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot shows that the district in yellow (district 2) is the one with\n", "largest greenspace, almost 10 times as much as the district in dark blue\n", "(district 4).\n", "\n", "Now let’s take a look at our `properties` data." ], "id": "6b42f69a-f12e-44c1-862f-a211932dbd94" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "head(properties)\n", "str(properties)" ], "id": "6610e21b-29b7-422d-a064-a66336621c98" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the dataset has attributes of the properties: `size`,\n", "`price`, `age`, …\n", "\n", "We can see that the CRS used for this dataset is also the Greek CRS.\n", "This is good news: both our datasets are referenced to the same grid.\n", "\n", "Let’s plot our `properties` dataset." ], "id": "6318bf19-eaaa-47e9-ab25-9d5d04c0d766" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot(properties[\"prpsqm\"])" ], "id": "db8124e8-51b9-43e6-a61a-6d2c97a34602" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot shows the distribution of properties with their color\n", "representing the price per square meter `prpsqm`. It’s hard to see a\n", "pattern from this chart, but it seems that most of the expensive\n", "properties are located towards the south.\n", "\n", "Let’s plot our property locations on top of our district boundaries." ], "id": "da16fc2b-6108-44a8-8184-c75b5d9e17b3" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot(depmunic[\"num_dep\"], reset = FALSE)\n", "plot(properties[\"price\"], add = TRUE, col = \"darkgreen\")" ], "id": "172d4cdc-2504-40c4-8589-3a2bbf64bf3f" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Empirical Methods\n", "\n", "Let’s try to put together a model to explain the price of properties in\n", "Athens using the hedonic pricing method. Our proposed specification will\n", "be: $$\n", "P_{i} = \\beta_0 + \\beta_1 Size_{i} + \\beta_2 Age_{i} + \\beta_{3} Dist_{i} + \\sum_{j=2}^{7} \\gamma_{j} D_{ij} + \\epsilon\n", "$$\n", "\n", "where\n", "\n", "- $P_{i}$ is the price (in euros) per square meter of apartment $i$,\n", "- $Size_{i}$ is the size in square meters of apartment $i$,\n", "- $Age_{i}$ is the age in years of apartment $i$,\n", "- $Dist_{i}$ is the distance between apartment $i$ and the closest\n", " metro in meters,\n", "- $D_{ij}$ are dummy variables indexing the districts $j$ for\n", " apartment $i$.\n", "\n", "The dummy variables are included to account for district-specific\n", "characteristics, such as the number of museums, parks, area of\n", "greenspace, and other district specific factors that we don’t have data\n", "on.\n", "\n", "Looking at our `properties` dataset, we see that we already have data on\n", "all of our regressors except district location of the apartments. We\n", "need to use what we learned about geospatial operations to assign a\n", "categorical variable indicating district location to each property.\n", "\n", "Below we use the function `st_join()` to merge both datasets based on\n", "the intersection of the spatial elements, and pass the attribute\n", "`num_dep` (the district number) to the merged dataset." ], "id": "77c7a054-b2c0-4add-892d-b3ca69445781" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# find location districts for each property\n", "merged_data <- st_join(properties, depmunic[\"num_dep\"])\n", "merged_data" ], "id": "e00e3b44-2fa0-487b-9561-2a611c90591b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that `st_join()` dropped the district boundaries. That is fine,\n", "since `num_dep` already indicates where each property is located. As a\n", "matter of fact, we don’t even need the points anymore. We only needed\n", "them to connect the apartments to their district! Let’s clean our data\n", "before running our model.\n", "\n", "> **Note**: we should drop spatial objects with the function\n", "> `st_drop_geometry()`, as opposed to `select()`." ], "id": "d39bd114-3eb9-4133-b60f-611e8536f1aa" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_data <- merged_data %>%\n", " select(size, prpsqm, age, dist_metro, num_dep)%>%\n", " mutate(num_dep = as_factor(num_dep))%>% # so that R understand that this is a categorical variable\n", " st_drop_geometry()\n", "str(model_data)" ], "id": "b199e23d-e6e8-42ad-8b8a-4b61fd5e4e5d" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results and Analysis\n", "\n", "Now, we’re ready to run our hedonic pricing model. Let’s do that using\n", "`lm()`. R will automatically create dummy variables for each district." ], "id": "0c369660-4190-4ca3-bb69-1e15bf559d52" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1 <- lm(prpsqm ~ size + age + dist_metro + num_dep, data = model_data)\n", "coeftest(model1, vcov = vcovHC)" ], "id": "5901e15d-1b34-4510-8122-94fbf4b0f450" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The most significant factor contributing to the price of an apartment in\n", "Athens appears to be district location. Although age and distance to\n", "metro are statistically significant, their effects are quite small (as\n", "long as the property is not *super* old or *super* far from the closest\n", "metro). However, we do see a wide range of effects for district\n", "locations. Once we control for size, age, and distance from metro, the\n", "mere fact of being in district 4 decreases the price per square meter by\n", "over 1,000 euros! Maybe it’s the lack of greenspace that we have noted\n", "earlier…\n", "\n", "Our model breaks down the price of an apartment as follows:\n", "\n", "$$\n", "P_{i} = 2352 + 4.3 Size_{i} -20 Age_{i} -0.1 Dist_{i} + \\sum_{i=2}^{7} \\gamma_{i} D_{i} + e\n", "$$\n", "\n", "In which the added effects for the dummies are:\n", "\n", "- District 2: $-334$ euros\n", "- District 3: $-490$ euros\n", "- District 4: $-1031$ euros\n", "- District 5: $-822$ euros\n", "- District 6: $-918$ euros\n", "- District 7: $-653$ euros\n", "\n", "District 1 is the base case (i.e., when all dummies equal zero), and is\n", "the most accretive location to the price of a property. That makes\n", "sense, since District 1 is the historic center of Athens.\n", "\n", "Now, let’s look at how well our model actually predicts the price per\n", "square meter of properties. This can shed light on the appropriateness\n", "of using hedonic pricing strategies for properties in Athens." ], "id": "64a7842f-5a9f-438c-8d65-1e0d8db99bff" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "summary(model1)$adj.r.squared" ], "id": "839eb9b4-471b-40cc-939e-d3f301f0fa4e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our adjusted $R^2$ says that our model only explains ~35% of the\n", "variation in the price of properties in Athens. The cause of low\n", "predictive power could be an error of methodology (prices of properties\n", "in Athens might not be the sum of their component parts!), or a missing\n", "variable problem.\n", "\n", "> **Think Deeper**: what component parts of properties might we be\n", "> missing in our hedonic pricing model?" ], "id": "5a628720-a863-46cc-b019-66aaa86716e2" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "ir", "display_name": "R", "language": "r" } } }