Source code for dea_tools.app.deacoastlines

"""
Digital Earth Australia Coastline widget, which can be used to 
interactively extract shoreline data using transects.
"""

# Import required packages
import fiona
import os
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 geopandas as gpd
from io import BytesIO
import ipywidgets as widgets

import dea_tools.app.widgetconstructors as deawidgets
from dea_tools.coastal import get_coastlines, transect_distances


WMS_ADDRESS = "https://geoserver.dea.ga.gov.au/geoserver/wms"


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"),
    )


[docs] class transect_app(HBox): def __init__(self): super().__init__() ###################### # INITIAL ATTRIBUTES # ###################### self.output_name = "example_output" self.export_csv = False self.export_plot = False self.product_list = [ ("ESRI World Imagery", "none"), ("Open Street Map", "open_street_map"), ] self.product = self.product_list[0][1] self.mode_list = [('Distance', 'distance'), ('Width', 'width')] self.mode = self.mode_list[0][1] self.target = None self.action = None self.gdf_drawn = None self.gdf_uploaded = None ################## # HEADER FOR APP # ################## # Create the Header widget header_title_text = "<h3>Digital Earth Australia Coastlines shoreline transect extraction</h3>" instruction_text = "Select parameters and draw a transect on the map to extract shoreline data. <b>In distance mode</b>, draw a transect line starting from land that crosses multiple shorelines. <br><b>In width mode</b>, draw a transect line that intersects shorelines at least twice. Alternatively, <b>upload an vector file</b> to extract shoreline data for multiple existing transects." self.header = deawidgets.create_html( f"{header_title_text}<p>{instruction_text}</p>") 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): # Remove previously uploaded data if present self.gdf_uploaded = None fileupload_transects._counter = 0 # Get data from action self.action = action # 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 Albers and compute area gdf_drawn_albers = gdf.copy().to_crs("EPSG:3577") m2_per_km2 = 10**6 area = gdf_drawn_albers.envelope.area.values[0] / m2_per_km2 polyarea_label = 'Total area of DEA Coastlines data to extract' polyarea_text = f"<b>{polyarea_label}</b>: {area:.2f} km<sup>2</sup>" # Test area size if area <= 50000: confirmation_text = '<span style="color: #33cc33"> <b>(Area to extract falls within recommended limit; click "Extract shoreline data" to continue)</b></span>' self.header.value = header_title_text + polyarea_text + confirmation_text self.gdf_drawn = gdf else: warning_text = '<span style="color: #ff5050"> <b>(Area to extract is too large, please select a smaller transect)</b></span>' self.header.value = header_title_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 = ['polyline'] draw_control = deawidgets.create_drawcontrol(desired_drawtools) # Load DEACoastLines WMS deacl_url = WMS_ADDRESS deacl_layer = 'dea:DEACoastlines' deacoastlines = WMSLayer( url=deacl_url, layers=deacl_layer, format='image/png', transparent=True, attribution='DEA Coastlines © 2020 Geoscience Australia') # Begin by displaying an empty layer group, and update the group with desired WMS on interaction. self.map_layers = LayerGroup(layers=(deacoastlines,)) self.map_layers.name = 'Map Overlays' # Create map widget self.m = deawidgets.create_map(map_center=(-28, 135), zoom_level=4, basemap=basemaps.Esri.WorldImagery) self.m.layout = make_box_layout() # Add tools to map widget self.m.add_control(draw_control) self.m.add_layer(self.map_layers) # Store current basemap for future use self.basemap = self.m.basemap ############################ # WIDGETS FOR APP CONTROLS # ############################ # Create parameter widgets text_output_name = deawidgets.create_inputtext(self.output_name, self.output_name) checkbox_csv = deawidgets.create_checkbox(self.export_csv, 'Distance table (.csv)') checkbox_plot = deawidgets.create_checkbox(self.export_plot, 'Figure (.png)') deaoverlay_dropdown = deawidgets.create_dropdown( self.product_list, self.product_list[0][1]) mode_dropdown = deawidgets.create_dropdown(self.mode_list, self.mode_list[0][1]) run_button = create_expanded_button("Extract shoreline data", "info") fileupload_transects = widgets.FileUpload(accept='', multiple=True) #################################### # UPDATE FUNCTIONS FOR EACH WIDGET # #################################### # Run update functions whenever various widgets are changed. text_output_name.observe(self.update_text_output_name, "value") checkbox_csv.observe(self.update_checkbox_csv, "value") checkbox_plot.observe(self.update_checkbox_plot, "value") deaoverlay_dropdown.observe(self.update_deaoverlay, "value") mode_dropdown.observe(self.update_mode, "value") run_button.on_click(self.run_app) draw_control.on_draw(update_geojson) fileupload_transects.observe(self.update_fileupload_transects, "value") ################################## # COLLECTION OF ALL APP CONTROLS # ################################## parameter_selection = VBox([ HTML("<b>Output name:</b>"), text_output_name, HTML( '<b>Transect extraction mode:</b><br><img src="https://i.imgur.com/9fdTH9C.png">' ), mode_dropdown, HTML("<b></br>Output files:</b>"), checkbox_plot, checkbox_csv, HTML( "</br><i><b>Advanced</b></br>Upload a GeoJSON or ESRI " "Shapefile (<5 mb) containing one or more transect lines.</i>"), fileupload_transects ]) map_selection = VBox([ HTML("</br><b>Map overlay:</b>"), deaoverlay_dropdown, ]) parameter_selection.layout = make_box_layout() map_selection.layout = make_box_layout() ############################### # SPECIFICATION OF APP LAYOUT # ############################### # 0 1 2 3 4 5 6 7 8 9 # --------------------------------------------- # 0 | Header | Map sel. | # --------------------------------------------- # 1 | Params | | # 2 | | | # 3 | | | # 4 | | Map | # 5 | | | # ---------- | # 6 | Run | | # --------------------------------------------- # 7 | Status info | # --------------------------------------------- # 8 | | # 9 | Output/figure | # 10 | | # 11 | ------------------------------------------| # Create the layout #[rowspan, colspan] grid = GridspecLayout(12, 10, height="1350px", width="auto") # Header and controls grid[0, :8] = self.header grid[0, 8:] = map_selection grid[1:6, 0:2] = parameter_selection grid[6, 0:2] = run_button # Status info, map and plot grid[1:7, 2:] = self.m # map grid[7:8, :] = self.status_info grid[8:, :] = self.output_plot # Display using HBox children attribute self.children = [grid] ###################################### # DEFINITION OF ALL UPDATE FUNCTIONS # ###################################### # Set the output csv def update_fileupload_transects(self, change): # Clear any drawn data if present self.gdf_drawn = None # Temporary compatibility fix for ipywidget > 8.0 # TODO: Update code to use new fileupload API documented here: # https://ipywidgets.readthedocs.io/en/latest/user_migration_guides.html#fileupload uploaded_data = {f["name"]: {"content": f.content.tobytes()} for f in change.new} # Save to file for uploaded_filename in uploaded_data.keys(): with open(uploaded_filename, "wb") as output_file: content = uploaded_data[uploaded_filename]["content"] output_file.write(content) with self.status_info: try: print('Loading vector data...', end='\r') valid_files = [ file for file in uploaded_data.keys() if file.lower().endswith(('.shp', '.geojson')) ] valid_file = valid_files[0] transect_gdf = (gpd.read_file(valid_file).to_crs( "EPSG:4326").explode(index_parts=True).reset_index(drop=True)) # Use ID column if it exists if 'id' in transect_gdf: transect_gdf = transect_gdf.set_index('id') print(f"Uploaded '{valid_file}'; automatically labelling " "transects using column 'id'.") else: print( f"Uploaded '{valid_file}'; no 'id' column detected, " f"labelling transects from 0 to {len(transect_gdf.index) - 1}." ) # Create a geodata geodata = GeoData(geo_dataframe=transect_gdf, style={ 'color': 'black', 'weight': 3 }) # Add to map xmin, ymin, xmax, ymax = transect_gdf.total_bounds self.m.fit_bounds([[ymin, xmin], [ymax, xmax]]) self.m.add_layer(geodata) # If completed, add to attribute self.gdf_uploaded = transect_gdf except IndexError: print( "Cannot read uploaded files. Please ensure that data is " "in either GeoJSON or ESRI Shapefile format.", end='\r') self.gdf_uploaded = None except fiona.errors.DriverError: print( "Shapefile is invalid. Please ensure that all shapefile " "components (e.g. .shp, .shx, .dbf, .prj) are uploaded.", end='\r') self.gdf_uploaded = None # Set output name def update_text_output_name(self, change): self.output_name = change.new # Output CSV def update_checkbox_csv(self, change): self.export_csv = change.new # Output plot def update_checkbox_plot(self, change): self.export_plot = change.new # Set mode def update_mode(self, change): self.mode = change.new # Update product def update_deaoverlay(self, change): self.product = change.new # Load DEACoastLines WMS deacl_url = WMS_ADDRESS deacl_layer = "dea:DEACoastlines" deacoastlines = WMSLayer( url=deacl_url, layers=deacl_layer, format="image/png", transparent=True, attribution="DEA Coastlines © 2023 Geoscience Australia") if self.product == "none": self.map_layers.clear_layers() self.map_layers.add_layer(deacoastlines) elif self.product == "open_street_map": self.map_layers.clear_layers() layer = basemap_to_tiles(basemaps.OpenStreetMap.Mapnik) self.map_layers.add_layer(layer) self.map_layers.add_layer(deacoastlines) def run_app(self, change): # Clear progress bar and output areas before running self.status_info.clear_output() self.output_plot.clear_output() # Run DEA Coastlines analysis with self.status_info: warnings.filterwarnings("ignore") # Load transects from either map or uploaded files if self.gdf_uploaded is not None: transect_gdf = self.gdf_uploaded run_text = 'uploaded file' elif self.gdf_drawn is not None: transect_gdf = self.gdf_drawn transect_gdf.index = [self.output_name] run_text = 'selected transect' else: print(f'No transect drawn or uploaded. Please select a transect on the map, or upload a GeoJSON or ESRI Shapefile.', end='\r') transect_gdf = None # If valid data was returned, load DEA Coastlines data if transect_gdf is not None: # Load Coastlines data from WFS deacl_gdf = get_coastlines(bbox=transect_gdf) # Test that data was correctly returned if len(deacl_gdf.index) > 0: # Dissolve by year to remove duplicates, then sort by date deacl_gdf = deacl_gdf.dissolve(by='year', as_index=False) deacl_gdf['year'] = deacl_gdf.year.astype(int) deacl_gdf = deacl_gdf.sort_values('year') deacl_gdf = deacl_gdf.set_index('year') else: print( "No annual shoreline data was found near the " "supplied transect. Please draw or select a new " "transect.", end='\r') deacl_gdf = None # If valid DEA Coastlines data returned, calculate distances if deacl_gdf is not None: print(f'Analysing transect distances using "{self.mode}" mode...', end='\r') dist_df = transect_distances( transect_gdf.to_crs('EPSG:3577'), deacl_gdf, mode=self.mode) # If valid data was produced: if dist_df.any(axis=None): # Successful output print(f'DEA Coastlines data successfully extracted for {run_text}.') # Export distance data if self.export_csv: # Create folder if required and set path out_dir = 'deacoastlines_outputs' os.makedirs(out_dir, exist_ok=True) csv_filename = f"{out_dir}/{self.output_name}.csv" # Export to file dist_df.to_csv(csv_filename, index_label="Transect") print(f'Distance data exported to "{csv_filename}".') # Generate plot with self.output_plot: fig, ax = plt.subplots(constrained_layout=True, figsize=(15, 5.5)) dist_df.T.plot(ax=ax, linewidth=3) ax.legend(frameon=False, ncol=3, title='Transect') ax.set_title(f"Digital Earth Australia Coastlines transect extraction - {self.output_name}") ax.set_ylabel(f"Along-transect {self.mode} (m)") ax.set_xlim(dist_df.T.index[0], dist_df.T.index[-1]) # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') plt.show() # Export plot with self.status_info: if self.export_plot: # Create folder if required and set path out_dir = 'deacoastlines_outputs' os.makedirs(out_dir, exist_ok=True) figure_filename = f"{out_dir}/{self.output_name}.png" # Export to file fig.savefig(figure_filename) print(f'Figure exported to "{figure_filename}".') else: print( "No valid shoreline data intersects with the " "supplied transect. This can occur if:\n\n" " - the transect does not intersect with any shorelines\n" " - the transect intersects with shorelines more than once in 'distance' mode\n" " - the transect intersects with shorelines only once in 'width' mode\n\n" "Please draw or upload a new transect.", end='\r')