02_analyze.py

import base64
import geopandas as gpd
import h3
import json
import plotly.express as px
from pyproj import Geod
import shapely


# Load Things
# ---

my_logo = base64.b64encode(open('logo-dark.png', 'rb').read())
my_logo_b64 = 'data:image/png;base64,{}'.format(my_logo.decode())

gdf = gpd.read_file('census_data_converted.geojson')
zone_df = gpd.read_file('boundaries/zones.geojson')


TOTAL_TREES = gdf.shape[0]
print(f'TOTAL TREES: {TOTAL_TREES}')

ZONE_NAME_TO_ID = {
    'Bommanahalli': 'bommanahalli',
    'East': 'east',
    'West': 'west',
    'South': 'south',
    'Raja Rajeswari Nagar': 'rr_nagar',
    'Mahadevpura': 'mahadevapura',
    'Yelahanka': 'yelahanka',
    'Dasarahalli': 'dasarahalli'
}

ZONE_SHAPES = {}
ZONE_AREAS = {}

for idx, row in zone_df.iterrows():
    geod = Geod(ellps="WGS84")
    zone_id = ZONE_NAME_TO_ID[row['Zone']]
    ZONE_SHAPES[zone_id] = row['geometry']
    ZONE_AREAS[zone_id] = abs(geod.geometry_area_perimeter(row['geometry'])[0]/1000000)


# Common Functions
# ---

def reverse_coord_order(line_string):
    return [p[::-1] for p in line_string]

def gdf_to_cells(gdf):
    polygons = []
    for geom in list(gdf['geometry']):
        obj = json.loads(shapely.to_geojson(geom))
        if obj['type']=='Polygon':
            polygons.append(obj['coordinates'])
        if obj['type']=='GeometryCollection':
            polygons += [g['coordinates'] for g in obj['geometries']]
    
    cells = set()
    for polygon in polygons:
        shell = reverse_coord_order(polygon[0])
        holes = []
        if len(polygon) > 1:
            holes = [reverse_coord_order(ls) for ls in polygon[1:]]
        h3shape = h3.LatLngPoly(shell, *holes)
        cells.update(h3.polygon_to_cells(h3shape, 8))
    
    return list(cells)
            

def idx_to_polygon(idx):
    coords = h3.cell_to_boundary(idx)
    coords = [coord[::-1] for coord in coords]
    return shapely.Polygon(coords)


def draw_choropleth(df, color_column, title, filename):
    geometries = df.geometry.apply(fix_geometry_type)
    fig = px.choropleth(
        df,
        geojson=geometries,
        locations=df.index,
        color=color_column,
        color_continuous_scale=px.colors.sequential.BuGn,
    )

    fig.update_layout(
        font_family = 'Jost',
        title_text = title,
        title_xref='paper', title_yref='paper',
        title_x=0.05, title_y=0.95,
        title_yanchor='middle',
        title_font_family = 'Jost',
        title_font_weight = 500,
    )

    fig.layout.images = [dict(
        source=my_logo_b64,
        xref="paper", yref="paper",
        x=0.0, y=0.95,
        sizex=0.1, sizey=0.1,
        xanchor="center", yanchor="middle"
    )]

    fig.add_annotation(
        xref='paper', yref='paper',
        yanchor='middle',
        x=0.0, y=0.05,
        text='© Karthik Devan (kdqed.com)',
        showarrow=False,
        font=dict(
            family='Jost',
            size=16,
        )
    )
    fig.update_geos(fitbounds="locations", visible=False)
    fig.write_image(filename)

def fix_geometry_type(geom):
    if geom.geom_type == 'GeometryCollection':
        return shapely.MultiPolygon(geom.geoms)
    else:
        return geom
    
# Coverage Map
# ---

cells = gdf_to_cells(zone_df)
h3df = gpd.GeoDataFrame({
    'H3Index8': cells,
    'geometry': list(map(idx_to_polygon, cells)),
    'Tree Count': [0] * len(cells)
}, geometry='geometry')

cell_counts = gdf[['H3Index8']].groupby(by='H3Index8').value_counts().reset_index(name='Count')

def update_count(r):
    h3df.loc[ h3df['H3Index8'] == r['H3Index8'], 'Tree Count' ] = r['Count']        
cell_counts.apply(update_count, axis=1)

