{"cells":[{"metadata":{},"cell_type":"markdown","source":"# **Tracking Tree Canopy Equality with Green Coverage Model**"},{"metadata":{},"cell_type":"markdown","source":"![](https://www.nationalgeographic.com/content/dam/science/2020/07/06/crown-shyness/crownshyness_mm9404_200619_00101.ngsversion.1594031243744.adapt.1900.1.jpg)\n*image source*: https://www.nationalgeographic.com/science/2020/07/tree-crown-shyness-forest-canopy/"},{"metadata":{},"cell_type":"markdown","source":"There is a multitude of benefits to tree coverage. Trees regulate extreme temperatures and manage flash flooding. They also absorb carbon dioxide from the atmosphere, helping to offset greenhouse gas emissions.\n\nGiven these benefits, it is clear that trees have an implicit value. In this notebook, we explore tree coverage inequality within cities and how that perpetuates existing imbalances between the rich and the poor. The imbalances we discuss are public health outcomes related to extreme heat and higher vulnerability to flooding. No other category of hazardous weather event in the United States has caused more fatalities over the last few decades than extreme heat.\n\nThe **KPI** we suggest to CDP is the **Gini Coefficient** for a 20 km x 20 km region around a central coordinate. We calculate this by querying satellite images for the area of interest, then extract a subset of the green spectrum to find the density of green coverage. We then analyse how evenly the coverage is spread to calculate the distribution of greenery."},{"metadata":{},"cell_type":"markdown","source":"\n# Executive Summary\n\n<a></a>\n\n* Have created an automated tool which takes a city name as input to output a green coverage map and Lorenz Curve\n\n* Provides actionable insights for any city by showing which areas are underrepresented with tree coverage and whether the tree distribution is equitable\n\n* Creates demand for low skilled work. A group severely impacted by the Covid-19 pandemic\n\n* Is an effective infrastructure project that can be started with relative ease to kick start an economy \n\n* Removes systematic bias towards minorities and lower-income groups\n\n* As global temperatures and extreme weather rise, improves public health outcomes for low socio-economic groups\n\nThe kernel is ordered as follows. We start with an exploratory analysis of the CDP disclosure data. Next, a detailed explanation of our method alongside examples of the application being used with Los Angeles and Sydney to illustrate the robustness of the product. We then elaborate on the advantages and shortcomings of the approach and discuss further iterations of the tool."},{"metadata":{},"cell_type":"markdown","source":"# Exploratory Data Analysis on CDP Disclosure Data\n\nTo understand what KPIs would add value to cities, we explored the disclosure data provided by CDP."},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# Library Imports\n\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom IPython.core.display import HTML\nfrom IPython.display import  Markdown\nimport seaborn as sns\nimport random\nfrom sklearn.cluster import KMeans\nfrom sklearn.preprocessing import StandardScaler\nimport geopandas as gpd\nfrom shapely.geometry import Point\nfrom  sklearn.metrics import davies_bouldin_score\nfrom sklearn.decomposition import PCA\nfrom matplotlib.path import Path\nfrom matplotlib.spines import Spine\nfrom matplotlib.projections.polar import PolarAxes\nfrom matplotlib.projections import register_projection\nfrom matplotlib.patches import Circle, RegularPolygon\nfrom matplotlib.path import Path\nfrom matplotlib.projections.polar import PolarAxes\nfrom matplotlib.projections import register_projection\nfrom matplotlib.spines import Spine\nfrom matplotlib.transforms import Affine2D\nfrom sklearn import tree\nimport statsmodels.api as sm\nfrom scipy import stats\nimport json\nfrom sklearn.feature_extraction.text import CountVectorizer\nfrom sklearn.feature_extraction.text import TfidfVectorizer\nimport plotly.graph_objects as go\nimport json\nimport torch\nimport torchvision.transforms as transforms\nimport os\nimport requests\nimport shutil\nfrom PIL import Image\nfrom tqdm import tqdm\nimport cv2\nimport folium\nfrom folium.plugins import HeatMap\nfrom IPython.display import display\nfrom kaggle_secrets import UserSecretsClient\n\ngoogle_api_key = UserSecretsClient().get_secret(\"google_api_key\")\n\n#centering figures\nHTML(\"\"\"<style> .output_png { display: table-cell; text-align: center; vertical-align: middle;}</style>\"\"\")\n\n#my colors\ncolors= ['#458B00','#629632','#397D02','#567E3A','#A6D785','#687E5A','#8AA37B','#476A34','#7BBF6A','#3D8B37','#426F42','#215E21']\n\n#remove pandas display limit\npd.options.display.max_colwidth = None","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# Data Imports\n\ncities_2018 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2018_Full_Cities_Dataset.csv\",\n                         usecols=[c for c in list(pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2018_Full_Cities_Dataset.csv\", nrows =1)) if c not in ['Questionnaire','Year Reported to CDP','File Name','Last update','Comments']])\ncities_2019 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2019_Full_Cities_Dataset.csv\",\n                         usecols=[c for c in list(pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2019_Full_Cities_Dataset.csv\", nrows =1)) if c not in ['Questionnaire','Year Reported to CDP','File Name','Last update','Comments']])\ncities_2020 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2020_Full_Cities_Dataset.csv\",\n                         usecols=[c for c in list(pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Responses/2020_Full_Cities_Dataset.csv\", nrows =1)) if c not in ['Questionnaire','Year Reported to CDP','File Name','Last update','Comments']])\n\ncompanies_climate_change_2018 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2018_Full_Climate_Change_Dataset.csv\",\n                         usecols=[c for c in list(pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2018_Full_Climate_Change_Dataset.csv\", nrows =1)) if c not in ['Questionnaire','Year Reported to CDP','File Name','Last update','Comments']])\ncompanies_climate_change_2019 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2019_Full_Climate_Change_Dataset.csv\",\n                         usecols=[c for c in list(pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2019_Full_Climate_Change_Dataset.csv\", nrows =1)) if c not in ['Questionnaire','Year Reported to CDP','File Name','Last update','Comments']])\ncompanies_climate_change_2020 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2020_Full_Climate_Change_Dataset.csv\",\n                         usecols=[c for c in list(pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Corporations/Corporations Responses/Climate Change/2020_Full_Climate_Change_Dataset.csv\", nrows =1)) if c not in ['Questionnaire','Year Reported to CDP','File Name','Last update','Comments']])\n\ngeo_cities_2020 = pd.read_csv(\"/kaggle/input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2020_Cities_Disclosing_to_CDP.csv\", usecols=[\"Account Number\",\"City Location\", \"Organization\", \"City\"]) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First, we investigated the most significant climate hazards facing cities today. As shown below, flooding, extreme heat and extreme precipitation account for over 50% of the hazards reported."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Isolating Cities Question\n\nquestion_3_responses = cities_2020[\n    (cities_2020['Question Number'] == '2.1') &\n    (cities_2020['Column Name'] == 'Climate Hazards')\n]\n\nquestion_3_responses = question_3_responses['Response Answer'].str.split(\">\", n = 1, expand = True).loc[:,0]\nquestion_3_response_count = question_3_responses.value_counts()\n\n# Isolating Coorperations Question\n\nquestion_2_responses = companies_climate_change_2020[\n    (companies_climate_change_2020['question_number'] == 'C2.3a') &\n    (companies_climate_change_2020['column_number'] == 3)\n]\n\nrisk_types = [\n    'Acute physical',\n    'Emerging regulation',\n    'Current regulation',\n    'Chronic physical',\n    'Market',\n    'Reputation',\n    'Technology'\n]\n\nquestion_2_responses = question_2_responses[~question_2_responses['response_value'].isin(risk_types)]\nquestion_2_response_count = question_2_responses['response_value'].value_counts().head(10)\n\n# Graphing\n\nfig = plt.figure(figsize=(15,9))\nax1 = fig.add_subplot(121)\n\nquestion_3_response_count.plot.pie(\n    textprops={'color':\"w\"},\n    pctdistance=0.7,\n    autopct='%.2f%%',\n    colors=colors, \n    labels=None,\n    ax=ax1,\n    ylabel=\"\"\n)\n\n\nax1.title.set_text(\"Climate Hazards Identified By Cities\")\nax1.legend(\n    question_3_response_count.index,\n    loc=\"lower center\", \n    bbox_to_anchor=(0.5, -0.4)\n)\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our next step was to understand what actions and goals cities are implementing to increase resilience to these hazards. Vectorising responses to Question 3.3 (2020)\n\n“Please describe the main goals of your city’s adaptation efforts and the metrics / KPIs for each goal.”\n\ngave us the results you see below. Interestingly, the two most common actions were ‘tree planting creation green’ and ‘planting creation green space’. To understand why these actions are so common, we then evaluated which climate issues the discussed actions addressed:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\npd.set_option('mode.chained_assignment', None)\n\nadaptation_actions_ques_3 = cities_2020[\n    (cities_2020['Question Number'] == '3.0') &\n    (cities_2020['Column Number'] == 2)\n]\n\nadaptation_actions_ques_3.loc[:, ['Response Answer']] = adaptation_actions_ques_3.loc[:, ['Response Answer']].fillna('No Response')\n\ndef get_top_n_bigram(corpus, n=None):\n    vec = TfidfVectorizer(ngram_range=(4,4), stop_words='english').fit(corpus)\n    bag_of_words = vec.transform(corpus)\n    sum_words = bag_of_words.sum(axis=0) \n    words_freq = [(word, sum_words[0, idx]) for word, idx in vec.vocabulary_.items()]\n    words_freq =sorted(words_freq, key = lambda x: x[1], reverse=True)\n    return words_freq[:n]\ncommon_words = get_top_n_bigram(adaptation_actions_ques_3['Response Answer'], 20)\n\nadaptation_actions_ques_3 = pd.DataFrame(common_words, columns = ['word' , 'count'])\n  \nadaptation_actions_ques_3  = adaptation_actions_ques_3.set_index('word')[:7]  ## Taking the first 5\n\n\n\nadaptation_goals_ques_3 = cities_2020[\n    (cities_2020['Question Number'] == '3.3') &\n    (cities_2020['Column Number'] == 1) \n]\nadaptation_goals_ques_3.loc[:, ['Response Answer']] = adaptation_goals_ques_3.loc[:, ['Response Answer']].fillna('No Response')\n\ndef get_top_n_bigram(corpus, n=None):\n    vec = TfidfVectorizer(ngram_range=(4,4), stop_words='english').fit(corpus)\n    bag_of_words = vec.transform(corpus)\n    sum_words = bag_of_words.sum(axis=0) \n    words_freq = [(word, sum_words[0, idx]) for word, idx in vec.vocabulary_.items()]\n    words_freq =sorted(words_freq, key = lambda x: x[1], reverse=True)\n    return words_freq[:n]\ncommon_words = get_top_n_bigram(adaptation_goals_ques_3['Response Answer'], 20)\n# for word, freq in common_words:\n#     print(word, freq)\nadaptation_goals_ques_3 = pd.DataFrame(common_words, columns = ['word' , 'count'])\n\nadaptation_goals_ques_3.loc[adaptation_goals_ques_3.word=='reforzar el sistema salud','word']= 'make a stronger health care'\nadaptation_goals_ques_3.loc[adaptation_goals_ques_3.word=='monitoramento risco em tempo','word']= 'monitoring risk in real time'\n\nadaptation_goals_ques_3  = adaptation_goals_ques_3.set_index('word')[:6]  ## Taking the first 5\n\n# Graphing\n\nfig = plt.figure(figsize=(15,9))\nax1 = fig.add_subplot(121)\nax2 = fig.add_subplot(122)\n\nadaptation_actions_ques_3['count'].plot.pie(\n    textprops={'color':\"w\"},\n    pctdistance=0.7,\n    explode=(0.1,0.1,0,0,0,0,0),\n    autopct='%.2f%%',\n    colors=colors, \n    labels=None,\n    ax=ax1,\n    ylabel=\"\"\n)\n\nadaptation_goals_ques_3['count'].plot.pie(\n    textprops={'color':\"w\"},\n    explode=(0.1,0,0,0,0,0),\n    pctdistance=0.7,\n    autopct='%.2f%%',\n    colors=colors, \n    labels=None,\n    ax=ax2,\n    ylabel=\"\"\n)\n\nax1.title.set_text(\"Actions - Present\") ## present \nax2.title.set_text(\"Goals/KPI's - Future\") ## future\nax1.legend(\n    adaptation_actions_ques_3.index,\n    loc=\"lower center\", \n    bbox_to_anchor=(0.5, -0.4)\n)\nax2.legend(\n    adaptation_goals_ques_3.index,\n    loc=\"lower center\", \n    bbox_to_anchor=(0.5, -0.4)\n)\n\nplt.show()\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As the above visualisation shows, improved green coverage increases resilience towards each of the top three risks we discovered earlier, clearly justifying its popularity."},{"metadata":{},"cell_type":"markdown","source":"By selecting these actions, we were then able to observe the metrics used to measure their implementation and success. As shown below, this analysis showed us that cities do not have a common KPI to measure the same objective.\n\nUltimately, our processing and visualisation of disclosure data revealed the demand for a unified green coverage KPI that could be easily used by any city, anywhere in the world. Such a tool would allow cities around the world to understand their urban greenery and compare their tree equality distribution with other cities. As a global organisation who has set themselves apart in disclosure, this KPI is perfect for CDP."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"question_3_responses = cities_2020[cities_2020['Question Number'] == '3.3'].copy()\n\nquestion_3_3_responses = cities_2020[(cities_2020['Question Number'] == '3.3') &\n                                     (cities_2020['Column Number'] == 1)].copy()\n\nquestion_3_3_responses['Response Answer'] = question_3_3_responses['Response Answer'].fillna('No Response')\n \n# Filter question 3.3 to cities that mention green spaces in adaptation goals\n\n# list_of_strings = ['tree', 'plant', 'park', 'planting', 'canopy', 'canopies']\nlist_of_strings = ['tree']\n\n\npattern = '|'.join(list_of_strings)\n\nquestion_3_3_filtered = question_3_3_responses[question_3_3_responses['Response Answer'].str.contains(pattern)]\n\n# Retrieve all columns for the rows that have the list of strings in thier adaptation goals.\n\nall_green_responses = pd.DataFrame()\n\nfor index, row in question_3_3_filtered.iterrows():\n\n    all_green_responses = all_green_responses.append(question_3_responses[(question_3_responses['Account Number'] == row['Account Number']) &\n                                          (question_3_responses['Row Number'] == row['Row Number'])], ignore_index=True)\n    \nall_green_responses_2 = all_green_responses['Response Answer'][all_green_responses['Column Number'] == 2]\n\nall_green_responses_2 = all_green_responses_2.str.split(\">\", n=1, expand=True).loc[:, 0].value_counts()\n\nfig = go.Figure(data=[go.Sankey(\n    node = dict(\n      pad = 15,\n      thickness = 20,\n      line = dict(color = \"#215E21\", width = 0.5),\n      label = [\"Tree\", \"Extreme hot temperature\", \"Flood and sea level rise\", \"Extreme Precipitation\", \"Storm and wind\", \"Water Scarcity\", \"Chemical change\", \"Extreme cold temperature\", \"Biological hazards\", \"Mass movement\"],\n      color = \"#4BB74C\"\n    ),\n    link = dict(\n      source = [0, 0, 0, 0, 0, 0, 0, 0, 0], # indices correspond to labels, eg A1, A2, A1, B1, ...\n      target = [1, 2, 3, 4, 5, 6, 7, 8, 9],\n      value = [20, 10, 8, 6, 5, 3, 2, 2, 1],\n    color = \"#C1FFC1\"\n  ))])\n\nfig.update_layout(title_text=\"Hazard related to Tree\", font_size=15)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n# The intersection between the tree coverage and social inequality\n\n## Health outcomes \n\nMedical and health researchers have shown that fatalities during heatwaves are most commonly due to respiratory and cardiovascular diseases, primarily from heats negative effect on the cardiovascular system. In an attempt to control one's internal temperature, the body's instinct is to circulate large quantities of blood to the skin. However, to perform this protective measure against overheating actually harms the body by inducing extra strain on the heart. This excess strain has the potential to trigger a cardiac event in those with chronic health problems, such as the elderly (Cui et al.). Frumkin showed that the relationship between mortality and temperature creates a J-shaped function, showing a steeper slope at higher temperatures. Records show that more casualties have resulted from heat waves than hurricanes, floods, and tornadoes together. The significance of this statistic is that extreme heat events are so deadly and are becoming more prevalent (Stone et al.).\n\nA study held from 1989 to 2000 have also recorded a rise of 5.7% in mortality during heatwaves. It revealed that Rome's elderly population endures a higher mortality rate during heat waves, at 8% excess for the 65–74 age group and 15% for above 74 (Schifano et al.). Another study found that in French cities during the 2003 heatwave, small towns saw an average excess mortality rate of 40%, while Paris witnessed an increase of 141%. During this period, a 0.5 °C increase above the average minimum nighttime temperature doubled the risk of death in the elderly (Dousset et al.).\n\nSince the air temperature of urban areas with more trees can be around 4C cooler than those without. On a more local level, the air temperature on a treeless residential street can be 10C higher than a nearby shady street (Sinfield et al.).\n\nA city’s poorest areas tend to have less tree canopy than wealthier areas, a pattern that is especially pronounced on the concrete-dense neighbourhoods where temperature regulation is most important (Ready et al.).\n\nThe areas which have faced systematic inequality with tree coverage are now seeing poorer public health outcomes because of it. Through no fault of their own, at-risk people are now more vulnerable to extreme heat because of where they live.\n\nThe best way to counter this issue is to plant trees in the areas which need it most. By planting trees, city councils can effectively regulate temperatures and therefore improve public health outcomes for vulnerable members of their society. Since global temperatures are rising, this is an issue which will only become more prevalent as time goes on. It is critical for councils to act now, which is what has inspired our KPI.\n\n\n## Danger from flooding\n\nThere are several types of flooding that affect communities:\n- River Flooding: When the amount of water entering a river exceeds its holding capacity and overflows its banks\n- Surface water flooding – when heavy rainfall runs on hard or saturated surfaces without getting absorbed into the ground and collects in low areas damaging properties and misplacing communities\n- Drain and sewer flooding – When heavy rain causes drain and sewers to be blocked\n- Coastal flooding – due to climate change, weather and tidal conditions cause an increase in sea levels affecting coastal neighbourhoods.\n\nPlanting trees and extending greenspaces help by intercepting rainfall which significantly slows the rains speed. Tree canopies can capture about 30% of the rainfall which then evaporates back into the atmosphere without reaching the ground. This can even occur in winter where trees intercept and re-evaporate rainfall. The root systems of trees also allow water to penetrate deeper into the soil, decreasing surface run-off while also increasing the water storage capacity of the soil. In urban spaces, an increase in impermeable surfaces such as roads, pavements and driveways has led to increased surface water run-off. Trees reduce surface run-off by 80% compared to asphalt. Incorporating green spaces in urban spaces would drastically reduce run-off leading to a decreased risk of flooding (Woodland).\n\n\n## How can this pull cities out of a recession without perpetuating socical inequities?\n\nTree planting projects can help pull cities out of a recession by mobilising the hardest-hit workers from the Covid-19 pandemic, low skilled workers. By getting this group employed, the economy can pull itself out of the recession as it will put money into the hands of people who need it most. This will stimulate local economies by allowing consumer spending to increase. Being a labour intensive project that requires practically no qualifications, it is a great way to get the masses employed in any city in the world. On top of this, trees often need at least two years of constant care (Ready et al.). Which means these employees can stay hired until the economy has rebounded, and the labour market is thriving again. Now, contrast this with investment in renewable energy. Of course, renewables are fantastic, but the workers involved are typically highly skilled and therefore are most likely still in the job market. If they are not, this is not a big issue as there are plenty of opportunities available to them. By investing in renewables, we are only worsening the income inequality as the low skilled workers are left behind.\n\n\n## How can this be done in a socialy equitable way with the backdrop of a global pandemic?\n\nAs mentioned above, by mobilising low skilled workers for these projects, cities can put money in the hands of people who need it most. This contributes to making society more equitable since it balances income distribution. With industries and hospitality practically reaching a standstill with the pandemic, we can help these newly unemployed workers find meaningful work again.\n\n\n## How can corporations help solve this problem?\n\nCities could outsource the labour requirement to corporations as it is not a councils expertise to directly hire and manage labour oriented workforces.\n\n\n## How does this measure the intersection in the context of resiliency?\n\nThe intersection allows cities to be resilient in a myriad of different ways. By planting more trees, cities are defending themselves against rising temperatures, they are protecting the health of their most vulnerable citizens and are mitigating flash flooding events.\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Tree Canopy Equality KPI\n\nOur KPI allows cities around the world to evaluate their green coverage consistently and reliably. We have created an interface for collecting and assessing satellite images for a given city to understand the equality of their greenery distribution.\n\n## Method:\n\nWe follow the following steps:\n1. Take a city name as input (must be the same as that sent to CDP)\n2. Retrieve coordinates of the city centre from CDP data\n3. Fetch zoomed out satellite images\n4. Determine city boundaries using a custom ML model\n5. Collect thousands of zoomed-in images within the city bounds and asses green coverage\n6. Combine results to form a Lorenz curve and calculate the Gini coefficient\n7. Visualise areas with least tree coverage on a heat map\n\n\nCode:\n\nWe have encapsulated this logic in the following blocks and included them below:\n1. The geographical grid search system\n2. Machine learning-based city detection\n3. Satellite imagery pipeline\n4. Complete pipeline\n"},{"metadata":{},"cell_type":"markdown","source":"## Green Coverage Tool"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# 1. Geographical grid search system \n\ndef getBoundingBoxes(bounds, zoom=14):\n    \n    resolution_map = {\n        14: 0.054885,\n        18: 0.003251\n    }\n    resolution = resolution_map[zoom]\n\n    xSorted = sorted(bounds, key=lambda p: p[0])\n    minX = xSorted[0][0]\n    maxX = xSorted[-1][0]\n\n    ySorted = sorted(bounds, key=lambda p: p[1])\n    minY = ySorted[0][1]\n    maxY = ySorted[-1][1]\n\n\n    gridWidth = int((maxX - minX) / resolution)\n    gridHeight = int((maxY - minY) / resolution)\n    return [\n        [\n            [minX + x * resolution, minY + y * resolution],\n            [minX + (x + 1) * resolution, minY + y * resolution],\n            [minX + (x + 1) * resolution, minY + (y + 1) * resolution],\n            [minX + x * resolution, minY + (y + 1) * resolution],\n            [minX + x * resolution, minY + y * resolution],\n        ]\n        for x in range(gridWidth)\n        for y in range(gridHeight)\n    ]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# 2. Machine learning based city detection\n# Only the functions needed for inference are included within this document\n# All training was completed in a seperate directory to avoid complexity balooning\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nMODEL_PATH = \"/kaggle/input/urban-classification-model/urban_classifier.pth\"\nmodel=torch.load(MODEL_PATH, map_location=torch.device('cpu'))\nmodel.eval()\n\n\nclass BasicClassificationDataset(torch.utils.data.Dataset):\n    def __init__(self, dir=None, image_paths=None):\n        self.dir = dir\n        self.image_paths = image_paths\n        self.transforms = transforms.Compose(\n            [\n                transforms.ToTensor(),\n                transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5)),\n                # transforms.RandomHorizontalFlip(0.5),\n                # transforms.RandomVerticalFlip(0.5),\n            ]\n        )\n        if dir:\n            self.labels = os.listdir(dir)\n            self.images = []\n\n            for i in range(len(self.labels)):\n                label = self.labels[i]\n                images = os.listdir(os.path.join(dir, label))\n                self.images += [\n                    [os.path.join(dir, label, image), i] for image in images\n                ]\n        if image_paths:\n            self.labels = [\"None\"]\n            self.images = list(map(lambda p: [p, 0], self.image_paths))\n\n    def __getitem__(self, idx):\n        image_path, image_class = self.images[idx]\n        img = Image.open(image_path).convert(\"RGB\")\n        img = self.transforms(img)\n        return img, image_class\n\n    def __len__(self):\n        return len(self.images)\n        \ndef urban_classification(image_paths):\n    dataset = BasicClassificationDataset(image_paths=image_paths)\n    loader = torch.utils.data.DataLoader(\n        dataset, batch_size=1, shuffle=False, num_workers=0\n    )\n\n    classifications = []\n\n    with torch.no_grad():\n        for i, data in enumerate(loader, 0):\n            images, labels = data\n            images = images.to(device)\n            outputs = model(images)\n            _, predicted = torch.max(outputs.data, 1)\n            classifications += predicted.tolist()\n\n    return classifications\n\n## Basic color detection for trees\n## If given more time we would train a segmentation model for the task\n\ndef get_green_mask(image_path):\n    img = Image.open(image_path).convert('RGB')\n    arr = np.array(np.asarray(img))\n\n    R = [(36,86),(0,255),(0,255)]\n    red_range = np.logical_and(R[0][0] < arr[:,:,0], arr[:,:,0] < R[0][1])\n    green_range = np.logical_and(R[1][0] < arr[:,:,0], arr[:,:,0] < R[1][1])\n    blue_range = np.logical_and(R[2][0] < arr[:,:,0], arr[:,:,0] < R[2][1])\n    valid_range = np.logical_and(red_range, green_range, blue_range)\n    mask = valid_range\n    # print(valid_range)\n    green_score = mask.sum() / (600 * 600)\n    return (img, mask, green_score)\n\ndef display_mask(image_path):\n    img, mask, green_score = get_green_mask(image_path)\n\n    print(green_score)\n    fig = plt.figure(figsize=(15, 9))\n    ax1 = fig.add_subplot(121)\n    ax2 = fig.add_subplot(122)\n\n    ax1.imshow(mask)\n    ax2.imshow(img)\n    plt.show()\n    ","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# 3. Satelite imagery pipeline\n\ndef satelite_download(tile, output_dir, zoom):\n    size = [600, 600]\n    zoom = zoom\n\n    centerX = (tile[0][0] + tile[2][0]) / 2\n    centerY = (tile[0][1] + tile[2][1]) / 2\n\n    image_url = f\"http://maps.googleapis.com/maps/api/staticmap?center={centerY},{centerX}&zoom={zoom}&size={size[0]}x{size[1]}&maptype=satellite&key={google_api_key}\"\n    filename = f\"{centerX},{centerY}.png\"\n    filename = os.path.join(output_dir, filename)\n    \n    #     print(image_url)\n\n    if os.path.exists(filename):\n        return filename\n    else:\n        # Open the url image, set stream to True, this will return the stream content.\n        r = requests.get(image_url, stream=True)\n\n        # Check if the image was retrieved successfully\n        if r.status_code == 200:\n            # Set decode_content value to True, otherwise the downloaded image file's size will be zero.\n            r.raw.decode_content = True\n\n            # Open a local file with wb ( write binary ) permission.\n            with open(filename, \"wb\") as f:\n                shutil.copyfileobj(r.raw, f)\n\n#             print(\"Image sucessfully Downloaded: \", filename)\n            return filename\n        else:\n#             print(\"Image Couldn't be retreived\")\n            return None","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# 4. Complete pipeline\n\ndef get_green_scores(city_name, data_dir, output_path):\n\n    if (google_api_key=='') or (google_api_key==None):\n        print('Google Maps API Key Needed')\n        return {}\n    \n    if (city_name not in list(geo_cities_2020[\"City\"])):\n        print('City name not found in CDP dataset')\n        print('May have different name like City of {}'.format(city_name))\n        return {}\n    \n    print(list(geo_cities_2020[geo_cities_2020[\"City\"] == city_name][\"City Location\"]))\n#     if (list(geo_cities_2020[geo_cities_2020[\"City\"] == city_name][\"City Location\"])):\n#         print('City name not found in CDP dataset')\n#         print('May have different name like City of {}'.format(city_name))\n#         return {}\n    \n    longitude, latitude  = geo_cities_2020[geo_cities_2020[\"City\"] == city_name]\\\n        [\"City Location\"].iloc[0].replace('POINT (', '').replace(')', '').split(\" \")\n    longitude, latitude = float(longitude), float(latitude)\n\n    width = 0.7\n    city_surrounding_area = [\n        [longitude - width/2,latitude - width/2],\n        [longitude + width/2,latitude + width/2]\n    ]\n    \n    #     print(city_surrounding_area)\n    #     print(\"city_surrounding_area\")\n    large_grids = getBoundingBoxes(city_surrounding_area, zoom=14)\n    large_image_paths = [satelite_download(grid, data_dir, 14) for grid in tqdm(large_grids, desc=\"Downloading large images\")]\n    \n    #     print(\"large_image_paths\")\n    #     print(large_image_paths)\n    urban_grid_mask = urban_classification(large_image_paths)\n    urban_grids = [large_grids[i] for i in range(len(urban_grid_mask)) if (urban_grid_mask[i]==0)]\n\n    small_grids = []\n    for grid in urban_grids:\n        small_grids += getBoundingBoxes(grid, zoom=18)\n    small_image_paths = [satelite_download(grid, data_dir, 18) for grid in tqdm(small_grids, desc=\"Downloading small images\")]\n    \n    mapped_data = {}\n        \n    count = 0\n    for path in tqdm(small_image_paths, desc=\"Segmenting images\"):\n        try:\n            file_name = path.split(\"/\")[-1]\n        except:\n            print(\"Something went wrong\")\n\n        if (count % 10 == 0):\n            with open(output_path, \"w\") as outfile:\n                json.dump(mapped_data, outfile)\n                #                 print(\"saved\")\n\n\n        if file_name not in mapped_data:\n            centerX = float(file_name.split(\",\")[1].replace(\".png\", \"\"))\n            centerY = float(file_name.split(\",\")[0])\n            _, mask, green_score = get_green_mask(path)\n            mapped_data[file_name] = green_score\n            \n        count += 1\n\n    with open(output_path, \"w\") as outfile:\n        json.dump(mapped_data, outfile)\n\n    return mapped_data","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"def make_map(green_scores, city_name):\n\n    sorted_green_scores = sorted(list(green_scores.values()))\n    Q1 = sorted_green_scores[int(len(sorted_green_scores)/6)]\n\n    lowest_green_scores = []\n    for file_name in green_scores:\n        centerX = float(file_name.split(\",\")[1].replace(\".png\", \"\"))\n        centerY = float(file_name.split(\",\")[0])\n        if green_scores[file_name]<Q1:\n            lowest_green_scores.append([centerX, centerY, green_scores[file_name]])\n\n    longitude, latitude  = geo_cities_2020[geo_cities_2020[\"City\"] == city_name]\\\n        [\"City Location\"].iloc[0].replace('POINT (', '').replace(')', '').split(\" \")\n    longitude, latitude = float(longitude), float(latitude)\n\n    base_map = folium.Map(location=[latitude, longitude], control_scale=True, zoom_start=12, max_zoom=13, min_zoom=10)\n    HeatMap(data=lowest_green_scores, radius=10, max_zoom=20, gradient={0.0: 'red', 0.05: 'red'}).add_to(base_map)\n    return base_map","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def graph_lorenz(green_scores, city_name):\n    fig = plt.figure(figsize=(15,9))\n    ax1 = fig.add_subplot(111)\n    green_scores = list(green_scores.values())\n    y = np.array(green_scores)\n    lorenz_curve(y, ax1)\n    ax1.title.set_text(f\"{city_name} - Gini Coefficient = {str(gini(y))[:5]}\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Example: Los Angeles"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"city_name = 'Los Angeles'\n\nsaved_la_data_path = \"/kaggle/input/green-score-tree-canopy/green_scores_la.json\"\nwith open(saved_la_data_path) as json_file:\n    la_green_scores = json.load(json_file)\n\ndisplay(make_map(la_green_scores, city_name))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Example: Sydney"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"city_name = 'Sydney'\n\nwith open(\"/kaggle/input/green-score-tree-canopy/green_scores_sydney.json\") as json_file:\n    syd_mapped_data = json.load(json_file)\n    \nlowest_green_scores = []\nfor file_name in syd_mapped_data:\n    centerX = float(file_name.split(\",\")[1].replace(\".png\", \"\"))\n    centerY = float(file_name.split(\",\")[0])\n    if syd_mapped_data[file_name]<0.1:\n        lowest_green_scores.append([centerX, centerY, syd_mapped_data[file_name]])\n\nbase_map = folium.Map(location=[-33.9358211, 151.03409150000005], control_scale=True, zoom_start=12, max_zoom=13, min_zoom=10)\nHeatMap(data=lowest_green_scores, radius=10, max_zoom=20, gradient={0.0: 'red', 0.05: 'red'}).add_to(base_map)\nbase_map","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Lorenz Curve and Gini Coefficient\n\nTo allow cities to measure their overall Tree Canopy Equality, we calculate a Lorenz curve in realtime so councils can get a graphical representation of what their unique distribution is. This curve allows cities to see the proportion of overall tree canopy assumed by the bottom x% of the land area. *(An individual land area is regarded as a 50 m x 50 m block of land)*\n\nWe compute the Gini Coefficient by calculating the area between the line of equality and the Lorenz curve. This final value is the **KPI** which cities want to optimise as it means they have equitable tree coverage throughout their city.\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"with open(\"/kaggle/input/green-score-tree-canopy/green_scores_sydney.json\") as json_file:\n    syd_mapped_data = json.load(json_file)\nwith open(\"/kaggle/input/green-score-tree-canopy/green_scores_la.json\") as json_file:\n    la_mapped_data = json.load(json_file)\n\nsyd_ma = list(syd_mapped_data.values())\nla_ma = list(la_mapped_data.values())\n\ndef gini(arr):\n    ## first sort\n    sorted_arr = arr.copy()\n    sorted_arr.sort()\n    n = arr.size\n    coef_ = 2. / n\n    const_ = (n + 1.) / n\n    weighted_sum = sum([(i+1)*yi for i, yi in enumerate(sorted_arr)])\n    return coef_*weighted_sum/(sorted_arr.sum()) - const_\n\n\ndef lorenz_curve(arr, ax):\n    sorted_arr = arr.copy()\n    sorted_arr.sort()\n    X_lorenz = sorted_arr.cumsum() / sorted_arr.sum()\n    X_lorenz = np.insert(X_lorenz, 0, 0)\n    X_lorenz[0], X_lorenz[-1]\n    ## scatter plot of Lorenz curve\n    ax.scatter(np.arange(X_lorenz.size) / (X_lorenz.size - 1), X_lorenz,\n               marker='x', color='darkgreen', s=2)\n    ## line plot of equality\n    ax.plot([0, 1], [0, 1], color='k')\n    \nfig = plt.figure(figsize=(15,9))\nax1 = fig.add_subplot(121)\nax2 = fig.add_subplot(122)\n\n\nx = np.array(syd_ma)\n\nlorenz_curve(x, ax1)\n\ny = np.array(la_ma)\n\nlorenz_curve(y, ax2)\n\nax1.title.set_text(f\"Sydney - Gini Coefficient = {str(gini(x))[:5]}\") ## Sydney \nax2.title.set_text(f\"Los Angeles - Gini Coefficient = {str(gini(y))[:5]}\") ## LA","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Try it yourself\n\nEnter a city name to see a map of where trees are needed most and a Lorenz Curve to see the total equality. To use yourself, uncomment the code below and enter your googe api key in to your kaggle secret keys."},{"metadata":{"trusted":true},"cell_type":"code","source":"# city_name = 'London'\n# output_path = '../../green_scores2.json'\n\n# green_scores = get_green_scores(city_name, './', output_path)\n# if green_scores != {}:\n#     display(make_map(green_scores, city_name))\n#     graph_lorenz(green_scores, city_name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Advantages\n\nThe best benchmark for our Green Coverage Tool and Gini Coefficient is Treepedia by MIT Senseable City Lab. They have tried to capture the green canopy by utilising Google Street view. They scan each panorama and use a predictive model to identify tree canopy size from eye level. The outputs they provide is a map of the city where each street view location as a dot with a varying green intensity depending on how much tree was in the panorama (this shows a similar visualisation to our heat map). They also provide a Green View Index which shows the percentage of canopy coverage for the entire city (similar to our Gini coefficient).\n\nOur solution recognises considerably more greenery, including parklands and backyards. As discussed, these features increase resilience to extreme weather events and therefore, need to be taken into consideration. Our product also has the advantage of directly integrating with CDPs data and customers."},{"metadata":{},"cell_type":"markdown","source":"## Disadvantages\n\nAlthough our work represents a great start, there are still many improvements to be made. These include:\n\n1. Training an ML segmentation model for identifying greenery rather than a colour based approach\n2. Training another classification model to identify area type at a smaller scale\n3. Create an online interface that allows organisations reporting to CDP to easily view and interact with the data\n\nThese steps are possible and would significantly improve the value of our product, however, time constraints prevented us from pursuing them."},{"metadata":{},"cell_type":"markdown","source":"# Final Thoughts\n\nOverall this was a great exercise in predicting tree canopy inequality in cities. Although we see much more room for growth and improvement with our models and visualisations, we have created a great tool that is ready to be used by cities to see where to plant trees and how equal their current tree distribution is.\n\nWe believe the generic nature of the tool strongly complements CDP's mission for global disclosure as it can be applied to any city with relative ease."},{"metadata":{},"cell_type":"markdown","source":"\n## Our team\n\nWe had the pleasure of working with an experienced and diverse team. Scattered across the world, we connected over zoom, slack and git:\n\n<table>\n  <tr>\n    <td><img src=\"https://media-exp1.licdn.com/dms/image/C5603AQHbEQCGU89V6A/profile-displayphoto-shrink_400_400/0/1596957301931?e=1611792000&v=beta&t=fZn_O9amBy7bHGYLu9sd1XEIwjURey3TW-AWQ5Mtn1w\" width=\"150\"></td>\n      <td style='text-align: left;'><strong>Adrian Sarstedt</strong><br/><i>Data Engineer at the Florey Institute</i><br/><i>Australia</i></td>\n      <td><img src=\"https://media-exp1.licdn.com/dms/image/C5603AQHVTBlErFKspw/profile-displayphoto-shrink_400_400/0/1551742466431?e=1611792000&v=beta&t=wBO3P2BMWHC39dkXI3Kkjq5gKV4mZLoFeCZs20YaomQ\" width=\"150\"></td>\n    <td style='text-align: left;'><strong>Hamish Gunasekara</strong><br/><i>Data/Risk Analyst at Afterpay</i><br/><i>Australia</i></td>\n  </tr>\n    <tr>\n    <td><img src=\"https://media-exp1.licdn.com/dms/image/C4D03AQHeh45CPg4YsA/profile-displayphoto-shrink_400_400/0?e=1611792000&v=beta&t=IDlea0_mutyIpY6FKT4Btbtc3DtmAB8TFSoAkdfejWE\" width=\"150\"></td>\n    <td style='text-align: left;'><strong>Adham Al Hossary</strong><br/><i>Data Scientist at C-Capture</i><br/><i>England</i></td>\n      <td><img src=\"https://media-exp1.licdn.com/dms/image/C4D03AQFUX0Sfw87-JQ/profile-displayphoto-shrink_400_400/0?e=1611792000&v=beta&t=klWv7UqcQhpXRzP2tjgsLOWIdMYlFMd_ZHrCYyuutes\" width=\"150\"></td>\n    <td style='text-align: left;'><strong>Alexandra Golab</strong><br/><i>Business Analyst at CitiBank</i><br/><i>Spain</i></td>\n  </tr>\n</table>"},{"metadata":{},"cell_type":"markdown","source":"# References\n\nThe study on heat and redlining:\nThe Effects of Historical Housing Policies on Resident Exposure to Intra-Urban Heat: A Study of 108 US Urban Areas - Jeremy S. Hoffman, Vivek Shandas and Nicholas Pendleton\nhttps://www.mdpi.com/2225-1154/8/1/12...\n\nInteractive maps of neighborhood heat and redlining:\nhttps://www.arcgis.com/apps/dashboard...\n\nRobert K. Nelson, LaDale Winling, Richard Marciano, Nathan Connolly, et al., “Mapping Inequality,” American Panorama, ed. Robert K. Nelson and Edward L. Ayers, accessed August 4, 2020.\nhttps://dsl.richmond.edu/panorama/red...\n\nAll about Urban Heat Islands from Climate Central [PDF]:\nhttp://assets.climatecentral.org/pdfs...\n\nCan trees and woods reduce flooding:\nhttps://www.woodlandtrust.org.uk/trees-woods-and-wildlife/british-trees/flooding/\n\nSatellite monitoring of summer heat waves in the Paris metropolitan area:\nhttps://rmets.onlinelibrary.wiley.com/doi/10.1002/joc.2222\n\nUrban tree canopy governance and redlined neighborhoods: an analysis of five cities\nhttps://dspace.mit.edu/handle/1721.1/127588\n\nResidential housing segregation and urban tree canopy in 37 US Cities\nhttps://osf.io/preprints/socarxiv/97zcs/\n"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}