1.AmazonasLimits.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "o6EiLjE7tL3k"
   },
   "outputs": [],
   "source": [
    "import urbanpy as up\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import contextily as ctx\n",
    "import rioxarray\n",
    "import os\n",
    "from tqdm.auto import tqdm\n",
    "from shapely.geometry import Polygon, MultiPolygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas = gpd.read_file(\"inputs/Amazonas\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas.crs.to_string()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas = amazonas.to_crs(\"EPSG:4326\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas.crs.to_string()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 447
    },
    "id": "v07RB4VYlXkZ",
    "outputId": "03c18114-81fb-4ad6-93b4-5dd6293d4f22"
   },
   "outputs": [],
   "source": [
    "ax = amazonas.plot(facecolor=\"none\", edgecolor=\"black\", linewidth=0.7)\n",
    "ctx.add_basemap(ax, crs=amazonas.crs.to_string())\n",
    "ax.set_title(\"Amazonas Limits\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas[\"geometry\"].explode()[0].area.sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_gdf = gpd.GeoDataFrame(\n",
    "    geometry=[amazonas[\"geometry\"].explode()[0][122]], crs=amazonas.crs\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = amazonas_gdf.plot(facecolor=\"none\", edgecolor=\"black\", linewidth=0.7)\n",
    "ctx.add_basemap(ax, crs=amazonas_gdf.crs.to_string())\n",
    "ax.set_title(\"Amazonas Limits\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_3D_2D(geometry):\n",
    "    \"\"\"\n",
    "    Takes a GeoSeries of 3D Multi/Polygons (has_z) and returns a list of 2D Multi/Polygons\n",
    "    \"\"\"\n",
    "    new_geo = []\n",
    "    for p in geometry:\n",
    "        if p.has_z:\n",
    "            if p.geom_type == \"Polygon\":\n",
    "                lines = [xy[:2] for xy in list(p.exterior.coords)]\n",
    "                new_p = Polygon(lines)\n",
    "                new_geo.append(new_p)\n",
    "            elif p.geom_type == \"MultiPolygon\":\n",
    "                new_multi_p = []\n",
    "                for ap in p:\n",
    "                    lines = [xy[:2] for xy in list(ap.exterior.coords)]\n",
    "                    new_p = Polygon(lines)\n",
    "                    new_multi_p.append(new_p)\n",
    "                new_geo.append(MultiPolygon(new_multi_p))\n",
    "    return new_geo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_gdf_2d = convert_3D_2D(amazonas_gdf.geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_clean = gpd.GeoDataFrame(geometry=amazonas_gdf_2d, crs=amazonas.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_clean.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_clean.to_parquet(\"outputs/amazonas_clean.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# amazonas_hexs_1 = up.geom.gen_hexagons(resolution=1, city=amazonas_clean)\n",
    "# amazonas_hexs_1.drop(\"geometry\").to_csv(\"outputs/amazonas_hexs_1.csv\", index=False)\\\n",
    "\n",
    "amazonas_hexs_2 = up.geom.gen_hexagons(resolution=2, city=amazonas_clean)\n",
    "amazonas_hexs_2.drop(\"geometry\", axis=1).to_csv(\n",
    "    \"outputs/amazonas_hexs_2.csv\", index=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_3 = up.geom.gen_hexagons(resolution=3, city=amazonas_clean)\n",
    "amazonas_hexs_3.drop(\"geometry\", axis=1).to_csv(\n",
    "    \"outputs/amazonas_hexs_3.csv\", index=False\n",
    ")\n",
    "\n",
    "amazonas_hexs_4 = up.geom.gen_hexagons(resolution=4, city=amazonas_clean)\n",
    "amazonas_hexs_4.drop(\"geometry\", axis=1).to_csv(\n",
    "    \"outputs/amazonas_hexs_4.csv\", index=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "\n",
    "amazonas_gdf.plot(facecolor=\"none\", edgecolor=\"black\", linewidth=0.7, ax=ax)\n",
    "amazonas_hexs_4.plot(ax=ax, facecolor=\"none\", edgecolor=\"red\", linewidth=0.3)\n",
    "ctx.add_basemap(ax, crs=amazonas_gdf.crs.to_string())\n",
    "ax.set_axis_off()\n",
    "\n",
    "amazonas_gdf.plot(facecolor=\"none\", edgecolor=\"black\", linewidth=0.7, ax=ax1)\n",
    "amazonas_hexs_4.plot(ax=ax1, facecolor=\"none\", edgecolor=\"red\", linewidth=0.3)\n",
    "ctx.add_basemap(ax1, crs=amazonas_gdf.crs.to_string())\n",
    "ax1.set_axis_off()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\n",
    "    \"outputs/amazonas_hexs_4_2maps.png\", dpi=300, bbox_inches=\"tight\", format=\"png\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = amazonas_gdf.plot(\n",
    "    facecolor=\"none\", edgecolor=\"black\", linewidth=0.7, figsize=(10, 10)\n",
    ")\n",
    "amazonas_hexs_4.plot(ax=ax, facecolor=\"none\", edgecolor=\"red\", linewidth=0.3)\n",
    "ctx.add_basemap(ax, crs=amazonas_gdf.crs.to_string())\n",
    "ax.set_title(\"Amazonas Hexs Res 4\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "plt.savefig(\"outputs/amazonas_hexs_4.svg\", dpi=300, bbox_inches=\"tight\", format=\"svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_5 = up.geom.gen_hexagons(resolution=5, city=amazonas_clean)\n",
    "amazonas_hexs_5.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_5.to_parquet(\"outputs/amazonas_hexs_5.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_6 = up.geom.gen_hexagons(resolution=6, city=amazonas_clean)\n",
    "amazonas_hexs_6.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_6.to_parquet(\"outputs/amazonas_hexs_6.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_7 = up.geom.gen_hexagons(resolution=7, city=amazonas_clean)\n",
    "amazonas_hexs_7.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amazonas_hexs_7.to_parquet(\"outputs/amazonas_hexs_7.parquet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Takes too much time\n",
    "# amazonas_hexs_8 = up.geom.gen_hexagons(resolution=5, city=amazonas_clean)\n",
    "# amazonas_hexs_8.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country_boundaries = gpd.read_file(\"inputs/Cartographic Boundary Files/LAC/level 0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country_boundaries.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country_boundaries.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country_boundaries[\"ADM0_PCODE\"].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amzn_countries = country_boundaries[\n",
    "    country_boundaries[\"ADM0_PCODE\"].isin(\n",
    "        [\"BO\", \"BR\", \"CO\", \"EC\", \"GY\", \"PE\", \"SR\", \"VE\"]\n",
    "    )\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = amazonas_gdf.plot(facecolor=\"none\", edgecolor=\"black\", linewidth=0.7)\n",
    "amzn_countries.plot(\"ADM0_PCODE\", ax=ax, alpha=0.5)\n",
    "ctx.add_basemap(ax, crs=amazonas_gdf.crs.to_string())\n",
    "ax.set_title(\"Amazon Countries\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WoldPop - Age and sex structures\n",
    "\n",
    "- Resolution: 100m^2\n",
    "- Year: 2020\n",
    "- Classes: 5-year age groups + <1 year\n",
    "- Version: Constrained\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o5mA7V9-nK8-"
   },
   "source": [
    "# Raster based version of population data download\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "worldpop_data = \"inputs/WorldPop\"\n",
    "countries = [\"PER\", \"COL\", \"GUY\", \"SUR\", \"VEN\", \"BOL\", \"ECU\"]  # \"BRA\"\n",
    "age_groups = [5, 10, 15]\n",
    "genders = [\"m\", \"f\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_rio_hex(\n",
    "    hexs: gpd.GeoDataFrame,\n",
    "    clip_geometries: gpd.GeoDataFrame,\n",
    "    rio_filename: str,\n",
    "    data_name: str,\n",
    "    agg: str,\n",
    "    band: int = 1,\n",
    ") -> gpd.GeoDataFrame:\n",
    "    \"\"\"\n",
    "    Merge raster data with hexagons.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    hexs : gpd.GeoDataFrame\n",
    "        GeoDataFrame with hexagons.\n",
    "    rio_filename : str\n",
    "        Filename of the raster data.\n",
    "    agg : str\n",
    "        Aggregation method to use.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gpd.GeoDataFrame\n",
    "        GeoDataFrame with the merged data.\n",
    "    \"\"\"\n",
    "    rds = rioxarray.open_rasterio(rio_filename)\n",
    "    rds.rio.set_crs(hexs.crs.to_string())\n",
    "    clipped = rds.rio.clip(clip_geometries, rds.rio.crs)\n",
    "    rio_data = clipped.sel(band=band).drop_vars(\"band\")\n",
    "    rio_data.name = data_name\n",
    "    rio_data_df = rio_data.to_pandas().unstack().reset_index()\n",
    "    rio_data_df.columns = [\"x\", \"y\", data_name]\n",
    "    rio_data_gdf = gpd.GeoDataFrame(\n",
    "        rio_data_df,\n",
    "        geometry=gpd.points_from_xy(rio_data_df.x, rio_data_df.y),\n",
    "        crs=rio_data.rio.crs.to_string(),\n",
    "    )\n",
    "    rio_data_gdf = rio_data_gdf.dropna(\n",
    "        subset=[data_name]\n",
    "    )  # Drop NaN values for faster processing\n",
    "    agg_dict = {data_name: agg}\n",
    "    return up.geom.merge_shape_hex(hexs=hexs, shape=rio_data_gdf, agg=agg_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ecuador\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecuador_bounds = amzn_countries[amzn_countries[\"ADM0_PCODE\"] == \"EC\"]\n",
    "ecuador_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = ecuador_bounds.explode().loc[10]\n",
    "temp[\"area\"] = temp.area\n",
    "temp.plot(\"area\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter continental ecuador\n",
    "ecuador_gdf = temp.cx[-82.5:-75, -5:-0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecuador_hexs = up.geom.gen_hexagons(resolution=7, city=ecuador_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "xHuEfEc6lxKh",
    "outputId": "f48ad117-2eac-47fe-f92a-bf0c913ec620"
   },
   "outputs": [],
   "source": [
    "print(\"Number of hexagons:\", ecuador_hexs.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 447
    },
    "id": "p2jWd086ltqn",
    "outputId": "b8aba2cb-561d-4ddb-eb24-75dfeb01f44b"
   },
   "outputs": [],
   "source": [
    "ax = ecuador_hexs.plot(facecolor=\"none\", edgecolor=\"black\", linewidth=0.25)\n",
    "ctx.add_basemap(ax, crs=ecuador_hexs.crs.to_string())\n",
    "ax.set_axis_off()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "progress_bar = tqdm(total=len(genders) * len(age_groups))\n",
    "country_string = \"ecu\"\n",
    "year = \"2020\"\n",
    "for gender in genders:\n",
    "    for age in age_groups:\n",
    "        filename = f\"{worldpop_data}/{country_string}/{country_string}_{gender}_{age}_{year}.tif\"\n",
    "        assert os.path.exists(filename), f\"File {filename} does not exist\"\n",
    "        ecuador_hexs = merge_rio_hex(\n",
    "            hexs=ecuador_hexs,\n",
    "            clip_geometries=ecuador_gdf.geometry,\n",
    "            rio_filename=filename,\n",
    "            data_name=f\"{country_string}_{gender}_{age}_{year}\",\n",
    "            agg=\"sum\",\n",
    "            band=1,\n",
    "        )\n",
    "        progress_bar.update(1)\n",
    "progress_bar.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}