<jupyter_start><jupyter_text>[![Open In SageMaker Studio Lab](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/SatelliteVu/SatelliteVu-AWS-Disaster-Response-Hackathon/blob/main/fire_prediction_visualisation/animation_fire_warning_using_open_street_map.ipynb)<jupyter_code># you might need spatialindex for running osmnx
'''
on ubuntu
sudo apt update
sudo apt install libspatialindex-dev

'''
import osmnx

import os
import sys
import time
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import pyproj
import s3fs
import geojson
from shapely.geometry import Polygon,Point
from shapely.ops import transform
from IPython.display import clear_output
from ipywidgets import HTML
from ipyleaflet import Map, basemaps, GeoJSON, GeoData,ImageOverlay,Popup
from ipyleaflet import SplitMapControl,LayersControl,projections,LayerGroup,AwesomeIcon,Marker,LegendControl,FullScreenControl

PIXEL_RES = 500
EPSG_MAP = projections.EPSG3857
colour_gradient_fire_list = [" #A10100","#DA1F05","#F33C04","#FE650D","#FFC11F","#FFF75D"]

def get_geojson_bbox_crs_ID(geojson_path):
    with open(geojson_path) as f:
        gj = geojson.load(f)
    features = gj['features'][0]
    poly = Polygon(features["geometry"]["coordinates"][0])
    crs = gj['crs']['properties']['name']
    crs = crs.split("::")[1]
    return poly,crs

def polygon_utm_to_wgs(poly,crs,reverse=False):
    crs_wgs84 = pyproj.CRS('EPSG:4326')
    crs_utm   = pyproj.CRS(f'EPSG:{crs}')
    if reverse:
        project = pyproj.Transformer.from_crs(crs_wgs84,crs_utm,always_xy=True).transform
    else:
        project = pyproj.Transformer.from_crs(crs_utm,crs_wgs84,always_xy=True).transform
    poly_wgs84 = transform(project, poly)
    return poly_wgs84

def get_osm(poly_wgs84):
    # query overpass API for buildings
    osm_buildings = osmnx.geometries.geometries_from_polygon(poly_wgs84,
                                                          tags={'building': True})
    # clean loaded data to just keep polygons
    building_polygons = osm_buildings[osm_buildings['geometry'].geom_type == 'Polygon']
    return building_polygons

def get_fire_data(fire_npy_path,bbox_raw,crs):
    # read fire data npy data
    # return geopandas dataframe for each pixels
    fire_arr = np.load(fire_npy_path)
    npy_x_list,npy_y_list = np.where(fire_arr==True)
    minx, miny, maxx, maxy = bbox_raw.bounds
    fire_utm_y_list = minx + npy_y_list*PIXEL_RES
    fire_utm_x_list = maxy - npy_x_list*PIXEL_RES
    n_fire = len(fire_utm_y_list)
    fire_lon_list = np.zeros(n_fire)
    fire_lat_list = np.zeros(n_fire)
    
    #place holder for frp
    fire_frp_list = fire_arr[npy_x_list,npy_y_list]
    fire_frp_list = fire_frp_list+np.random.random(n_fire)*200
    
    for ind,(utm_x,utm_y) in enumerate(zip(fire_utm_y_list,fire_utm_x_list)):
        fire_point_wgs84 = polygon_utm_to_wgs(Point(utm_x,utm_y),crs)
        fire_lon_list[ind] = fire_point_wgs84.x
        fire_lat_list[ind] = fire_point_wgs84.y
    df = pd.DataFrame({
        'npy_x':npy_x_list,
        'npy_y':npy_y_list,
        'utm_x':fire_utm_y_list,
        'utm_y':fire_utm_x_list,
        'frp':fire_frp_list,
        'Longitude':fire_lon_list,
        'Latitude':fire_lat_list       
    })
    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude))
    return gdf

