GeoPython2022.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9a96e449-efc1-4a59-8ab4-13fe52d309bf",
   "metadata": {},
   "source": [
    "# GeoPython 2022 - Introduction to `dask-geopandas`\n",
    "\n",
    "**Martin Fleischmann, Joris van den Bossche**\n",
    "\n",
    "22/06/2022, Basel\n",
    "\n",
    "## Setup\n",
    "\n",
    "Follow the Readme to set-up the environment correctly. You should have these packages installed:\n",
    "\n",
    "```\n",
    "- geopandas\n",
    "- dask-geopandas\n",
    "- pyogrio\n",
    "- pyarrow\n",
    "- python-graphviz\n",
    "- esda\n",
    "- dask-labextension # optionally, if using JupyterLab\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d85aacd2-6dff-4e19-a1b4-ec6201456d5d",
   "metadata": {},
   "source": [
    "## GeoPandas refresh\n",
    "\n",
    "Let's start with a quick refresh of GeoPandas.\n",
    "\n",
    "### What is GeoPandas?\n",
    "\n",
    "**Easy, fast and scalable geospatial analysis in Python**\n",
    "\n",
    "From the docs:\n",
    "\n",
    "> The goal of GeoPandas is to make working with geospatial data in python easier. It combines the capabilities of pandas and shapely, providing geospatial operations in pandas and a high-level interface to multiple geometries to shapely. GeoPandas enables you to easily do operations in python that would otherwise require a spatial database such as PostGIS.\n",
    "\n",
    "A quick demo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51fea161-12a7-423e-a152-1cd07402993a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76618132-2d59-4cb0-87a4-43e194aa7ded",
   "metadata": {},
   "source": [
    "GeoPandas includes some built-in data, we can use them as an illustration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a342e516-76ad-44c0-b5d5-895ce2ebf91c",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = geopandas.datasets.get_path(\"naturalearth_lowres\")\n",
    "path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc58b219-b3fc-4082-84ac-c2272fb587d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "world = geopandas.read_file(path)\n",
    "world.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a47d5e4-95fe-4161-ae7e-cf93e65b83b2",
   "metadata": {},
   "source": [
    "For the sake of simplicity here, we can remove Antarctica and re-project the data to the EPSG 3857, which will not complain about measuring the area (but never use EPSG 3857 to measure the actual area as it is extremely skewed)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98a2f3a2-a480-4405-859e-539fd156f465",
   "metadata": {},
   "outputs": [],
   "source": [
    "world = world.query(\"continent != 'Antarctica'\").to_crs(3857)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a3e24d4-4e1b-4f8d-aed0-16a255643080",
   "metadata": {},
   "source": [
    "GeoPandas GeoDataFrame can carry one or more geometry columns and brings the support of geospatial operations on these columns. Like a creation of a convex hull."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b6be908-ecc9-48eb-aff9-f079550518a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "world['convex_hull'] = world.convex_hull"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e660148-41db-41a7-8f83-dde4fa7735e8",
   "metadata": {},
   "source": [
    "This is equal to the code above as GeoPandas exposes geometry methods of the active geometry column to the GeoDataFrame level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb9af8d4-88cc-4fe6-966a-e9efaf6701dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "world['convex_hull'] = world.geometry.convex_hull  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4653ba04-06e8-4684-bdd0-c4775cf94baa",
   "metadata": {},
   "source": [
    "Now you can see that we have two geometry columns stored in our `world` GeoDataFrame but only the original one is treated as an _active_ geometry (that is the one accessible directly, without getting the column first)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "386275f4-43a0-4ed2-a9b7-ebf53242f2ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "world.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2a7eba8-a4f3-4c47-b9d7-d7652678fac8",
   "metadata": {},
   "source": [
    "We can also plot the results. Both Russia and Fiji are a bit weird as they cross the anti-meridian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b33c8337-9bce-4478-aa8e-b041ded3b692",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = world.plot(figsize=(12, 12))\n",
    "world.convex_hull.plot(ax=ax, facecolor='none', edgecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa74cdbe-ddeb-42a5-b466-072c93a97e37",
   "metadata": {},
   "source": [
    "## What is Dask\n",
    "\n",
    "From the docs:\n",
    "\n",
    "> Dask provides advanced parallelism and distributed out-of-core computation with a dask.dataframe module designed to scale pandas. Since GeoPandas is an extension to the pandas DataFrame, the same way Dask scales pandas can also be applied to GeoPandas.\n",
    "\n",
    "We will cover the high-level API of Dask. For more, see the [Dask tutorial](https://tutorial.dask.org).\n",
    "\n",
    "Let's import `numpy` and `pandas` for a comparison and three high-level Dask modules - `bag`, `array`, and `dataframe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "900e0740-f592-438b-b0e3-ee20f26e0751",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import dask.dataframe as dd\n",
    "import dask.array as da\n",
    "import dask.bag as db"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6770195-2ad8-44f2-b3c3-d1c032ad78ed",
   "metadata": {},
   "source": [
    "Before we explore those, let's introduce the dask `Client` as it will allow us to see how dask manages all its tasks.\n",
    "\n",
    "Here we create a Client on top of a local (automatically created) cluster with 4 workers (the laptop we use has 4 performance cores)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a23995a-cf0b-47aa-9552-12ca69f0cbab",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client\n",
    "\n",
    "client = Client(n_workers=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cba16798-a0df-456f-8ff6-45bd6ce6348d",
   "metadata": {},
   "outputs": [],
   "source": [
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31300c5c-39f8-4997-b9c5-6a9fdefea050",
   "metadata": {},
   "source": [
    "We can open the Dask dashboard to watch what is happenning in real-time using the link above, in the Client details. But if you have a [JupyterLab extension for Dask](https://github.com/dask/dask-labextension), you can watch different components directly from the JupyterLab interface."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63f301b8-fe13-4be3-b75b-0541b91b7355",
   "metadata": {},
   "source": [
    "### dask.bag\n",
    "\n",
    "With the Client and cluster in place, we can properly explore Dask. Let's start with a `dask.bag`, the simplest of the objects. You can imagine it as a distributed list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65ae1ed7-cde5-444a-a0f0-fbb8dd4feaaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = db.from_sequence([1, 2, 3, 4, 5, 6, 2, 1], npartitions=2)\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf414b2e-405b-4d5a-ab10-24a9cd7b3c5a",
   "metadata": {},
   "source": [
    "Now, note that when we try to call `b`, we don't see its contents. \n",
    "\n",
    "Let's check what happens with `sum`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "011741f8-d169-4893-9a17-0f4f72f6383d",
   "metadata": {},
   "outputs": [],
   "source": [
    "b.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01f7f777-c7a1-428e-bc07-c4f97f2c2b21",
   "metadata": {},
   "source": [
    "Again, we don't see the answer, but some abstract `Item` object instead. That is because Dask usually runs all the operations lazily and waits for a `compute` call before it does the actual computation.\n",
    "\n",
    "Instead, it plans what it should do and create a task graph. That looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e12a2d2b-97f4-43cf-951e-5272bc33370a",
   "metadata": {},
   "outputs": [],
   "source": [
    "b.sum().visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4a4ae09-8f2e-4f1e-88ae-aad41ec2f7ec",
   "metadata": {},
   "source": [
    "We can see individual partitions (rectangles), operations (circles), and movement of data between partitions.\n",
    "\n",
    "When we call `compute`, this task graph is executed and Dask returns the expected value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2733578-508d-4b45-b976-3992d0442e33",
   "metadata": {},
   "outputs": [],
   "source": [
    "b.sum().compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89de6a07-0c94-4f18-9435-e54dd7751831",
   "metadata": {},
   "source": [
    "### dask.array\n",
    "\n",
    "Let's move onto an array. Where bag is partitioned along 1 dimension (the sequence is essentially cut into pieces), array is like a numpy array split along both dimension. In practice, each of the partitions is a numpy array and dask array combines them together. Each partition can be then processed separately. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b204f4c-522c-487b-9433-82d68a0c65df",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.arange(100_000).reshape(200, 500)\n",
    "a = da.from_array(data, chunks=(100, 100))\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e642e71-9cbc-4244-ad00-8a5413b94ac1",
   "metadata": {},
   "source": [
    "We see some dimensions, dtypes and sizes here but not the data. Beacause again, all is done lazily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff387bb9-d49c-4797-abf5-cb29272ea77d",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[:50, 200]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05fd4c16-13b9-4faa-b7df-96d7b36d1b7f",
   "metadata": {},
   "source": [
    "Even indexing requires `compute` to return values, otherwised it still give a dask array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20c93739-1a76-4b5d-a9ef-6b86d72050c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[:50, 200].compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "346a5b76-622a-41c5-8a86-1ec9e611fe31",
   "metadata": {},
   "source": [
    "Simlarly for `mean`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c85e22a-6696-4e78-8769-53a14392e68b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a.mean().compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9639daa-a1ff-4524-98b4-def76aa5cda9",
   "metadata": {},
   "source": [
    "Since the mean is not super straightforward to parallelise (you can't just do mean in partitions and then mean of that), we can check how dask implements its logic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab7ab907-24ff-4d67-939b-efc36a16d216",
   "metadata": {},
   "outputs": [],
   "source": [
    "a.mean().visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40a4fc45-224b-4002-aaa3-9732f0e06481",
   "metadata": {},
   "source": [
    "Quite complex, right? Let's compare it to the indexing we did before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0127d06f-6e2d-46cd-9bd5-e38f297fe0e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[:50, 200].visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "761e124c-f672-4731-a7ad-c5ed93a39f37",
   "metadata": {},
   "source": [
    "You can see that dask efficiently accesses only that one partition it needs at this moment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26d6f86-7813-4664-85fc-e8209919649d",
   "metadata": {},
   "source": [
    "### dask.dataframe\n",
    "\n",
    "Finally, we move to the parallelised DataFrame. It mirrors the logic of the array implementation, with a difference that individual partitions are pandas.DataFrames and partitioning happens along a single axis (rows).\n",
    "\n",
    "Dask.dataframe tries to mirror the pandas API. The same approach as we will later see with dask-geopandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6b42178-4360-49b5-aeb2-b7931308c7dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"data/airports.csv\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3108a874-251a-4e24-863b-cf3ab400a712",
   "metadata": {},
   "source": [
    "In this specific case (`head()`), dask actually reads those 5 rows and shows them but that tends to be an exception, likely because it is a very cheap operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a80274cc-d68b-4720-9bf3-6a5937d2f5c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = dd.read_csv(\"data/airports.csv\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee579606-8519-4d1e-a177-d4f1007aef18",
   "metadata": {},
   "source": [
    "If you try to show the whole DataFrame, you get a placeholder that tells you how many partitions you have, which columns and what are their dtypes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "532fd533-6952-4924-903a-de3bd1568387",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c925fb9-7bbf-46cc-b5ab-27903a9712de",
   "metadata": {},
   "source": [
    "Since the `airports.csv` is a single file on disk, dask gives us a single partion. To create a partitioned data frame, we can repartition it and even save to a partitioned CSV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6368cef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.repartition(4).to_csv(\"data/airports_csv/*.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7a2e2fc",
   "metadata": {},
   "source": [
    "When we have more of CSV files in a folder, typically one per month or a country, we can read each as a partition. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "201f882e-73aa-4ace-bded-cf9401b365af",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = dd.read_csv(\"data/airports_csv/*.csv\")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff0b7b08-a1f3-43cd-b3d9-bfaf40589501",
   "metadata": {},
   "source": [
    "As before, all the computation is done lazily. Take a look at the computation of mean elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65a5188e-8150-4535-b52a-203d97d5b0af",
   "metadata": {},
   "outputs": [],
   "source": [
    "elevation = df.elevation_ft.mean()\n",
    "elevation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba44ca44-fe9a-499a-843c-3beb62ea3844",
   "metadata": {},
   "source": [
    "We get a dask Scalar object here but as before, we don't get the value until we call `compute()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a08fc50-897a-4fc5-9712-79d3988ce9aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "elevation.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74b81842-28b9-433d-b7c6-9558362dd831",
   "metadata": {},
   "source": [
    "You can probably notice the similarity of the graph with the one calculating mean over an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81839c97-eea2-4969-b3a2-e25d50191382",
   "metadata": {},
   "outputs": [],
   "source": [
    "elevation.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08b11a5a-c9e9-4d23-983c-920a4c024ef3",
   "metadata": {},
   "source": [
    "We can also quickly compare it to a task graph for something easier, like a sum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c0a51c0-91c8-4cc3-92d9-5a2cf148610d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.elevation_ft.sum().visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c9341b2-9687-4d51-9d5b-5555be12c114",
   "metadata": {},
   "source": [
    "Not that it makes much sense to compute a sum of elevations, but we can do that and if you're checking the dashboard, you'll notice very little communication as the task is easy to parallelise and we need to gather the results only in the final step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1ea7346-17c1-402e-9e4b-76d6aa336c04",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.elevation_ft.sum().compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76750d1b-b79d-436d-ad8a-c725876bedcf",
   "metadata": {},
   "source": [
    "## Dask-GeoPandas\n",
    "\n",
    "Dask-GeoPandas follows exactly the same model as `dask.dataframe` adopted for scaling `pandas.DataFrame`. We have a single `dask_geopandas.GeoDataFrame`, composed of individual partitions where each is a `geopandas.GeoDataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c83221e2-8016-409e-be3d-85fe4b0857cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask_geopandas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fa39737-d827-4421-9159-5268028f562d",
   "metadata": {},
   "source": [
    "## Create dask GeoDataFrame\n",
    "\n",
    "We have a plenty of options how to build a `dask_geopandas.GeoDataFrame`. From in-memory `geopandas.GeoDataFrame`, reading the GIS file (using pyogrio under the hood), reading GeoParquet or Feather, or from dask.dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f51ae5-fc1e-4fbd-a28c-a04ebad8ea77",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf = dask_geopandas.from_geopandas(world, npartitions=4)\n",
    "world_ddf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39216390-63b5-4a87-897a-8e65fca05b0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf_file = dask_geopandas.read_file(path, npartitions=4)\n",
    "world_ddf_file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97f3020d-bb99-4416-a582-11da211a8fb8",
   "metadata": {},
   "source": [
    "### Partitioned IO\n",
    "\n",
    "Since we are working with individual partitions, it is useful to save the dataframe already partitioned. The ideal file format for that is a GeoParquet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ec8c53c-6cc0-40b6-b9e1-17b786c22fec",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world_ddf.to_parquet(\"data/world/\")\n",
    "world_ddf.to_crs(4326).to_parquet(\"data/world_4326/\")  # later we will need the dataset in EPSG:4326 so we can already prepare it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "811b8855-9298-4835-b15e-f34ca9715cd5",
   "metadata": {},
   "source": [
    "For more complex tasks, we recommend using Parquet IO as an intermediate step to avoid large task graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6121085-a81c-47dd-bf15-cae358fe91bc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world_ddf = dask_geopandas.read_parquet(\"data/world/\")\n",
    "world_ddf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55470865-1758-475a-af94-c26a222ce1a3",
   "metadata": {},
   "source": [
    "## Embarrassingly parallel computation\n",
    "\n",
    "The first type of operations where you can benefit from parallelisation is so-called embarrassingly parallel computation. That is a computation where we treat individual partitions or individual rows indenpendently of the other, meaning there is no inter-worker communication and no data need to be sent elsewhere.\n",
    "\n",
    "One example of that is a calculation of area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6edb7684-32e4-495f-a1ec-74aca19866da",
   "metadata": {},
   "outputs": [],
   "source": [
    "area = world_ddf.area\n",
    "area.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03224973-bfce-4d0b-bf1a-fc366284b233",
   "metadata": {},
   "source": [
    "Similar one, this time returing a `dask_geopandas.GeoSeries` instead of a `dask.dataframe.Series` would be a `convex_hull` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38da2d25-5df0-4903-a438-13d59897d59e",
   "metadata": {},
   "outputs": [],
   "source": [
    "convex_hull = world_ddf.convex_hull\n",
    "convex_hull.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31ee2c9f-9336-46da-9f6f-49e8f572a8d0",
   "metadata": {},
   "source": [
    "Since both are creating a series, we can assing both as individual columns. Let's see how that changes the task graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1e2d93d-8d04-4c80-bc00-e43b49b9a9a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf['area'] = world_ddf.area\n",
    "world_ddf['convex_hull'] = world_ddf.convex_hull\n",
    "world_ddf.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19bc9cf6-a2be-4674-926f-f8efe04cf18f",
   "metadata": {},
   "source": [
    "Finally, we can call `compute()` and get all the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "881ee02b-8928-46fd-ba8e-188705c6b758",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "r = world_ddf.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14a463d8-1f08-4995-9394-2112d5beeafe",
   "metadata": {},
   "source": [
    "## Spatial join\n",
    "\n",
    "If you have to deal with a large dataframes and need a spatial join, dask-geopandas can help. Let's try to illustrate the logic of spatial join on the partitioned data using the locations of airports from around the world."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13ff8202-2d2f-454e-805a-a091eb19d828",
   "metadata": {},
   "outputs": [],
   "source": [
    "airports = pd.read_csv(\"data/airports.csv\")\n",
    "airports.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "851ef74b-a48c-4cfd-96ab-b100c3c4584b",
   "metadata": {},
   "source": [
    "The data comes as a CSV, so we first need to create a GeoDataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eec8022-19cc-45b3-a1d4-00dbbe55e983",
   "metadata": {},
   "outputs": [],
   "source": [
    "airports = geopandas.GeoDataFrame(\n",
    "    airports,\n",
    "    geometry=geopandas.GeoSeries.from_xy(\n",
    "        airports[\"longitude_deg\"],\n",
    "        airports[\"latitude_deg\"],\n",
    "        crs=4326,\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a36783f2-eef2-4a81-b771-329cc5e2dc95",
   "metadata": {},
   "source": [
    "And from that, we can create a partitioned `dask_geopandas.GeoDataFrame`. Note that we could also read the CSV with dask.dataframe and create a GeoDataFrame from that using the `dask_geopandas.from_dask_dataframe` function and `dask_geopandas.points_from_xy` to create the geometry. But since it all comfortably fits in memory, we can pick whichever option we like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92417d31-ff01-4059-b166-88114baca621",
   "metadata": {},
   "outputs": [],
   "source": [
    "airports_ddf = dask_geopandas.from_geopandas(\n",
    "    airports,\n",
    "    npartitions=12\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca1701ef-d33c-44a4-b81a-c7abeb352503",
   "metadata": {},
   "source": [
    "We will join the point data of airports with the `naturalearth_lowres` dataset we have stored as an already partitioned parquet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ec0cb6f-1ae5-43b1-8f0a-77f18b1c68c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf = dask_geopandas.read_parquet(\"data/world_4326/\")\n",
    "world_ddf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cbce801-7bf6-40be-b431-d8b68b870ca3",
   "metadata": {},
   "source": [
    "The API of the `sjoin` is exactly the same as you know it from geopandas. Just in this case, it currently only creates a task graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae07e6dc-5025-4333-8d0f-95cf6f5df82e",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined = airports_ddf.sjoin(world_ddf, predicate=\"within\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38b8c371-2051-44a8-9d70-2470a0eed316",
   "metadata": {},
   "source": [
    "We started from 12 partitions of `airports_ddf` and 4 partitions of `world_ddf`. Since we haven't told Dask how are these partitions spatially distributed, it just plans to do the join of each partition from one dataframe to each partition form the other one. 12x4 = 48 partitions in the end. We can easily check that with the `npartitions` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5261936f-2929-4e6a-8637-6eaba0019363",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined.npartitions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f06bbc8-aa26-4bb3-88f2-804814f24209",
   "metadata": {},
   "source": [
    "The whole logic can also be represented by a task graph that illustates the inneficiency of such an approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d382eb3-2b68-46e3-93ec-f74cd64bdd14",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d097ceee-e0bf-4451-84a1-db9057dca405",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Spatial partitioning\n",
    "\n",
    "Luckily, dask-geopandas supports spatial partitioning. It means that it can calculate the spatial extent of each partition (as the overall convex hull of the partition) and use it internally to do smarter decisions when creating the task graph. \n",
    "\n",
    "But first, we need to calculate these paritions. This operation is done eagerly and involves immediate reading of all geometries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdcb758e-1ad0-4d04-ba58-26ce04dca646",
   "metadata": {},
   "outputs": [],
   "source": [
    "airports_ddf.calculate_spatial_partitions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac742c14-ee33-4df9-a037-cdbbdda6c16d",
   "metadata": {},
   "source": [
    "The resulting `spatial_partitions` attribute is a `geopandas.GeoSeries`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4c7de49-1798-4d3e-b479-f18f0fc149e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "airports_ddf.spatial_partitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ee6da69-67c8-4ec5-94c4-c8651b88a2da",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "airports_ddf.spatial_partitions.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "457e6d7a-b697-4215-8011-97faa89dcf48",
   "metadata": {},
   "source": [
    "As you can see from the plot above, our partitions are not very homogenous in terms of their spatial distribution. Each contains points from nearly whole world. And that does not help with simplification of a task graph."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5034996f-e493-4552-bf43-f7d30e49b06a",
   "metadata": {},
   "source": [
    "### The goal\n",
    "\n",
    "We need our partitions to be spatially coherent to minimise the amount of inter-worker communication. So we have to find a way of reshuffling the data in between workers.\n",
    "\n",
    "### Hilbert curve\n",
    "\n",
    "One way of doing so is to follow the Hilbert space-filling curve, which is a 2-dimensional curve like the one below along which we can map our geometries (usually points). The distance along the Hilbert curve then represents a spatial proximity. Two points with a similar Hilbert distance are therefore ensured to be close in space.\n",
    "\n",
    "![Hilbert](fig/Hilbert-curve_rounded-gradient-animated.gif)\n",
    "\n",
    "(Animation by Tim Sauder, https://en.wikipedia.org/wiki/Hilbert_curve#/media/File:Hilbert-curve_rounded-gradient-animated.gif)\n",
    "\n",
    "`dask-geopandas` (as of 0.1.0) implements Hilbert curve and two other methods based on a similar concept of space-filling (Morton curve and Geohash). You can either compute them directly or let `dask-geopandas` use them under the hood in a `spatial_shuffle` method that computes the distance along the curve and uses it to reshuffle the dataframe into spatially homogenous chunks. (Note that geometries are abstracted to the midpoint of their bounding box for the purpose of measuring the distance along the curve.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "695f0568-7fbb-4e30-b4e8-cd919c2edb89",
   "metadata": {},
   "outputs": [],
   "source": [
    "hilbert_distance = airports_ddf.hilbert_distance()\n",
    "hilbert_distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3466b3b7-c3f0-4087-aa3b-cd1412482d90",
   "metadata": {},
   "outputs": [],
   "source": [
    "hilbert_distance.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae0bd3e-8a27-4167-a066-6999b41c9496",
   "metadata": {},
   "source": [
    "`spatial_shuffle` uses by default `hilbert_distance` and partitions the dataframe based on this Series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e48c50ac-d090-4bde-90d9-c5d3df3e9774",
   "metadata": {},
   "outputs": [],
   "source": [
    "airports_ddf = airports_ddf.spatial_shuffle()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "357fbb71-e8dd-458a-bb5e-eb70ed6f53e9",
   "metadata": {},
   "source": [
    "We can now check how the new partitions look like in space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec44c64-d875-4d71-9256-7b2a4fc7908d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "airports_ddf.spatial_partitions.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2629430-2ef0-48b1-aeb3-b1fe67a032bc",
   "metadata": {},
   "source": [
    "When we are reading parquet file, its metadata already contain the information on the extent of each partition and we therefore don't have to calculate them by reading all the geometries. We can quickly check that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "214a5f75-8199-412a-a2ca-f03e27040a65",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf.spatial_partitions is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7902b24-d1b6-476d-8023-1129b7f71948",
   "metadata": {},
   "source": [
    "The world dataset has known partitions but is not spatially shuffled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0a7f683-7e16-45e6-9018-50fcdc0c888d",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf.spatial_partitions.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32806bca-85ba-473e-b397-87fb9d563c06",
   "metadata": {},
   "source": [
    "Even without doing that, we can already see that the resulting number of partitions is now 33, instead of 48 as some of the joins that would result in an empty dataframe are simply filtered out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "438faef8-657c-4c0b-94d4-7f7b4e4793f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined = airports_ddf.sjoin(world_ddf, predicate=\"within\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d4033ae-0a22-435c-9a3e-7409bbee9900",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined.npartitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f72dcd13-d4aa-4f02-9e22-07b27f84586d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "joined.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d812565a-a6df-41b7-883f-0bc2a74336ca",
   "metadata": {
    "tags": []
   },
   "source": [
    "### What about a larger problem?\n",
    "\n",
    "Dropping down from 48 to 33 partitions doesn't sound like a big deal. But you usually want to use dask-geopandas to tackle a bit larger problems. To simulate one, We can load a GADM dataset containd detailed administrative boundaries of the whole world (around 2GB GPKG) and join our airport data to that.\n",
    "\n",
    "_Note that this dataset is not part of the repository so you will not be able to run these two cells._\n",
    "\n",
    "_If you want to run the code you need to download the `gadm404.gpkg` from [GADM.org](https://gadm.org/download_world.html) and unzip it to the `data` folder._\n",
    "\n",
    "_To create the `gadm_spatial` used below, you need to read the GPKG, spatially shuffle it and save it as a GeoParquet:_\n",
    "\n",
    "```\n",
    "gadm_ddf = dask_geopandas.read_file('data/gadm404.gpkg', npartitions=64)\n",
    "gadm_ddf.spatial_shuffle().to_parquet(\"data/gadm_spatial/\")\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21dac6f2-8b94-4c73-b776-d5da8b390f35",
   "metadata": {},
   "outputs": [],
   "source": [
    "gadm_ddf = dask_geopandas.read_file('data/gadm404.gpkg', npartitions=64)\n",
    "joined = airports_ddf.sjoin(gadm_ddf, predicate=\"within\")\n",
    "joined.npartitions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cccff35-1237-4bf6-b11a-9363ffb9cc02",
   "metadata": {},
   "source": [
    "Without any spatial sorting of the GADM dataset, we have to do 12*64 joins.\n",
    "\n",
    "But we can load the same dataset that has been spatially sorted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d37393aa-764b-44eb-a1d1-10e6fd1ddb2c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gadm_sorted = dask_geopandas.read_parquet(\"data/gadm_spatial/\")\n",
    "joined = airports_ddf.sjoin(gadm_sorted, predicate=\"within\")\n",
    "joined.npartitions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c4da771-eab8-4f05-900c-64f7924f9590",
   "metadata": {},
   "source": [
    "The resulting number of partitions is 151, filtering out more than 80% of spatial joins that no longer need to be done."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ca7bb1b-1ebd-4aa0-ae7d-eb4bf934a4ca",
   "metadata": {},
   "source": [
    "## Aggregations with dissolve\n",
    "\n",
    "Dissolve is a typical operation when you usually need to do some shuffling of data between workers to ensure that all observations within the same category (specified using the `by` keyword) end up in the same partition so they can actually be dissolved into a single polygon. As you can imagine, proper spatial partitions may help as well but in this case, they help only in the computation, not in the task graph.\n",
    "\n",
    "Let's have a look at an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b90d0e1-6b9b-4281-ba2f-136a608ed4ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf = dask_geopandas.read_parquet(\"data/world/\")\n",
    "\n",
    "continents = world_ddf.dissolve('continent', split_out=6)\n",
    "continents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0df856a7-00a9-494e-9697-c9fa3449fe92",
   "metadata": {},
   "source": [
    "Above, we are using the API we know from GeoPandas with one new keyword - `split_out`. That specifies into how many partitions should we send the dissolved result. The whole method is based on the `groupby`, exactly as the original one, which returns a single partition by default. We rarely want that to happen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b308e36e-a6c4-4418-8f9f-eabdab337199",
   "metadata": {},
   "outputs": [],
   "source": [
    "continents.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c728957c-6199-43bf-98d4-7eb97d427fc6",
   "metadata": {},
   "source": [
    "The task graph shows exactly what happens. Since Dask doesn't know which categories are where, it designs the task graph to move potentially shuffle data from every original partition to every new one. In reality, some of these will be empty. And the better spatial partitions we have, the more of them will be empty, hence our operation will be cheaper."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c59ac64f-5a0c-407e-8c38-703888e6d900",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Custom functions with `map_partitions`\n",
    "\n",
    "Not every function you may need is built-in and you often need to apply a custom one to the partitioned dataframe. The most common way of doing that is the `map_partitions` method.\n",
    "\n",
    "The typical use case is below. We want to describe the shape of each polygon using the `shape_index` from the `esda` package. With geopandas, it would look like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41d81a3e-0f37-45e2-a94b-a9c1d2bd098d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from esda.shape import shape_index\n",
    "\n",
    "world['shape_idx'] = shape_index(world)\n",
    "world.explore('shape_idx')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a607879-a30c-436f-98d9-4b6d84ece225",
   "metadata": {},
   "source": [
    "But when you try to use the same code with a `dask_geopandas.GeoDataFrame`, it will not work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d0e3cc-133f-4816-974e-7ae46fb98922",
   "metadata": {},
   "outputs": [],
   "source": [
    "# THIS WILL FAIL\n",
    "world_ddf = dask_geopandas.read_parquet(\"data/world/\")\n",
    "world_ddf['shape_idx'] = shape_index(world_ddf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f588dd1c-ec69-4eb2-900c-e2d200d2ab49",
   "metadata": {},
   "source": [
    "In fact, it doesn't even give us a meaningful error message. Simply, because `esda` does not expect `dask_geopandas.GeoDataFrame` here, but expects `geopandas.GeoDataFrame` or some form of an in-memory array of geometries. Then it returns an array of floats.\n",
    "\n",
    "In this case, we can use `map_partitions` to _map_ the `shape_index()` function across individual partitions. A dummy code would look something like this:\n",
    "\n",
    "```py\n",
    "results = []\n",
    "\n",
    "for partition in ddf:\n",
    "    results.append(shape_index(partition))\n",
    "```\n",
    "\n",
    "The actual code is of course not a loop like this but principle is the same. We take each partition and apply the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5e3b06b-f2b6-4338-91ac-68f728f1a704",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf = dask_geopandas.read_parquet(\"data/world/\")\n",
    "shape_idx = world_ddf.map_partitions(shape_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cacb7765-be71-4945-9bde-a42f16e89c7f",
   "metadata": {},
   "source": [
    "The resulting task graph is the same as we have seen before in simple cases of embarrasingly parralel computation. `map_partitions` will always be embarrasingly parralel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9d310e9-3d5b-44af-abd1-f5e8ab193fec",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_idx.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1696723-3c49-4f2f-a371-4ee73349967c",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = shape_idx.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "420752cd-b2e5-4b45-8894-34f3e5c533d1",
   "metadata": {},
   "source": [
    "### Custom function\n",
    "\n",
    "We can also write our custom functions to be used with `map_partitions`. The only rule is that everything needs to happen within a single GeoDataFrame, i.e. within a single partition independently of the other. But we don't need to return a new column but we can also return a single value for each partition. Like a sum of area covered by polygons in each partition. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cd19de3-3ef0-4be7-af4a-110e9d2d8f2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_fn(gdf):\n",
    "    \"\"\"get a sum of area covered by polygons in a gdf\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    gdf : GeoDataFrame\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "    \n",
    "    \"\"\"\n",
    "    area = gdf.geometry.area\n",
    "    return sum(area)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd45ad4e-bb4d-4297-b7eb-7cb759b8eaee",
   "metadata": {},
   "source": [
    "We cannot assign this as a new column as we did above, but we also don't need it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "277bf13e-4c50-49fc-92f6-3d4cd679aec0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world_ddf.map_partitions(my_fn).compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c348112-6efa-4188-bcb7-96bff6ababaa",
   "metadata": {},
   "source": [
    "The task graph is the same as before, even though we return only a single value per partition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddb28ac8-c261-4042-98c4-cbee12a0a5a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf.map_partitions(my_fn).visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f61e4a0c-a697-4142-93e7-fa05f870ce99",
   "metadata": {},
   "source": [
    "If you want to assign the result as a new column of your dataframe, you need to ensure you return an array or a Series of the correct length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "940bdc65-cec9-416b-aa69-4e08090eb654",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_hull_area(gdf):\n",
    "    \"\"\"Get area of each convex hull and return pandas.Series\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    gdf : GeoDataFrame\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    pandas.Series\n",
    "    \"\"\"\n",
    "    \n",
    "    hulls = gdf.convex_hull\n",
    "    return hulls.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ab8b8f5-38e9-41c2-be31-08f660d8f7f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf['hull_area'] = world_ddf.map_partitions(get_hull_area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2f47751-0aad-4b9c-b4d0-3f9725d1726d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world_ddf.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0776eb2b-bbf3-4d0b-b11c-901406603a14",
   "metadata": {},
   "source": [
    "The task graph now includes the `assign` step, assigning the new column to the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "821b0352-d39f-4b2b-bab5-62634f48f925",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world_ddf.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3de6194f-824a-45fa-9bc4-45d2553c8d74",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Specifying meta-data\n",
    "\n",
    "To build a task graph, Dask doesn't need to see the data. But it needs to understand their general structure and dtypes. With simple `map_partitions` cases, you don't need to worry about that as Dask is often able to figure that out itself. But sometimes it struggles. \n",
    "\n",
    "Let's look a bit under the hood here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea6b4a4c-ef19-473d-a762-fa33b4396700",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf = dask_geopandas.read_parquet(\"data/world_4326/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c66c4e7-f754-40ff-8a6e-0d919a79d7e6",
   "metadata": {},
   "source": [
    "Each object contains its `_meta` data. For GeoDataFrame, that is usually an empty frame with columns and their dtypes set. Like this one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8270c25-c5d9-4d48-89af-91228629c07e",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf._meta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3ba021a-e725-490c-9f7e-7a28d7d614f4",
   "metadata": {},
   "source": [
    "You can see that all dtypes are set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9c5ab04-460a-4d75-9990-5d74aec3977b",
   "metadata": {},
   "outputs": [],
   "source": [
    "world_ddf._meta.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a4d36af-6f17-42cb-8d2d-7bd38c39d4e2",
   "metadata": {},
   "source": [
    "Now, we can try to implement our own version of `dissolve` that works well when all data fit in memory, to make the exmaple a bit more complicated. We need two steps:\n",
    "\n",
    "1. Shuffle the data into partitions based on the `continent` column. That ensures that all observations from the same continent are within a single partition.\n",
    "2. Use `map_partitions` to apply `dissolve` from geopandas.\n",
    "\n",
    "We can use the `shuffle` method to do the first step. Note that this is also a form of spatial partitioning and it may be useful to follow the attribute if you have one and the final partitions are of a roughly the same size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97a5a994-c2f9-480f-bd9d-787c7ee8753a",
   "metadata": {},
   "outputs": [],
   "source": [
    "shuffled = world_ddf.shuffle(\n",
    "    \"continent\", npartitions=7, ignore_index=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c65ccd1e-1529-476f-8695-6799508fe10b",
   "metadata": {},
   "outputs": [],
   "source": [
    "shuffled.calculate_spatial_partitions()\n",
    "shuffled.spatial_partitions.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78068a5a-83d5-4658-abb7-d3c4eca01d2e",
   "metadata": {},
   "source": [
    "Sometimes the inference of the meta DataFrame just fails, mostly because it is empty. If that happens, you can manually specify the meta data data frame and pass it to the `map_partitions`. Even if it does not fail, like in this case, it is often better to directly pass it as it can be chaper and you avoid potential issues that may come as a result of a wrong inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4071f1a6-6211-47b5-8f6d-02341242c30a",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta = world_ddf._meta.dissolve(by=\"continent\", as_index=False)\n",
    "meta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f11e880-898a-4255-bce1-af3eb815af42",
   "metadata": {},
   "source": [
    "With the `meta` defined, we can take the `geopandas.GeoDataFrame.dissolve` method and pass it to `map_partitions`. All `**kwargs` are passed as attributes to the method/function you are applying."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39e38f18-c95e-4a75-89da-519eaf9e0295",
   "metadata": {},
   "outputs": [],
   "source": [
    "dissolved = shuffled.map_partitions(\n",
    "    geopandas.GeoDataFrame.dissolve, by=\"continent\", as_index=False, meta=meta\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e27c056-5169-46c9-8498-0bad90f93894",
   "metadata": {},
   "outputs": [],
   "source": [
    "dissolved.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf2125b0-b23f-4e63-a623-b6e6fa0508d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "dissolved.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6e5baa6-3368-40b3-a350-cfbfe62bd556",
   "metadata": {},
   "source": [
    "When you are done, you can shut down the Dask client. If you want to do the exercises below, do not do that yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "139b7bd1-d45f-4a1f-ae16-4e3103979532",
   "metadata": {},
   "outputs": [],
   "source": [
    "client.shutdown()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6174c47-1287-4819-9ea2-a1cda6a6b298",
   "metadata": {},
   "source": [
    "## Limits and caveats\n",
    "\n",
    "Truth to be told, we are now playing with the version 0.1 of dask-geopandas and not everything is as polished as we would like it to be. So there are some things which are not yet fully supported.\n",
    "\n",
    "- **Overlapping computation** - With `dask.dataframe` and `dask.array` you can use `map_overlap` to do some overlapping computations for which you need observations from neighbouring partitions. With dask-geopandas, we would need this overlap to be spatial and that is not yet supported. That means that whatever depends on topology or similar operation is currently not very easy to parallelise.\n",
    "- **Spatial indexing** - While you can use spatial indexing over spatial partitions and then within individual partitions as we use it under the hood in `sjoin`, it requires a bit of low-level dask code to make it correctly run. We hope to make that easier at some point in the future. We also want to expand the use of the spatial partitioning information to more methods (currently only `sjoin` makes use of it).\n",
    "- **Memory management** - Even though Dask can work out-of-core and you may seem dask-geopadnas behaving that way sometimes, we still have some unresolved memory issues due to geometries being stored as C objects, hence their actual size is not directly visible to Dask.\n",
    "\n",
    "There are also the same, or at least very similar, rules when not to use dask-geopandas as they apply to vanilla dask.dataframe (from [Dask documentation](https://docs.dask.org/en/stable/dataframe-best-practices.html)):\n",
    "\n",
    "- For data that fits into RAM, geopandas can often be faster and easier to use than Dask. While “Big Data” tools can be exciting, they are almost always worse than normal data tools while those remain appropriate. But for embarrasinlgy parallel computation, it will often bring speedup with minimal overhead.\n",
    "- Similar to above, even if you have a large dataset there may be a point in your computation where you’ve reduced things to a more manageable level. You may want to switch to (geo)pandas at this point.\n",
    "\n",
    "```\n",
    "df = dd.read_parquet('my-giant-file.parquet')\n",
    "df = df[df.name == 'Alice']              # Select a subsection\n",
    "result = df.groupby('id').value.mean()   # Reduce to a smaller size\n",
    "result = result.compute()                # Convert to pandas dataframe\n",
    "result...                                # Continue working with pandas\n",
    "```\n",
    "\n",
    "- Usual pandas and GeoPandas performance tips like avoiding apply, using vectorized operations, using categoricals, etc., all apply equally to Dask DataFrame and dask-geopandas.\n",
    "\n",
    "See more best practices in the [Dask documentation](https://docs.dask.org/en/stable/dataframe-best-practices.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40524b9e-3b01-488b-83b9-b2c8aa9cf663",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "Use the `data/airports_csv` folder and try to do the following using Dask:\n",
    "\n",
    "- Read the contents as a `dask.dataframe`\n",
    "- Create a valid `dask_geopandas.GeoDataFrame`\n",
    "- Calculate and explore spatial partitions. If you think there's a need to spatially shuffle the data, do so.\n",
    "    - Try comparing different sorting methods (check the docs!). Which one is the best an why?\n",
    "- How many airports are there per continent? And how many per country?\n",
    "- Are there any points not falling onto ground? How many?\n",
    "- Where would be the ideal single airport in each country if you had to build only one?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geopython_tutorial",
   "language": "python",
   "name": "geopython_tutorial"
  },
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}