draw_choropleth(h3df, 'Tree Count', 
    'Density Map of Trees in Published Data (BBMP Tree Census)',
    'results/density_map.png'
)


# Zone & Ward-wise Count
# ---

zone_grouping = gdf[['Zone']].groupby(by='Zone').value_counts()
zonewise_counts = zone_grouping.reset_index(name='Total Trees')
zonewise_counts.loc[7] = ['dasarahalli', 0]
zonewise_counts['Trees Per Sq. km'] = zonewise_counts.apply(
    lambda row: round(row['Total Trees'] / ZONE_AREAS[row['Zone']]),
    axis = 1,
)
zonewise_counts.sort_values(by='Trees Per Sq. km', ascending=False, inplace=True)
zonewise_counts.to_markdown('results/zonewise_treecount.md', index=False)

zonewise_gdf = gpd.GeoDataFrame(
    zonewise_counts,
    geometry = zonewise_counts['Zone'].apply(lambda x: ZONE_SHAPES[x]),
)

draw_choropleth(zonewise_gdf, 'Trees Per Sq. km', 
    'Zone-wise Tree Density (BBMP Tree Census)',
    'results/zonewise_density_map.png'
)

ward_grouping = gdf[['WardNumber']].groupby(by='WardNumber').value_counts()
wardwise_counts = ward_grouping.sort_values(ascending=False).reset_index(name='Total Trees')
wardwise_counts.to_markdown('results/wardwise_treecount.md', index=False)


# Tree Frequency Count
# ---

tree_grouping = gdf[['TreeName']].groupby(by='TreeName').value_counts()
tree_counts = tree_grouping.sort_values(ascending=False).reset_index(name='Count')
tree_counts['Percentage'] = tree_counts['Count'].apply(
    lambda c: str(round(100*c/TOTAL_TREES)) + '%'
)
tree_counts.to_markdown('results/tree_name_count.md', index=False)


# Identified and Unknown Trees
# ---
igdf = gdf[gdf['TreeName'] != 'Others']
ugdf = gdf[gdf['TreeName'] == 'Others']

# Trees By Zone
# ---

grouped_counts = igdf[['Zone', 'TreeName']].groupby(by=['Zone', 'TreeName']).value_counts()
grouped_counts = grouped_counts.reset_index(name='Count').sort_values(
    by=['Zone', 'Count'],
    ascending=[True, False],
).groupby(by='Zone')

top3 = grouped_counts.head(3)
top3.to_markdown('results/zonewise_top_trees.md', index=False)

grouped_unknown = ugdf[['Zone']].groupby(by='Zone').value_counts()
unknown_counts = grouped_unknown.reset_index(name='Unknown')
unknown_counts = unknown_counts.merge(zonewise_counts, how='inner', on='Zone')
unknown_counts['Percent Unknown'] = round(100*unknown_counts['Unknown']/unknown_counts['Total Trees'], 2)

unknown_counts.sort_values(by='Percent Unknown', ascending=False, inplace=True)
unknown_counts['Percent Unknown'] = unknown_counts['Percent Unknown'].apply(
    lambda x: str(round(x)) + '%'
)
unknown_counts[['Zone', 'Unknown', 'Percent Unknown']].to_markdown('results/zonewise_unknowncount.md', index=False)

# Trees By Ward
# ---

grouped_counts = igdf[['WardNumber', 'TreeName']].groupby(by=['WardNumber', 'TreeName']).value_counts()
grouped_counts = grouped_counts.reset_index(name='Count').sort_values(
    by=['WardNumber', 'Count'],
    ascending=[True, False],
).groupby(by='WardNumber')

top1 = grouped_counts.head(1)
top1.to_markdown('results/wardwise_top_trees.md', index=False)

grouped_unknown = ugdf[['WardNumber']].groupby(by='WardNumber').value_counts()
unknown_counts = grouped_unknown.reset_index(name='Unknown')
unknown_counts = unknown_counts.merge(wardwise_counts, how='inner', on='WardNumber')
unknown_counts['Percent Unknown'] = round(100*unknown_counts['Unknown']/unknown_counts['Total Trees'], 2)

unknown_counts.sort_values(by='Percent Unknown', ascending=False, inplace=True)
unknown_counts.to_markdown('results/wardwise_unknowncount.md', index=False)