# -*- coding: utf-8 -*- """ Satellite imagery animation widget, which can be used to interactively produce animations for multiple DE Africa products. """ # Import required packages # Force GeoPandas to use Shapely instead of PyGEOS # In a future release, GeoPandas will switch to using Shapely by default. import os os.environ['USE_PYGEOS'] = '0' import fiona import sys import datacube import warnings import matplotlib.pyplot as plt from datacube.utils.geometry import CRS from ipyleaflet import ( WMSLayer, basemaps, basemap_to_tiles, Map, DrawControl, WidgetControl, LayerGroup, LayersControl, GeoData, ) from traitlets import Unicode from ipywidgets import ( GridspecLayout, Button, Layout, HBox, VBox, HTML, Output, ) import json import itertools import numpy as np import geopandas as gpd from io import BytesIO import ipywidgets as widgets import datetime from skimage import exposure from skimage.filters import unsharp_mask from datacube.utils import masking from datacube.utils.geometry import Geometry from datacube.utils.masking import mask_invalid_data import deafrica_tools.app.widgetconstructors as deawidgets from deafrica_tools.dask import create_local_dask_cluster from deafrica_tools.spatial import reverse_geocode from deafrica_tools.datahandling import pan_sharpen_brovey import warnings warnings.filterwarnings("ignore") # WMS params and satellite style bands sat_params = { "Landsat": { "products": ["ls5_sr", "ls7_sr", "ls8_sr", "ls9_sr"], "styles": { "True colour": ("true_colour", ["red", "green", "blue"]), "False colour": ( "false_colour", ["swir_1", "nir", "green"], ), }, }, "Sentinel-2": { "products": ["s2_l2a"], "styles": { "True colour": ("simple_rgb", ["red", "green", "blue"]), "False colour": ( "infrared_green", ["swir_2", "nir_1", "green"], ), }, }, } def make_box_layout(): return Layout( # border='solid 1px black', margin="0px 10px 10px 0px", padding="5px 5px 5px 5px", width="100%", height="100%", ) def create_expanded_button(description, button_style): return Button( description=description, button_style=button_style, layout=Layout(width="auto", height="auto"), ) def update_map_layers(self): """ Updates map to add new DE Africa layers, styles or basemap when selected using menu options. Triggers data reload by resetting load params and output arrays. """ # Clear data load params to trigger data re-load self.timeseries_ds = None self.load_params = None self.query_params = None # Clear all layers and add basemap self.map_layers.clear_layers() self.map_layers.add_layer(self.basemap) def extract_data(self): # Connect to datacube database dc = datacube.Datacube(app="Exporting satellite images") # Configure local dask cluster client = create_local_dask_cluster(return_client=True, display_client=True) # Convert to geopolygon geopolygon = Geometry(geom=self.gdf_drawn.geometry[0], crs=self.gdf_drawn.crs) # Create query. start_date = np.datetime64(self.start_date) end_date = np.datetime64(self.end_date) self.query_params = { "time": (str(start_date), str(end_date)), "geopolygon": geopolygon, } # Find matching datasets dss = [ dc.find_datasets(product=i, **self.query_params) for i in sat_params[self.dealayer]["products"] ] dss = list(itertools.chain.from_iterable(dss)) # If data is found if len(dss) > 0: # Get CRS crs = str(dss[0].crs) self.load_params = { "measurements": sat_params[self.dealayer]["styles"][self.style][1], "resolution": (-self.resolution, self.resolution), "output_crs": crs, "group_by": "solar_day", "dask_chunks": {"time": 1, "x": 2048, "y": 2048}, "resampling": {"*": "cubic", "oa_fmask": "nearest", "fmask": "nearest"}, } # Load data from deafrica_tools.datahandling import load_ard timeseries_ds = load_ard( dc=dc, products=sat_params[self.dealayer]["products"], min_gooddata=1.0 - (self.max_cloud_cover / 100), ls7_slc_off=False, mask_pixel_quality=self.cloud_mask, **self.load_params, **self.query_params, ) # Set invalid nodata pixels to NaN timeseries_ds = mask_invalid_data(timeseries_ds) # Else if no data is returned, return None else: timeseries_ds = None # Close down the dask client client.close() return timeseries_ds.compute() def plot_data(self, fname): # Data to plot to_plot = self.timeseries_ds # If rolling median specified if self.rolling_median: with self.status_info: print( f"\nApplying rolling median ({self.rolling_median_window} timesteps window)" ) to_plot = to_plot.rolling( time=int(self.rolling_median_window), center=True, min_periods=1 ).median() # If resampling freq specified if self.resample_freq: with self.status_info: print(f"\nResampling data to {self.resample_freq} frequency") to_plot = to_plot.resample(time=self.resample_freq).median() # Raise by power to dampen bright features and enhance dark. # Raise vmin and vmax by same amount to ensure proper stretch if self.power < 1.0: with self.status_info: print(f"\nApplying power transformation ({self.power})") to_plot = to_plot ** self.power # Apply unsharp masking to enhance overall dynamic range, # and improve fine scale detail if self.unsharp_mask: with self.status_info: print( f"\nApplying unsharp masking with {self.unsharp_mask_radius} " f"radius and {self.unsharp_mask_amount} amount" ) from skimage.exposure import rescale_intensity funcs_list = [ rescale_intensity, lambda x: unsharp_mask( x, radius=self.unsharp_mask_radius, amount=self.unsharp_mask_amount ), ] else: funcs_list = None from deafrica_tools.plotting import xr_animation xr_animation( output_path=fname, ds=to_plot.dropna(dim="time", how="all"), show_text="", bands=sat_params[self.dealayer]["styles"][self.style][1], interval=self.interval, width_pixels=self.width, show_gdf=deacoastlines_overlay(to_plot) if self.deacoastlines else None, gdf_kwargs={"linewidth": 3}, percentile_stretch=(self.vmin, self.vmax), image_proc_funcs=funcs_list, show_date="%Y" if self.resample_freq == "1Y" else "%b %Y", annotation_kwargs={"fontsize": 75}, ) # Add plot preview below map and finish plt.show() with self.status_info: print(f"\nImage successfully exported to:\n{fname}.") def deacoastlines_overlay(ds): import geopandas as gpd import pandas as pd import matplotlib from shapely.geometry import box, Point from deafrica_tools.coastal import get_coastlines # Get bounding box of data xmin, ymin, xmax, ymax = ds.geobox.geographic_extent.boundingbox bounds = [xmin, ymin, xmax, ymax] # Load data deacl_gdf = get_coastlines(bbox=bounds) # Clip to extent of satellite data bbox = gpd.GeoDataFrame(geometry=[ds.geobox.extent.geom], crs=ds.geobox.crs) deacl_gdf = gpd.overlay(deacl_gdf, bbox.to_crs(deacl_gdf.crs)) deacl_gdf = deacl_gdf.dissolve("year") # values("year", ascending=True) # Apply colours norm = matplotlib.colors.Normalize(vmin=0, vmax=len(deacl_gdf.index)) cmap = matplotlib.cm.get_cmap("inferno") rgba = cmap(norm(deacl_gdf.reset_index().index)) deacl_gdf["color"] = list(rgba) deacl_gdf["start_time"] = pd.to_datetime(deacl_gdf.index) + pd.DateOffset(months=0) deacl_gdf = deacl_gdf.sort_index() if len(deacl_gdf.index) > 0: return deacl_gdf else: return None class animation_app(HBox): def __init__(self): super().__init__() ###################### # INITIAL ATTRIBUTES # ###################### # Basemap self.basemap_list = [ ("ESRI World Imagery", basemap_to_tiles(basemaps.Esri.WorldImagery)), ("Open Street Map", basemap_to_tiles(basemaps.OpenStreetMap.Mapnik)), ] self.basemap = self.basemap_list[0][1] # Satellite data end_date = datetime.datetime.today() start_date = datetime.datetime( year=end_date.year - 3, month=end_date.month, day=end_date.day ) self.start_date = start_date.strftime("%Y-%m-%d") self.end_date = end_date.strftime("%Y-%m-%d") self.dealayer_list = [ ("Landsat", "Landsat"), ("Sentinel-2", "Sentinel-2"), ] self.dealayer = self.dealayer_list[0][1] # Styles self.styles_list = ["True colour", "False colour"] self.style = self.styles_list[0] # Analysis params self.resolution = 30 self.vmin = 0.01 self.vmax = 0.99 self.power = 1.0 self.output_list = [("MP4", "mp4"), ("GIF", "gif")] self.output_format = self.output_list[0][1] self.rolling_median = False self.rolling_median_window = 20 self.unsharp_mask = False self.unsharp_mask_radius = 20 self.unsharp_mask_amount = 0.3 self.max_size = False self.width = 900 self.interval = 100 self.cloud_mask = False self.max_cloud_cover = 20 self.resample_list = [ ("None", False), ("Monthly", "1M"), ("Quarterly", "Q-DEC"), ("Yearly", "1Y"), ] self.resample_freq = self.resample_list[0][1] self.deacoastlines = False # Drawing params self.target = None self.action = None self.gdf_drawn = None # Data load params self.timeseries_ds = None self.load_params = None self.query_params = None ################## # HEADER FOR APP # ################## # Create the Header widget header_title_text = ( "
Select the desired satellite data, imagery date range " "and image style, then zoom in and draw a rectangle to " "select an area export as a satellite imagery time-series " "animation.
" ) self.header = deawidgets.create_html(f"{header_title_text}{instruction_text}") self.header.layout = make_box_layout() ##################################### # HANDLER FUNCTION FOR DRAW CONTROL # ##################################### # Define the action to take once something is drawn on the map def update_geojson(target, action, geo_json): # Get data from action self.action = action # Clear data load params to trigger data re-load self.timeseries_ds = None self.load_params = None self.query_params = None # Convert data to geopandas json_data = json.dumps(geo_json) binary_data = json_data.encode() io = BytesIO(binary_data) io.seek(0) gdf = gpd.read_file(io) gdf.crs = "EPSG:4326" # Convert to WGS 84 / NSIDC EASE-Grid 2.0 Global and compute area gdf_drawn_nsidc = gdf.copy().to_crs("EPSG:6933") m2_per_ha = 10000 area = gdf_drawn_nsidc.area.values[0] / m2_per_ha polyarea_label = "Total area of satellite data to extract" polyarea_text = f"{polyarea_label}: {area:.2f} ha" # Test area size if self.max_size: confirmation_text = ( ' ' "(Overriding maximum size limit; use with caution as may lead to memory issues)" ) self.header.value = ( header_title_text + instruction_text + polyarea_text + confirmation_text ) self.gdf_drawn = gdf elif area <= 50000: confirmation_text = ( ' ' "(Area to extract falls within " "recommended 50000 ha limit)" ) self.header.value = ( header_title_text + instruction_text + polyarea_text + confirmation_text ) self.gdf_drawn = gdf else: warning_text = ( ' ' "(Area to extract is too large, " "please select an area less than 50000 " "ha)" ) self.header.value = ( header_title_text + instruction_text + polyarea_text + warning_text ) self.gdf_drawn = None ########################### # WIDGETS FOR APP OUTPUTS # ########################### self.status_info = Output(layout=make_box_layout()) self.output_plot = Output(layout=make_box_layout()) ######################################### # MAP WIDGET, DRAWING TOOLS, WMS LAYERS # ######################################### # Create drawing tools desired_drawtools = ["rectangle"] draw_control = deawidgets.create_drawcontrol(desired_drawtools) # Begin by displaying an empty layer group, and update the group with desired WMS on interaction. self.map_layers = LayerGroup(layers=()) self.map_layers.name = "Map Overlays" # Create map widget self.m = deawidgets.create_map(map_center=(5.65, 26.17), zoom_level=13) self.m.layout = make_box_layout() # Add tools to map widget self.m.add_control(draw_control) self.m.add_layer(self.map_layers) # Update all maps to starting defaults update_map_layers(self) ############################ # WIDGETS FOR APP CONTROLS # ############################ # Create parameter widgets dropdown_basemap = deawidgets.create_dropdown( self.basemap_list, self.basemap_list[0][1] ) dropdown_dealayer = deawidgets.create_dropdown( self.dealayer_list, self.dealayer_list[0][1] ) dropdown_output = deawidgets.create_dropdown( self.output_list, self.output_list[0][1] ) date_picker_start = deawidgets.create_datepicker( value=start_date, ) date_picker_end = deawidgets.create_datepicker( value=end_date, ) dropdown_styles = deawidgets.create_dropdown( self.styles_list, self.styles_list[0] ) slider_percentile = widgets.FloatRangeSlider( value=[0.01, 0.99], min=0, max=1, step=0.001, description="", layout={"width": "85%"}, ) run_button = create_expanded_button("Generate animation", "info") floatslider_max_cloud_cover = widgets.IntSlider( value=20, min=0, max=100, step=1, description="", layout={"width": "85%"}, ) checkbox_rolling_median = deawidgets.create_checkbox( self.rolling_median, "Apply rolling median