<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/fire_overlay_burned_area_sentinel2.ipynb)Fire prediction data is overlaid on Sentinel 2 imagery<jupyter_code>import os
import sys
import time
import glob
import geojson
import requests
import datetime
import numpy as np
import pandas as pd
import geopandas as gpd
import pyproj
import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon
from shapely.ops import transform
# Polygon function conflict with shapely Polygon, just import ipyleaflet here to avoid it
import ipyleaflet 
# from ipyleaflet import Map, LegendControl,basemaps, GeoJSON, GeoData,ImageOverlay,Heatmap,LocalTileLayer
# from ipyleaflet import ZoomControl,SplitMapControl,LayersControl,projections,FullScreenControl,Polygon
from ipywidgets import Layout

IMAGE_SIZE = 64
GRID_RES   = 500 #m, only works for UTM
EPSG_MAP = ipyleaflet.projections.EPSG3857
PIXEL_RES = 500

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_fire_data(fire_npy_path,bbox_raw,crs):
    # read fire data npy data
    # return geopandas dataframe for each pixels
    if 'probabilities' in fire_npy_path:
        fire_thres = 0.5
        fire_arr = np.load(fire_npy_path)[:,:,0]
    else:
        fire_thres = 0
        fire_arr = np.load(fire_npy_path)
    npy_x_list,npy_y_list = np.where(fire_arr>fire_thres)
    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]
    polygon_list  = []
    
    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
        polygon = []
        four_corners = [[utm_x,utm_y],
                       [utm_x+PIXEL_RES,utm_y],
                       [utm_x+PIXEL_RES,utm_y+PIXEL_RES],
                       [utm_x,utm_y+PIXEL_RES],
                       [utm_x,utm_y]]
        for ix,iy in four_corners:
            point_wgs84 = polygon_utm_to_wgs(Point(ix,iy),crs)
            polygon.append((point_wgs84.y,point_wgs84.x))
        polygon_list.append(polygon)
    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,
        'polygon':polygon_list
    })
    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude))
    return gdf

def get_bbox_from_central_latlon(lat,lon,box_boundary = 0.01):
    left   = lon-box_boundary
    right  = lon+box_boundary
    bottom = lat-box_boundary
    top    = lat+box_boundary
    bbox = f"{left},{bottom},{right},{top}"
    return bbox

def roi_bbox_to_geojson(bbox,buffer = 0):
    left, bottom, right, top = [float(i) for i in bbox.split(',')]
    left   = left-buffer
    right  = right+buffer
    bottom = bottom-buffer
    top    = top+buffer
    roi_poly = [[left,top],
                [right,top],
                [right,bottom],
                [left,bottom],
                [left,top]]
    roi_poly = {'type': 'Polygon', 'coordinates': [roi_poly]}
    return roi_poly
case_list =[37557,37571] # 20211006-20211007
# load predicted data
for ind,test_id in enumerate(case_list):
    geojson_path = os.path.join(str(test_id),"bbox.geojson")
    fire_npy_path = os.path.join(str(test_id),"probabilities.npy")
    poly_raw,crs    = get_geojson_bbox_crs_ID(geojson_path)
    poly_wgs84  = polygon_utm_to_wgs(poly_raw,crs)
    central_wgs84 = poly_wgs84.centroid
    lon = central_wgs84.x
    lat = central_wgs84.y
    gdf_fire = get_fire_data(fire_npy_path,poly_raw,crs)
    if ind == 0:
        gdf_fire_predict = gdf_fire
    else:
        gdf_fire_predict = gdf_fire_predict.append(gdf_fire, ignore_index=True)
print (len(gdf_fire_predict))
# get s2 before and after burning image

bbox = get_bbox_from_central_latlon(lat,lon,0.02)
roi_poly = roi_bbox_to_geojson(bbox)


date_str_preburn  = f"2021-08-01T00:00:00Z/2021-08-20T00:00:00Z"
date_str_postburn = f"2021-10-06T00:00:00Z/2021-10-09T00:00:00Z"

stac_items_preburn = requests.post(f'https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/',
                        json={'intersects': roi_poly, 
#                               'query': {'eo:cloud_cover': {'lt': 5, 'gt':0}},
                             'datetime': date_str_preburn}).json()
preburn_url = stac_items_preburn['features'][0]['assets']['visual']['href']


stac_items_postburn = requests.post(f'https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/',
                        json={'intersects': roi_poly, 
#                               'query': {'eo:cloud_cover': {'lt': 5, 'gt':0}},
                             'datetime': date_str_postburn}).json()
postburn_url = stac_items_postburn['features'][0]['assets']['visual']['href']
print (preburn_url)
print (postburn_url)
titiler_url = "https://titiler.xyz/cog/tiles/WebMercatorQuad/{z}/{x}/{y}"


m = ipyleaflet.Map(center=(lat,lon), 
    zoom = 12, 
    basemap= ipyleaflet.basemaps.OpenStreetMap.Mapnik,
#         basemap= basemaps.Esri.WorldImagery,
    layout=Layout(width='1280px', height='720px'),
    crs = EPSG_MAP)

for polygon in gdf_fire_predict['polygon']:
    multipolygon = ipyleaflet.Polygon(
        locations=polygon,
        color="black",
        fill_color="red",
        fill_opacity = 0.1,
        fill = True,
        weight = 3
    )

    m.add_layer(multipolygon)

tilelayer = ipyleaflet.LocalTileLayer(
    path=f'{titiler_url}@1x?url={preburn_url}')
tilelayer2 = ipyleaflet.LocalTileLayer(
    path=f'{titiler_url}@1x?url={postburn_url}')

right_layer = tilelayer2
left_layer = tilelayer

control = ipyleaflet.SplitMapControl(left_layer=left_layer, right_layer=right_layer)
m.add_control(control)
m.add_control(ipyleaflet.FullScreenControl())

legend = ipyleaflet.LegendControl({"Predicted Fire":"red"},name='' , position="topright")
m.add_control(legend)

m<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11891137473831123, 'fraction_numerical': 0.02512212142358688, 'mean_word_length': 3.602440590879897, 'pattern_counts': {'":': 2, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '534619', 'license': 'Apache-2.0', 'n_tokens_mistral': 2899, 'n_tokens_neox': 2749, 'n_words': 406, 'path': 'fire_prediction_visualisation/fire_overlay_burned_area_sentinel2.ipynb', 'repo_name': 'SatelliteVu/SatelliteVu-AWS-Disaster-Response-Hackathon', 'size': '7165'}