def calculate_closest_distance(building_utm_x,building_utm_y,fire_utm_x_list,fire_utm_y_list):
    distance_list = ((fire_utm_x_list-building_utm_x)**2+(fire_utm_y_list-building_utm_y)**2)**0.5
    return np.min(distance_list)

def run_warning_system():
    import random
    random_colour_list = ['#4D8C57','#78A161','#A3B56B','#CDCA74','#F8DE7E']


    for index,row in gdf_osm.iterrows(): # Looping over all points
        building_type = row['building']
        if building_type == 'yes':
            building_type = 'building'
        if pd.notna(row['name']):
            building_type = row['name']
        central_point = row['geometry'].centroid
        central_point_utm = polygon_utm_to_wgs(central_point,crs,True)
        building_utm_x = central_point_utm.x
        building_utm_y = central_point_utm.y
        distance = calculate_closest_distance(building_utm_x,building_utm_y,fire_utm_x_list,fire_utm_y_list)
        clear_output(wait=True)
        thereshold = 10000
        if distance> thereshold:
            continue
        text = f"WARNING, {building_type} at lat/long {round(central_point.y,3)}/{round(central_point.x,3)} is {round(distance/1000,1)}km away from fire!"
        
        print (text)
        
        
        message = HTML()
        message.value = text
        popup = Popup(
            location=(central_point.y+0.03,central_point.x),
            child=message,
            close_button=False,
            auto_close=False,
            close_on_escape_key=False
            )
        m.add_layer(popup)

        # use awesomeicon instead
        icon1 = AwesomeIcon(
            name='building-o',
            marker_color='blue',
            icon_color='white',
            spin=False
        )
        geo_data_osm = Marker(icon=icon1, location=(central_point.y,central_point.x))
        m.add_layer(geo_data_osm)
        
        time.sleep(1)
        m.remove_layer(popup)

# get data
        
fire_id = 2285
geojson_path = os.path.join(str(fire_id),"bbox.geojson")
poly_raw,crs    = get_geojson_bbox_crs_ID(geojson_path)
buffer_zone = 5000 # get anything within 5km
poly = Polygon(poly_raw.buffer(buffer_zone))
poly_wgs84  = polygon_utm_to_wgs(poly,crs)
central_wgs84 = poly_wgs84.centroid
lon = central_wgs84.x
lat = central_wgs84.y
gdf_osm = get_osm(poly_wgs84)

fire_npy_path = os.path.join(str(fire_id),"todays_fires.npy")
gdf_fire = get_fire_data(fire_npy_path,poly_raw,crs)
fire_utm_x_list = gdf_fire['utm_x']
fire_utm_y_list = gdf_fire['utm_y']

print ('Done')

m = Map(center=(lat,lon), 
    zoom = 11, 
    basemap= basemaps.Esri.WorldTopoMap,
    crs = EPSG_MAP)
geo_data_fire = GeoData(geo_dataframe = gdf_fire,
    style={'color': 'black', 'radius':8, 'fillColor': '#A10100', 'opacity':0.5, 'weight':1.9, 'dashArray':'2', 'fillOpacity':0.6},
    hover_style={'fillColor': 'red' , 'fillOpacity': 0.2},
    point_style={'radius': 5, 'color': 'red', 'fillOpacity': 0.8, 'fillColor': 'red', 'weight': 3},
#         name = 'Release'
                  )
m.add_layer(geo_data_fire)

m.add_control(FullScreenControl())
m
# time.sleep(20)
run_warning_system()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11329974041838449, 'fraction_numerical': 0.021224614444953428, 'mean_word_length': 3.770575382374363, 'pattern_counts': {'":': 1, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '462948', 'license': 'Apache-2.0', 'n_tokens_mistral': 2591, 'n_tokens_neox': 2461, 'n_words': 375, 'path': 'fire_prediction_visualisation/animation_fire_warning_using_open_street_map.ipynb', 'repo_name': 'SatelliteVu/SatelliteVu-AWS-Disaster-Response-Hackathon', 'size': '6549'}