deforestation-sample.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "#import geopandas as gpd\n",
    "import requests\n",
    "import rasterio as rio\n",
    "from rasterio.plot import show\n",
    "from rasterio.warp import transform\n",
    "from rasterio.transform import (xy, rowcol)\n",
    "from rasterio.windows import Window\n",
    "import itertools as it\n",
    "import rioxarray as rx\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Rectangle\n",
    "import numpy as np\n",
    "from pathlib import Path\n",
    "import time\n",
    "import timeit\n",
    "import os.path\n",
    "import math\n",
    "from string import Template\n",
    "\n",
    "from typing import Tuple, Optional, List\n",
    "\n",
    "from leaf.deforestation import (\n",
    "    cache_earthenginepartners_hansen,\n",
    "    earthenginepartners_hansen)\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/assets_for_deforestation.csv', sep='\\t')\n",
    "df.country.value_counts().head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# spanning the range 180W-180E and 80N-60S\n",
    "# 80N, 70N, ... 10N, 00N, 10S, 20S, ... 50S - 180W\n",
    "# 80N - 180W, 170W, ... 010W, 000E, 010E.... 170E\n",
    "\n",
    "# https://worldpopulationreview.com/country-rankings/deforestation-rates-by-country\n",
    "# https://gist.github.com/graydon/11198540\n",
    "# long lat bottom left, long lat top right\n",
    "country_bounding_boxes = {\n",
    "    'BR': ('Brazil', (-73.9872354804, -33.7683777809, -34.7299934555, 5.24448639569)),\n",
    "    'ID': ('Indonesia', (95.2930261576, -10.3599874813, 141.03385176, 5.47982086834)),\n",
    "    'CD': ('Congo (Kinshasa)', (12.1823368669, -13.2572266578, 31.1741492042, 5.25608775474)),\n",
    "    'AO': ('Angola', (11.6400960629, -17.9306364885, 24.0799052263, -4.43802336998)),\n",
    "    #Tazmania\n",
    "    'MM': ('Myanmar', (92.3032344909, 9.93295990645, 101.180005324, 28.335945136)),\n",
    "    #'PY': ('Paraguay', (-62.6850571357, -27.5484990374, -54.2929595608, -19.3427466773)),\n",
    "    #'BO': ('Bolivia', (-69.5904237535, -22.8729187965, -57.4983711412, -9.76198780685)),\n",
    "    'MZ': ('Mozambique', (30.1794812355, -26.7421916643, 40.7754752948, -10.3170960425)),\n",
    "    'AR': ('Argentina', (-73.4154357571, -55.25, -53.628348965, -21.8323104794)),\n",
    "    'AU': ('Australia', (113.338953078, -43.6345972634, 153.569469029, -10.6681857235)),\n",
    "    'VN': ('Vietnam', (102.170435826, 8.59975962975, 109.33526981, 23.3520633001)),\n",
    "    'PE': ('Peru', (-81.4109425524, -18.3479753557, -68.6650797187, -0.0572054988649)),\n",
    "    'PG': ('Papua New Guinea', (141.000210403, -10.6524760881, 156.019965448, -2.50000212973)),\n",
    "    'PH': ('Philippines', (117.17427453, 5.58100332277, 126.537423944, 18.5052273625)),\n",
    "    'CO': ('Colombia', (-78.9909352282, -4.29818694419, -66.8763258531, 12.4373031682)),\n",
    "    # Note bene: Low in ranking but high in assets...\n",
    "    'CN': ('China', (73.6753792663, 18.197700914, 135.026311477, 53.4588044297)),\n",
    "    'US': ('United States', (-171.791110603, 18.91619, -66.96466, 71.3577635769)),\n",
    "}\n",
    "\n",
    "country_bounding_boxes_adds = {\n",
    "}\n",
    "\n",
    "def ranges(bboxes, country_code: str, step = 10) -> Tuple[range, range]:\n",
    "    country = bboxes[country_code]\n",
    "    temp = tuple(map(lambda x: math.ceil(x), country[1]))\n",
    "    latitudes = range(temp[1], temp[3], step)\n",
    "    longitudes = range(temp[0], temp[2], step)\n",
    "    ranges = (latitudes, longitudes)\n",
    "    return ranges\n",
    "\n",
    "def rectangle(bboxes, country_code: str) -> Rectangle:\n",
    "    country = bboxes[country_code]\n",
    "    temp = tuple(map(lambda x: math.ceil(x), country[1]))\n",
    "    rectangle = Rectangle(xy=(temp[0], temp[1]),\n",
    "                          width=abs(temp[2] - temp[0]), height=abs(temp[3] - temp[1]),\n",
    "                          linewidth=1, color='red', fill=False)\n",
    "    return rectangle\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a map of areas of interest\n",
    "fig = plt.figure(figsize=(15, 10))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "\n",
    "#ax.set_extent([-73, -35, -33, 6])  # Define the extent of the map (lon_min, lon_max, lat_min, lat_max)\n",
    "\n",
    "# Plot assets\n",
    "\n",
    "for sector_main in df.sector_main.unique():\n",
    "    foo = df[df.sector_main == sector_main]\n",
    "    ax.scatter(foo['longitude'], foo['latitude'],\n",
    "        marker='o', \n",
    "        s=2, alpha = 0.5, label=sector_main,\n",
    "        transform=ccrs.PlateCarree())\n",
    "\n",
    "codes = [code for code, _ in country_bounding_boxes.items()]\n",
    "for code in codes:\n",
    "    ax.add_patch(rectangle(country_bounding_boxes, code))\n",
    "\n",
    "# Add map features\n",
    "ax.coastlines()\n",
    "\n",
    "# Add country borders\n",
    "ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=1)\n",
    "\n",
    "# Add title and labels\n",
    "plt.title('Assets from GEM data')\n",
    "plt.xlabel('Longitude')\n",
    "plt.ylabel('Latitude')\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "# save the plot\n",
    "fig.savefig('graphs/map_assets_all.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# L807043,  -8.2664, -42.1428\n",
    "\n",
    "asset_ranges = ([-9, -8, 10], [-43, -42, 10])\n",
    "earthenginepartners_hansen('data/assets_for_deforestation.csv', '\\t', \n",
    "                           asset_ranges[0], asset_ranges[1],\n",
    "                           'data/L807043_with_deforestation.csv', verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nota bene: Robert Norris - this will re-create assets_with_deforestation.csv from \n",
    "# assets_for_deforestation.csv, using whatever countries have been included in \n",
    "# country_bounding_boxes. Any additionally required GeoTIFF files from the Hansen \n",
    "# dataset will be downloaded and cached. It uses the default offset of 16 to inspect \n",
    "# an region 990m x 990m centred on each asset.\n",
    "\n",
    "brazil = ranges(country_bounding_boxes, 'BR')\n",
    "earthenginepartners_hansen('data/assets_for_deforestation.csv', '\\t', \n",
    "                           brazil[0], brazil[1],\n",
    "                           'data/assets_with_deforestation.csv', verbose=False)\n",
    "\n",
    "remaining_country_codes = [code for code in country_bounding_boxes.keys() if code != 'BR']\n",
    "for code in remaining_country_codes:\n",
    "    country_ranges = ranges(country_bounding_boxes, code)\n",
    "    earthenginepartners_hansen('data/assets_with_deforestation.csv', '\\t', \n",
    "                               country_ranges[0], country_ranges[1],\n",
    "                               'data/assets_with_deforestation.csv', verbose=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "country_codes = [code for code in country_bounding_boxes_adds.keys()]\n",
    "for code in country_codes:\n",
    "    country_ranges = ranges(country_bounding_boxes, code)\n",
    "    earthenginepartners_hansen('data/assets_with_deforestation.csv', '\\t', \n",
    "                               country_ranges[0], country_ranges[1],\n",
    "                               'data/assets_with_deforestation.csv', verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nota bene: Robert Norris - this will re-create assets_with_deforestation_n.csv \n",
    "# where n is the dimension in metres, from assets_for_deforestation.csv, using whatever countries have \n",
    "# been included in country_bounding_boxes. Any additionally required GeoTIFF files from the Hansen \n",
    "# dataset will be downloaded and cached. Increasing the size may ignore more assets at the borders\n",
    "# of each GeoTIFF but more of the surrounding area will be considered per asset.\n",
    "\n",
    "for offset in [32, 48, 80, 128]:\n",
    "\n",
    "    dimension = (offset * 2 + 1) * 30 # each pixel is ~30m x 30m\n",
    "    GEMFile = f'data/assets_with_deforestation_{dimension}.csv'\n",
    "\n",
    "    brazil = ranges(country_bounding_boxes, 'BR')\n",
    "    earthenginepartners_hansen('data/assets_for_deforestation.csv', '\\t', \n",
    "                            brazil[0], brazil[1],\n",
    "                            GEMFile, offset, verbose=False)\n",
    "\n",
    "    remaining_country_codes = [code for code in country_bounding_boxes.keys() if code != 'BR']\n",
    "    for code in remaining_country_codes:\n",
    "        country_ranges = ranges(country_bounding_boxes, code)\n",
    "        earthenginepartners_hansen(GEMFile, '\\t', \n",
    "                                country_ranges[0], country_ranges[1],\n",
    "                                GEMFile, offset, verbose=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import missingno as msno\n",
    "\n",
    "df = pd.read_csv('data/assets_with_deforestation.csv', sep='\\t')\n",
    "msno.bar(df)\n",
    "# len(f\"Number of valid observations: {df[df['2022'] >= 0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries = [country[0] for code, country in country_bounding_boxes.items()]\n",
    "for country in countries:\n",
    "    value_counts = df[df.treecover2000.notna() & (df.country == country)].sector_main.value_counts()\n",
    "    print(f'{country}\\n{value_counts}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lossyears = [str(lossyear) for lossyear in range(2001, 2023)]\n",
    "countries = [country[0] for code, country in country_bounding_boxes.items()]\n",
    "with pd.option_context('expand_frame_repr', False):\n",
    "    for country in countries:\n",
    "        print(country)\n",
    "        print(df[(df[lossyears] > 0).any(axis=1) & (df.country == country)].head(2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[(df[lossyears] > 0).any(axis=1)].country.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msno.heatmap(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/regression_sample.csv', sep='\\t')\n",
    "#df.country.value_counts().head(10)\n",
    "#df[df.country == 'Brazil'].head()\n"
   ]
  }
 ],
 "metadata": {
  "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}