matreed_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MaTREEd data preprocessing notebook\n",
    "This notebook has been developed for preparing the data exposed by the MaTREEd web application (https://gisdevio.github.io/MaTREEd). More info about the project and the data preprocessing operations can be found at https://github.com/GISdevio/MaTREEd/blob/main/project-description.md."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "0) Import necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Import raw data \n",
    "\n",
    "Original data links: \n",
    "- trees: https://challenge.greemta.eu/dataset/trees\n",
    "- neighborhoods: https://challenge.greemta.eu/dataset/administrativeunits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create dataframes from raw data (trees and neighborhoods)\n",
    "trees_csv = pd.read_csv(r'raw_data/trees_madrid.zip', compression='zip', sep=',', quotechar='\"', error_bad_lines=False, encoding='utf8')\n",
    "neigh = gpd.read_file(r'raw_data/neighborhoods_madrid.geojson')\n",
    "\n",
    "# Make trees_csv a GeoDataframe\n",
    "trees_geo = gpd.GeoDataFrame(trees_csv, geometry=gpd.points_from_xy(trees_csv['longitude'], trees_csv['latitude']), crs='EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check raw tree data\n",
    "trees_csv.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check raw neighborhoods data\n",
    "neigh.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Preliminary operations on raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add CODBARRIO to trees_geo\n",
    "trees_join = gpd.sjoin(trees_geo, neigh, how='left', op='within')\n",
    "\n",
    "# Remove some columns to trees_geo and create the trees geodataframe for processing\n",
    "trees = trees_join[['objectid',\n",
    " 'specie_name',\n",
    " 'senescence',\n",
    " 'crown_diameter',\n",
    " 'height',\n",
    " 'trunk_girth',\n",
    " 'geometry',\n",
    " 'CODBARRIO']]\n",
    "\n",
    "# Add area to neighbourhoods [Km2]\n",
    "neigh['area_km2']= neigh['geometry'].to_crs({'init': 'epsg:32630'}).map(lambda p: p.area / 10**6)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Data Cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove null values\n",
    "\n",
    "# Remove trees with NaN values in columns senescence, crown_diameter, height and trunk_girth\n",
    "df_nan = trees.dropna(subset = ['senescence', 'crown_diameter', 'height', 'trunk_girth'])\n",
    "\n",
    "# Remove possible outliers (based on manual data exploration)\n",
    "# Keep only remaining trees with values in column crown_diameter between the 1st and 99.99th percentile of the column\n",
    "df_dc = df_nan.drop(df_nan[np.logical_or((df_nan.crown_diameter < df_nan.crown_diameter.quantile(0.01)),(df_nan.crown_diameter > df_nan.crown_diameter.quantile(0.9999)))].index)\n",
    "# Keep only remaining trees with values in column height between the 1st and 99.999th percentile of the column \n",
    "df_a = df_dc.drop(df_dc[np.logical_or((df_dc.height < df_dc.height.quantile(0.01)),(df_dc.height > df_dc.height.quantile(0.99999)))].index)\n",
    "# Keep only remaining trees with values in column trunk_girth between the 1st and 99.99th percentile of the column\n",
    "df_p = df_a.drop(df_a[np.logical_or((df_a.trunk_girth < df_a.trunk_girth.quantile(0.01)),(df_a.trunk_girth > df_a.trunk_girth.quantile(0.9999)))].index)\n",
    "\n",
    "# save the clean trees dataframe\n",
    "trees_clean = df_p\n",
    "\n",
    "# flush no more useful dataframes from memory\n",
    "del[df_nan, df_dc, df_a, df_p, trees_join, trees_csv,trees_geo]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Compute trees statistics for neighbours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add tree counts and percentages to neighborhoods (total and after cleaning)\n",
    "count_tot = trees.CODBARRIO.value_counts().to_frame()\n",
    "count_tot.rename(columns={'CODBARRIO':'count_tot'}, inplace=True)\n",
    "\n",
    "count_clean = trees_clean.CODBARRIO.value_counts().to_frame()\n",
    "count_clean.rename(columns={'CODBARRIO':'count_clean'}, inplace=True)\n",
    "\n",
    "neigh['tree_tot']= neigh.set_index('CODBARRIO').join(count_tot).set_index(neigh.index)['count_tot']\n",
    "neigh['tree_clean']= neigh.set_index('CODBARRIO').join(count_clean).set_index(neigh.index)['count_clean']\n",
    "\n",
    "#  Add the percentage of complete tree records out of the total available (used to compute stats\n",
    "neigh['tree_valid']= round((neigh['tree_clean'] / neigh['tree_tot'])*100, 2)\n",
    "\n",
    "# Add perennifolio_count and caducifolio_count to neighbours (after cleaning)\n",
    "trees_perennifolio = trees_clean[trees_clean[\"senescence\"] == 'PERENNIFOLIO']\n",
    "count_perennifolio = trees_perennifolio.CODBARRIO.value_counts().to_frame()\n",
    "count_perennifolio.rename(columns={'CODBARRIO':'count_perennifolio'}, inplace=True)\n",
    "\n",
    "trees_caducifolio = trees_clean[trees_clean[\"senescence\"] == 'CADUCIFOLIO']\n",
    "count_caducifolio = trees_caducifolio.CODBARRIO.value_counts().to_frame()\n",
    "count_caducifolio.rename(columns={'CODBARRIO':'count_caducifolio'}, inplace=True)\n",
    "\n",
    "neigh['count_perennifolio']= neigh.set_index('CODBARRIO').join(count_perennifolio).set_index(neigh.index)['count_perennifolio']\n",
    "neigh['count_caducifolio']= neigh.set_index('CODBARRIO').join(count_caducifolio).set_index(neigh.index)['count_caducifolio']\n",
    "neigh['perc_perennifolio']= round((neigh['count_perennifolio'] / neigh['tree_clean'])*100, 2)\n",
    "neigh['perc_caducifolio']= round((neigh['count_caducifolio'] / neigh['tree_clean'])*100, 2)\n",
    "\n",
    "# Add avg_trunk_girth, avg_crown_diameter, sum_diammetro_copa, avg_height, tree_cover to neighbours (after cleaning)\n",
    "avg_crown_diameter=trees_clean.groupby('CODBARRIO')['crown_diameter'].mean().to_frame()\n",
    "avg_crown_diameter.rename(columns={'crown_diameter':'avg_crown_diameter'}, inplace=True)\n",
    "\n",
    "trees_clean['AREA_COPA']=pow(trees_clean['crown_diameter'],2)*np.pi/4\n",
    "sum_area_copa=trees_clean.groupby('CODBARRIO')['AREA_COPA'].sum().to_frame()\n",
    "sum_area_copa.rename(columns={'AREA_COPA':'sum_area_copa'}, inplace=True)\n",
    "\n",
    "avg_height=trees_clean.groupby('CODBARRIO')['height'].mean().to_frame()\n",
    "avg_height.rename(columns={'height':'avg_height'}, inplace=True)\n",
    "\n",
    "avg_trunk_girth=trees_clean.groupby('CODBARRIO')['trunk_girth'].mean().to_frame()\n",
    "avg_trunk_girth.rename(columns={'trunk_girth':'avg_trunk_girth'}, inplace=True)\n",
    "\n",
    "neigh['avg_crown_diameter']= neigh.set_index('CODBARRIO').join(avg_crown_diameter).set_index(neigh.index)['avg_crown_diameter']\n",
    "neigh['sum_area_copa']= neigh.set_index('CODBARRIO').join(sum_area_copa).set_index(neigh.index)['sum_area_copa']\n",
    "neigh[\"tree_cover\"]=round(neigh[\"sum_area_copa\"]/(neigh[\"area_km2\"]*10**6)*100,2)\n",
    "neigh['avg_height']= neigh.set_index('CODBARRIO').join(avg_height).set_index(neigh.index)['avg_height']\n",
    "neigh['avg_trunk_girth']= neigh.set_index('CODBARRIO').join(avg_trunk_girth).set_index(neigh.index)['avg_trunk_girth']\n",
    "\n",
    "\n",
    "# Add name and percentage of the three most abundant species in each neighborhood\n",
    "# Create a new DataFrame containing name and count of the three most abundant species for each neighborhood\n",
    "data=[]\n",
    "for index,row in neigh.iterrows():\n",
    "    esp_trees = trees_clean[trees_clean[\"CODBARRIO\"] == row.CODBARRIO]\n",
    "    list_trees = esp_trees.specie_name.value_counts()[0:3].to_frame()\n",
    "    list_trees.rename(columns={'specie_name':'count_especie'}, inplace=True)\n",
    "    if len(list_trees)==1:\n",
    "        data.append({'CODBARRIO':row.CODBARRIO,'especie_1':list_trees.index[0], 'especie_1_count':list_trees.count_especie[0] , 'especie_2':'', 'especie_2_count':None, 'especie_3':'', 'especie_3_count':None})\n",
    "    elif len(list_trees)==2:\n",
    "        data.append({'CODBARRIO':row.CODBARRIO,'especie_1':list_trees.index[0], 'especie_1_count':list_trees.count_especie[0] , 'especie_2':list_trees.index[1], 'especie_2_count':list_trees.count_especie[1], 'especie_3':'', 'especie_3_count':None})\n",
    "    elif len(list_trees)==3:\n",
    "        data.append({'CODBARRIO':row.CODBARRIO,'especie_1':list_trees.index[0], 'especie_1_count':list_trees.count_especie[0] , 'especie_2':list_trees.index[1], 'especie_2_count':list_trees.count_especie[1], 'especie_3':list_trees.index[2], 'especie_3_count':list_trees.count_especie[2] })\n",
    "    else:\n",
    "        data.append({'CODBARRIO':row.CODBARRIO,'especie_1':'', 'especie_1_count':None, 'especie_2':'', 'especie_2_count':None, 'especie_3':'', 'especie_3_count':None})\n",
    "\n",
    "agg_list_tree = pd.DataFrame(data)\n",
    "\n",
    "\n",
    "# Merge the new DataFrame and the neighbor DataFrame according to CODBARRIO column\n",
    "neigh = pd.merge(neigh,agg_list_tree,on='CODBARRIO')\n",
    "del agg_list_tree\n",
    "\n",
    "# Compute the percentage of the species\n",
    "neigh[\"especie_1_perc\"]=round((neigh[\"especie_1_count\"]/neigh['tree_clean'])*100,2)\n",
    "neigh[\"especie_2_perc\"]=round((neigh[\"especie_2_count\"]/neigh['tree_clean'])*100,2)\n",
    "neigh[\"especie_3_perc\"]=round((neigh[\"especie_3_count\"]/neigh['tree_clean'])*100,2)\n",
    "neigh[\"other_perc\"]=100-neigh[\"especie_1_perc\"]-neigh[\"especie_2_perc\"]-neigh[\"especie_3_perc\"]\n",
    "\n",
    "# flush no more useful dataframes from memory\n",
    "del[avg_height, avg_crown_diameter, avg_trunk_girth, count_caducifolio, count_perennifolio, count_clean, count_tot,\n",
    "    esp_trees, sum_area_copa, trees_caducifolio, trees_perennifolio]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5) Compute carbon dioxide stock and sequestration rate for trees and neighbours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''COMPUTE CARBON DIOXIDE INDICATORS'''\n",
    "\n",
    "# function to return Carbon Dioxide stock [Kg CO2]\n",
    "def Cstock(df):\n",
    "    trees = df \n",
    "    c_stock_list = []\n",
    "    for row in trees.iterfeatures():\n",
    "        sen = row['properties']['senescence']\n",
    "        d = (row['properties']['trunk_girth']/np.pi)*100\n",
    "        if sen == \"PERENNIFOLIO\":\n",
    "            c_stock = (0.16155*(d)**2.310647)*0.5*3.67\n",
    "        else :\n",
    "            c_stock = (0.035702*(d)**2.580671)*0.5*3.67\n",
    "        c_stock_list.append(c_stock)\n",
    "    return c_stock_list \n",
    "\n",
    "# function to return Carbon Dioxide sequestration rate [Kg CO2/year]\n",
    "def Cseq(df): \n",
    "    trees = df\n",
    "    c_seq_list = []\n",
    "    for row in trees.iterfeatures():\n",
    "        sen = row['properties']['senescence']\n",
    "        c_stock = row['properties'][\"c_stock\"]\n",
    "        d = (row['properties']['trunk_girth']/np.pi)*100\n",
    "        d_t1 = (-0.5425 + 0.3189*np.log(d)) + d \n",
    "        if sen == \"PERENNIFOLIO\":\n",
    "            c_stock_t1 = (0.16155*(d_t1)**2.310647)*0.5*3.67\n",
    "        else :\n",
    "            c_stock_t1 = (0.035702*(d_t1)**2.580671)*0.5*3.67\n",
    "        c_seq = c_stock_t1 - c_stock\n",
    "        if c_seq < 0.0:\n",
    "            c_seq = 0.0\n",
    "            \n",
    "        c_seq_list.append(c_seq)\n",
    "    return c_seq_list \n",
    "\n",
    "# Apply Carbon functions\n",
    "trees_clean[\"c_stock\"] = Cstock(trees_clean) \n",
    "trees_clean[\"c_seq\"] = Cseq(trees_clean) \n",
    "\n",
    "# Add cumulative Carbon Dioxide stock and Carbon Dioxide sequestration rate to neighborhoods\n",
    "sum_c_stock=trees_clean.groupby('CODBARRIO')['c_stock'].sum().to_frame()\n",
    "sum_c_stock.rename(columns={'c_stock':'sum_c_stock'}, inplace=True)\n",
    "\n",
    "sum_c_seq=trees_clean.groupby('CODBARRIO')['c_seq'].sum().to_frame()\n",
    "sum_c_seq.rename(columns={'c_seq':'sum_c_seq'}, inplace=True)\n",
    "\n",
    "neigh['sum_c_stock']= neigh.set_index('CODBARRIO').join(sum_c_stock).set_index(neigh.index)['sum_c_stock']\n",
    "neigh['sum_c_seq']= neigh.set_index('CODBARRIO').join(sum_c_seq).set_index(neigh.index)['sum_c_seq']\n",
    "\n",
    "# normalize stats on neighborhood area\n",
    "neigh['tree_tot_n'] = neigh['tree_tot']/neigh[\"area_km2\"]\n",
    "neigh['sum_c_stock_n'] = neigh['sum_c_stock']/neigh[\"area_km2\"]\n",
    "neigh['sum_c_seq_n'] = neigh['sum_c_seq']/neigh[\"area_km2\"]\n",
    "\n",
    "# Set NaN values of columns of interest to 0\n",
    "cols=[\"tree_tot\",\"avg_trunk_girth\",\"avg_height\",\"avg_crown_diameter\",\"tree_cover\",\"perc_perennifolio\",\"perc_caducifolio\",\"especie_1_perc\",\"especie_2_perc\",\"especie_3_perc\",\"other_perc\",\"tree_valid\",\"sum_c_stock\",\"sum_c_seq\"]\n",
    "\n",
    "for c in cols:\n",
    "     neigh[c] = neigh[c].fillna(0)\n",
    "     \n",
    "# flush no more useful dataframes and variables from memory \n",
    "del[sum_c_seq, sum_c_stock, row, c, cols, data, index, list_trees]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check final tree data\n",
    "trees_clean.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check final neighborhood data\n",
    "neigh.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6) Export processed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to save outputs \n",
    "#trees_clean.to_file('.../trees_madrid_clean.geojson', driver='GeoJSON')\n",
    "#neigh.to_file('.../neighborhoods_madrid_stats.geojson', driver='GeoJSON')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}