diff --git a/Real_world_examples/Forest_monitoring.ipynb b/Real_world_examples/Forest_monitoring.ipynb index 21c03504..630585d3 100644 --- a/Real_world_examples/Forest_monitoring.ipynb +++ b/Real_world_examples/Forest_monitoring.ipynb @@ -26,14 +26,14 @@ "source": [ "## Background\n", "\n", - "Forests worldwide are in a state of flux, with accelerating losses in some regions and gains in others [(Hansen et al., 2013)](https://doi.org/10.1126/science.1244693). The Global Forest Change 2000-2021 dataset characterizes the global forest extent and change from 2000 to 2021. \n", + "Forests worldwide are in a state of flux, with accelerating losses in some regions and gains in others [(Hansen et al., 2013)](https://doi.org/10.1126/science.1244693). The Global Forest Change 2000-2023 dataset characterizes the global forest extent and change from 2000 to 2023. \n", "\n", "For this dataset: \n", "- Forest loss is defined as stand-replacement disturbance, or a change from a forest to non-forest state\n", "- Forest gain is defined as the inverse of as the inverse of loss, or a non-forest to forest change entirely within the study period\n", "- Tree cover is defined as canopy closure for all vegetation taller than 5m in height.\n", "\n", - "The Year of gross forest cover loss event (`lossyear`) layer shows the forest loss during the period 2000 to 2021. Forest loss is encoded as either 0 (no loss) or else a value in the range 1-20, representing loss detected primarily in the year 2001-2021, respectively.\n", + "The Year of gross forest cover loss event (`lossyear`) layer shows the forest loss during the period 2000 to 2023. Forest loss is encoded as either 0 (no loss) or else a value in the range 1-20, representing loss detected primarily in the year 2001-2023, respectively.\n", "\n", "\n", "The Tree canopy cover for year 2000 (`treecover2000`) layer shows the tree cover in the year 2000. \n", @@ -111,10 +111,10 @@ "\n", "| Global Forest Change layers available | |\n", "| --- | --- |\n", - "| **Year of gross forest cover loss event** | Shows forest loss during the period 2000–2021 | \n", + "| **Year of gross forest cover loss event** | Shows forest loss during the period 2000–2023 | \n", "| **Global forest cover gain 2000–2012** | Shows forest gain during the period 2000-2012 |\n", "| **Tree canopy cover for year 2000** | Show the tree canopy cover for the year 2000 |\n", - "| **All layers** | Show the forest loss during the period 2000–2021, forest gain during the period 2000-2012 and the tree canopy cover for the year 2000 |\n", + "| **All layers** | Show the forest loss during the period 2000–2023, forest gain during the period 2000-2012 and the tree canopy cover for the year 2000 |\n", "\n", "Use the `Forest Cover Loss Time Range` slider to set the time range for which to load the **Year of gross forest cover loss event** layer. Tick the `Override maximum size limit` box to override the app's default 500 square kilometres area limit. This can be used to load larger areas of imagery, but should be used with caution as it can lead to memory issues or crashes.\n", "\n", @@ -135,7 +135,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "9e01e78f8dd84328a90629807a49409d", + "model_id": "d0223d060cdd4162b17fed65c8f9112c", "version_major": 2, "version_minor": 0 }, @@ -222,7 +222,7 @@ { "data": { "text/plain": [ - "'2023-08-14'" + "'2024-09-05'" ] }, "execution_count": 4, diff --git a/Tools/deafrica_tools/__init__.py b/Tools/deafrica_tools/__init__.py index 149c0ec7..4c3d12d0 100644 --- a/Tools/deafrica_tools/__init__.py +++ b/Tools/deafrica_tools/__init__.py @@ -1,6 +1,6 @@ __locales__ = __path__[0] + "/locales" -__version__ = "2.4.8" +__version__ = "2.4.9" def set_lang(lang=None): diff --git a/Tools/deafrica_tools/app/changefilmstrips.py b/Tools/deafrica_tools/app/changefilmstrips.py index 31fb2d87..563102d3 100644 --- a/Tools/deafrica_tools/app/changefilmstrips.py +++ b/Tools/deafrica_tools/app/changefilmstrips.py @@ -3,37 +3,31 @@ inside the Real_world_examples folder. """ -# Load modules -import os -import dask -import datacube import warnings + +import datacube +import matplotlib.pyplot as plt import numpy as np import pandas as pd import xarray as xr -import matplotlib.pyplot as plt +from datacube.utils.geometry import CRS, assign_crs +from ipyleaflet import basemap_to_tiles, basemaps from odc.algo import geomedian_with_mads from odc.ui import select_on_a_map -from dask.utils import parse_bytes -from datacube.utils.geometry import CRS, assign_crs -from datacube.utils.rio import configure_s3_access -from datacube.utils.dask import start_local_dask -from ipyleaflet import basemaps, basemap_to_tiles -# Load utility functions -from deafrica_tools.datahandling import load_ard, mostcommon_crs from deafrica_tools.dask import create_local_dask_cluster +from deafrica_tools.datahandling import load_ard, mostcommon_crs def run_filmstrip_app( - output_name, - time_range, - time_step, - tide_range=(0.0, 1.0), - resolution=(-30, 30), - max_cloud=0.5, - ls7_slc_off=False, - size_limit=10000, + output_name: str, + time_range: tuple, + time_step: dict[str, int], + tide_range: tuple[float] = (0.0, 1.0), + resolution: tuple[int] = (-30, 30), + max_cloud: float = 0.5, + ls7_slc_off: bool = False, + size_limit: int = 10000, ): """ An interactive app that allows the user to select a region from a @@ -91,7 +85,7 @@ def run_filmstrip_app( size_limit : int, optional An optional integer (in hectares) specifying the size limit for the data query. Queries larger than this size will receive - a warning that he data query is too large (and may + a warning that the data query is too large (and may therefore result in memory errors). @@ -117,10 +111,7 @@ def run_filmstrip_app( # Plot interactive map to select area basemap = basemap_to_tiles(basemaps.Esri.WorldImagery) - geopolygon = select_on_a_map(height="600px", - layers=(basemap,), - center=centre_coords, - zoom=14) + geopolygon = select_on_a_map(height="600px", layers=(basemap,), center=centre_coords, zoom=14) # Set centre coords based on most recent selection to re-focus # subsequent data selections @@ -129,12 +120,14 @@ def run_filmstrip_app( # Test size of selected area msq_per_hectare = 10000 area = geopolygon.to_crs(crs=CRS("epsg:6933")).area / msq_per_hectare - radius = np.round(np.sqrt(size_limit), 1) + radius = np.round(np.sqrt(size_limit), 1) # noqa F841 if area > size_limit: - print(f"Warning: Your selected area is {area:.00f} hectares. " - f"Please select an area of less than {size_limit} hectares." - f"\nTo select a smaller area, re-run the cell " - f"above and draw a new polygon.") + print( + f"Warning: Your selected area is {area:.00f} hectares. " + f"Please select an area of less than {size_limit} hectares." + f"\nTo select a smaller area, re-run the cell " + f"above and draw a new polygon." + ) else: @@ -147,12 +140,9 @@ def run_filmstrip_app( client = create_local_dask_cluster(return_client=True) # Obtain native CRS - crs = mostcommon_crs(dc=dc, - product="ls8_sr", - query={ - "time": "2014", - "geopolygon": geopolygon - }) + crs = mostcommon_crs( + dc=dc, product="ls8_sr", query={"time": "2014", "geopolygon": geopolygon} + ) # Create query based on time range, area selected, custom params query = { @@ -160,10 +150,7 @@ def run_filmstrip_app( "geopolygon": geopolygon, "output_crs": crs, "resolution": resolution, - "dask_chunks": { - "x": 3000, - "y": 3000 - }, + "dask_chunks": {"x": 3000, "y": 3000}, "align": (resolution[1] / 2.0, resolution[1] / 2.0), } @@ -179,39 +166,46 @@ def run_filmstrip_app( ) # Optionally calculate tides for each timestep in the satellite - # dataset and drop any observations out side this range + # dataset and drop any observations outside this range if tide_range != (0.0, 1.0): from deafrica_tools.coastal import tidal_tag + ds = tidal_tag(ds=ds, tidepost_lat=None, tidepost_lon=None) min_tide, max_tide = ds.tide_height.quantile(tide_range).values - ds = ds.sel(time=(ds.tide_height >= min_tide) & - (ds.tide_height <= max_tide)) + ds = ds.sel(time=(ds.tide_height >= min_tide) & (ds.tide_height <= max_tide)) ds = ds.drop("tide_height") - print(f" Keeping {len(ds.time)} observations with tides " - f"between {min_tide:.2f} and {max_tide:.2f} m") + print( + f" Keeping {len(ds.time)} observations with tides " + f"between {min_tide:.2f} and {max_tide:.2f} m" + ) # Create time step ranges to generate filmstrips from - bins_dt = pd.date_range(start=time_range[0], - end=time_range[1], - freq=pd.DateOffset(**time_step)) + bins_dt = pd.date_range( + start=time_range[0], end=time_range[1], freq=pd.DateOffset(**time_step) + ) # Bin all satellite observations by timestep. If some observations # fall outside the upper bin, label these with the highest bin labels = bins_dt.astype("str") - time_steps = (pd.cut(ds.time.values, bins_dt, - labels=labels[:-1]).add_categories( - labels[-1]).fillna(labels[-1])) + time_steps = ( + pd.cut(ds.time.values, bins_dt, labels=labels[:-1]) + .add_categories(labels[-1]) + .fillna(labels[-1]) + ) - time_steps_var = xr.DataArray(time_steps, [("time", ds.time.values)], - name="timestep") + time_steps_var = xr.DataArray(time_steps, [("time", ds.time.values)], name="timestep") # Resample data temporally into time steps, and compute geomedians - ds_geomedian = (ds.groupby(time_steps_var).apply( + ds_geomedian = ds.groupby(time_steps_var).apply( lambda ds_subset: geomedian_with_mads( - ds_subset, compute_mads=False, compute_count=False))) + ds_subset, compute_mads=False, compute_count=False + ) + ) - print("\nGenerating geomedian composites and plotting " - "filmstrips... (click the Dashboard link above for status)") + print( + "\nGenerating geomedian composites and plotting " + "filmstrips... (click the Dashboard link above for status)" + ) ds_geomedian = ds_geomedian.compute() # Reset CRS that is lost during geomedian compositing @@ -230,16 +224,14 @@ def run_filmstrip_app( # and aspect ratio n_obs = output_array.sizes["timestep"] ratio = output_array.sizes["x"] / output_array.sizes["y"] - fig, axes = plt.subplots(1, - n_obs + 1, - figsize=(5 * ratio * (n_obs + 1), 5)) + fig, axes = plt.subplots(1, n_obs + 1, figsize=(5 * ratio * (n_obs + 1), 5)) fig.subplots_adjust(wspace=0.05, hspace=0.05) # Add timesteps to the plot, set aspect to equal to preserve shape for i, ax_i in enumerate(axes.flatten()[:n_obs]): - output_array.isel(timestep=i).plot.imshow(ax=ax_i, - vmin=percentiles[0], - vmax=percentiles[1]) + output_array.isel(timestep=i).plot.imshow( + ax=ax_i, vmin=percentiles[0], vmax=percentiles[1] + ) ax_i.get_xaxis().set_visible(False) ax_i.get_yaxis().set_visible(False) ax_i.set_aspect("equal") @@ -248,11 +240,12 @@ def run_filmstrip_app( # by first taking the log of the array (so change in dark areas # can be identified), then computing standard deviation between # all timesteps - (np.log(output_array).std(dim=["timestep"]).mean( - dim="variable").plot.imshow(ax=axes.flatten()[-1], - robust=True, - cmap="magma", - add_colorbar=False)) + ( + np.log(output_array) + .std(dim=["timestep"]) + .mean(dim="variable") + .plot.imshow(ax=axes.flatten()[-1], robust=True, cmap="magma", add_colorbar=False) + ) axes.flatten()[-1].get_xaxis().set_visible(False) axes.flatten()[-1].get_yaxis().set_visible(False) axes.flatten()[-1].set_aspect("equal") diff --git a/Tools/deafrica_tools/app/crophealth.py b/Tools/deafrica_tools/app/crophealth.py index f6fd091f..675d6565 100644 --- a/Tools/deafrica_tools/app/crophealth.py +++ b/Tools/deafrica_tools/app/crophealth.py @@ -1,49 +1,44 @@ # crophealth.py -''' +""" Functions for loading and interacting with data in the crop health notebook, inside the Real_world_examples folder. -''' +""" # Load modules # 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' - -from ipyleaflet import ( - Map, - GeoJSON, - DrawControl, - basemaps -) + +os.environ["USE_PYGEOS"] = "0" + import datetime as dt +import json +import warnings +from io import BytesIO + import datacube -from osgeo import ogr +import geopandas as gpd +import ipywidgets as widgets import matplotlib as mpl import matplotlib.pyplot as plt -import rasterio -from rasterio.features import geometry_mask import xarray as xr +from ipyleaflet import DrawControl, GeoJSON, Map, basemaps from IPython.display import display -import warnings -import ipywidgets as widgets -import json -import geopandas as gpd -from io import BytesIO +from osgeo import ogr -# Load utility functions +from deafrica_tools.bandindices import calculate_indices from deafrica_tools.datahandling import load_ard from deafrica_tools.spatial import xr_rasterize -from deafrica_tools.bandindices import calculate_indices -def load_crophealth_data(lat, lon, buffer, date): +def load_crophealth_data(lat: float, lon: float, buffer: float, date: str) -> xr.Dataset: """ Loads Sentinel-2 analysis-ready data (ARD) product for the crop health case-study area over the last two years. Last modified: April 2020 - + Parameters ---------- lat: float @@ -51,7 +46,7 @@ def load_crophealth_data(lat, lon, buffer, date): lon: float The central longitude to analyse buffer: - The number of square degrees to load around the central latitude and longitude. + The number of square degrees to load around the central latitude and longitude. For reasonable loading times, set this as `0.1` or lower. date: The most recent date to show data for. @@ -59,17 +54,17 @@ def load_crophealth_data(lat, lon, buffer, date): Returns ---------- - ds: xarray.Dataset + ds: xarray.Dataset data set containing combined, masked data Masked values are set to 'nan' """ - + # Suppress warnings - warnings.filterwarnings('ignore') + warnings.filterwarnings("ignore") # Initialise the data cube. 'app' argument is used to identify this app - dc = datacube.Datacube(app='Crophealth-app') - + dc = datacube.Datacube(app="Crophealth-app") + # Define area to load latitude = (lat - buffer, lat + buffer) longitude = (lon - buffer, lon + buffer) @@ -84,20 +79,14 @@ def load_crophealth_data(lat, lon, buffer, date): # Construct the data cube query products = ["s2_l2a"] - + query = { - 'x': longitude, - 'y': latitude, - 'time': time, - 'measurements': [ - 'red', - 'green', - 'blue', - 'nir', - 'swir_2' - ], - 'output_crs': 'EPSG:6933', - 'resolution': (-20, 20) + "x": longitude, + "y": latitude, + "time": time, + "measurements": ["red", "green", "blue", "nir", "swir_2"], + "output_crs": "EPSG:6933", + "resolution": (-20, 20), } # Load the data and mask out bad quality pixels @@ -106,22 +95,22 @@ def load_crophealth_data(lat, lon, buffer, date): # Calculate the normalised difference vegetation index (NDVI) across # all pixels for each image. # This is stored as an attribute of the data - ds = calculate_indices(ds, index='NDVI', satellite_mission='s2') + ds = calculate_indices(ds, index="NDVI", satellite_mission="s2") # Return the data - return(ds) + return ds -def run_crophealth_app(ds, lat, lon, buffer): +def run_crophealth_app(ds: xr.Dataset, lat: float, lon: float, buffer: float): """ Plots an interactive map of the crop health case-study area and allows the user to draw polygons. This returns a plot of the average NDVI value in the polygon area. Last modified: January 2020 - + Parameters ---------- - ds: xarray.Dataset + ds: xarray.Dataset data set containing combined, masked data Masked values are set to 'nan' lat: float @@ -129,17 +118,17 @@ def run_crophealth_app(ds, lat, lon, buffer): lon: float The central longitude corresponding to the area of loaded ds buffer: - The number of square degrees to load around the central latitude and longitude. + The number of square degrees to load around the central latitude and longitude. For reasonable loading times, set this as `0.1` or lower. """ - + # Suppress warnings - warnings.filterwarnings('ignore') + warnings.filterwarnings("ignore") # Update plotting functionality through rcParams - mpl.rcParams.update({'figure.autolayout': True}) - - # Define polygon bounds + mpl.rcParams.update({"figure.autolayout": True}) + + # Define polygon bounds latitude = (lat - buffer, lat + buffer) longitude = (lon - buffer, lon + buffer) @@ -150,60 +139,43 @@ def run_crophealth_app(ds, lat, lon, buffer): "properties": { "style": { "stroke": True, - "color": 'red', + "color": "red", "weight": 4, "opacity": 0.8, "fill": True, "fillColor": False, "fillOpacity": 0, "showArea": True, - "clickable": True + "clickable": True, } }, "geometry": { "type": "Polygon", "coordinates": [ [ - [ - longitude[0], - latitude[0] - ], - [ - longitude[1], - latitude[0] - ], - [ - longitude[1], - latitude[1] - ], - [ - longitude[0], - latitude[1] - ], - [ - longitude[0], - latitude[0] - ] + [longitude[0], latitude[0]], + [longitude[1], latitude[0]], + [longitude[1], latitude[1]], + [longitude[0], latitude[1]], + [longitude[0], latitude[0]], ] - ] - } + ], + }, } - + # Create a map geometry from the geom_obj dictionary # center specifies where the background map view should focus on # zoom specifies how zoomed in the background map should be - loadeddata_geometry = ogr.CreateGeometryFromJson(str(geom_obj['geometry'])) + loadeddata_geometry = ogr.CreateGeometryFromJson(str(geom_obj["geometry"])) loadeddata_center = [ loadeddata_geometry.Centroid().GetY(), - loadeddata_geometry.Centroid().GetX() + loadeddata_geometry.Centroid().GetX(), ] loadeddata_zoom = 16 # define the study area map studyarea_map = Map( - center=loadeddata_center, - zoom=loadeddata_zoom, - basemap=basemaps.Esri.WorldImagery + center=loadeddata_center, zoom=loadeddata_zoom, basemap=basemaps.Esri.WorldImagery ) # define the drawing controls @@ -216,82 +188,76 @@ def run_crophealth_app(ds, lat, lon, buffer): ) # add drawing controls and data bound geometry to the map - studyarea_map.add_control(studyarea_drawctrl) - studyarea_map.add_layer(GeoJSON(data=geom_obj)) + studyarea_map.add(studyarea_drawctrl) + studyarea_map.add(GeoJSON(data=geom_obj)) # Index to count drawn polygons polygon_number = 0 # Define widgets to interact with - instruction = widgets.Output(layout={'border': '1px solid black'}) + instruction = widgets.Output(layout={"border": "1px solid black"}) with instruction: - print("Draw a polygon within the red box to view a plot of " - "average NDVI over time in that area.") + print( + "Draw a polygon within the red box to view a plot of " + "average NDVI over time in that area." + ) - info = widgets.Output(layout={'border': '1px solid black'}) + info = widgets.Output(layout={"border": "1px solid black"}) with info: print("Plot status:") - fig_display = widgets.Output(layout=widgets.Layout( - width="50%", # proportion of horizontal space taken by plot - )) + fig_display = widgets.Output( + layout=widgets.Layout( + width="50%", # proportion of horizontal space taken by plot + ) + ) with fig_display: plt.ioff() fig, ax = plt.subplots(figsize=(8, 6)) ax.set_ylim([0, 1]) - colour_list = plt.rcParams['axes.prop_cycle'].by_key()['color'] + colour_list = plt.rcParams["axes.prop_cycle"].by_key()["color"] # Function to execute each time something is drawn on the map def handle_draw(self, action, geo_json): nonlocal polygon_number # Execute behaviour based on what the user draws - if geo_json['geometry']['type'] == 'Polygon': + if geo_json["geometry"]["type"] == "Polygon": info.clear_output(wait=True) # wait=True reduces flicker effect - + # Save geojson polygon to io temporary file to be rasterized later jsonData = json.dumps(geo_json) binaryData = jsonData.encode() io = BytesIO(binaryData) io.seek(0) - + # Read the polygon as a geopandas dataframe gdf = gpd.read_file(io) gdf.crs = "EPSG:4326" # Convert the drawn geometry to pixel coordinates - xr_poly = xr_rasterize(gdf, ds.NDVI.isel(time=0), crs='EPSG:6933') + xr_poly = xr_rasterize(gdf, ds.NDVI.isel(time=0), crs="EPSG:6933") # Construct a mask to only select pixels within the drawn polygon masked_ds = ds.NDVI.where(xr_poly) - - masked_ds_mean = masked_ds.mean(dim=['x', 'y'], skipna=True) + + masked_ds_mean = masked_ds.mean(dim=["x", "y"], skipna=True) colour = colour_list[polygon_number % len(colour_list)] # Add a layer to the map to make the most recently drawn polygon # the same colour as the line on the plot - studyarea_map.add_layer( + studyarea_map.add( GeoJSON( data=geo_json, - style={ - 'color': colour, - 'opacity': 1, - 'weight': 4.5, - 'fillOpacity': 0.0 - } + style={"color": colour, "opacity": 1, "weight": 4.5, "fillOpacity": 0.0}, ) ) # add new data to the plot - xr.plot.plot( - masked_ds_mean, - marker='*', - color=colour, - ax=ax - ) + xr.plot.plot(masked_ds_mean, marker="*", color=colour, ax=ax) # reset titles back to custom ax.set_title("Average NDVI from Sentinel-2") @@ -302,7 +268,7 @@ def handle_draw(self, action, geo_json): fig_display.clear_output(wait=True) # wait=True reduces flicker effect with fig_display: display(fig) - + with info: print("Plot status: polygon sucessfully added to plot.") @@ -312,8 +278,10 @@ def handle_draw(self, action, geo_json): else: info.clear_output(wait=True) with info: - print("Plot status: this drawing tool is not currently " - "supported. Please use the polygon tool.") + print( + "Plot status: this drawing tool is not currently " + "supported. Please use the polygon tool." + ) # call to say activate handle_draw function on draw studyarea_drawctrl.on_draw(handle_draw) @@ -331,7 +299,5 @@ def handle_draw(self, action, geo_json): # +-----------+-----------+ # | info | # +-----------------------+ - ui = widgets.VBox([instruction, - widgets.HBox([studyarea_map, fig_display]), - info]) + ui = widgets.VBox([instruction, widgets.HBox([studyarea_map, fig_display]), info]) display(ui) diff --git a/Tools/deafrica_tools/app/deacoastlines.py b/Tools/deafrica_tools/app/deacoastlines.py index 2b00a486..550a24a3 100644 --- a/Tools/deafrica_tools/app/deacoastlines.py +++ b/Tools/deafrica_tools/app/deacoastlines.py @@ -1,5 +1,5 @@ """ -Digital Earth Africa Coastline widget, which can be used to +Digital Earth Africa Coastline widget, which can be used to interactively extract shoreline data using transects. """ @@ -8,51 +8,32 @@ # 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, -) +os.environ["USE_PYGEOS"] = "0" + import json -import geopandas as gpd +import warnings +from datetime import datetime from io import BytesIO + +import fiona +import geopandas as gpd import ipywidgets as widgets +import matplotlib.pyplot as plt +from ipyleaflet import GeoData, LayerGroup, WMSLayer, basemap_to_tiles, basemaps +from ipywidgets import HTML, Button, GridspecLayout, HBox, Layout, Output, VBox import deafrica_tools.app.widgetconstructors as deawidgets from deafrica_tools.coastal import get_coastlines, transect_distances -from owslib.wms import WebMapService + def make_box_layout(): return Layout( # border='solid 1px black', - margin='0px 10px 10px 0px', - padding='5px 5px 5px 5px', - width='100%', - height='100%', + margin="0px 10px 10px 0px", + padding="5px 5px 5px 5px", + width="100%", + height="100%", ) @@ -81,7 +62,7 @@ def __init__(self): ("Open Street Map", "open_street_map"), ] self.product = self.product_list[0][1] - self.mode_list = [('Distance', 'distance'), ('Width', 'width')] + self.mode_list = [("Distance", "distance"), ("Width", "width")] self.mode = self.mode_list[0][1] self.target = None self.action = None @@ -94,9 +75,14 @@ def __init__(self): # Create the Header widget header_title_text = "

Digital Earth Africa Coastlines shoreline transect extraction

" - instruction_text = "Select parameters and draw a transect on the map to extract shoreline data. In distance mode, draw a transect line starting from land that crosses multiple shorelines.
In width mode, draw a transect line that intersects shorelines at least twice. Alternatively, upload an vector file to extract shoreline data for multiple existing transects." - self.header = deawidgets.create_html( - f"{header_title_text}

{instruction_text}

") + instruction_text = ( + "Select parameters and draw a transect on the map to extract shoreline data. " + "In distance mode, draw a transect line starting from land that crosses " + "multiple shorelines.
In width mode, draw a transect line that intersects " + "shorelines at least twice. Alternatively, upload an vector file to extract " + "shoreline data for multiple existing transects." + ) + self.header = deawidgets.create_html(f"{header_title_text}

{instruction_text}

") self.header.layout = make_box_layout() ##################################### @@ -125,16 +111,22 @@ def update_geojson(target, action, geo_json): gdf_drawn_nsidc = gdf.copy().to_crs("EPSG:6933") m2_per_km2 = 10**6 area = gdf_drawn_nsidc.envelope.area.values[0] / m2_per_km2 - polyarea_label = 'Total area of DE Africa Coastlines data to extract' + polyarea_label = "Total area of DE Africa Coastlines data to extract" polyarea_text = f"{polyarea_label}: {area:.2f} km2" # Test area size if area <= 50000: - confirmation_text = ' (Area to extract falls within recommended limit; click "Extract shoreline data" to continue)' + confirmation_text = ( + ' (Area to extract falls within recommended ' + 'limit; click "Extract shoreline data" to continue)' + ) self.header.value = header_title_text + polyarea_text + confirmation_text self.gdf_drawn = gdf else: - warning_text = ' (Area to extract is too large, please select a smaller transect)' + warning_text = ( + ' (Area to extract is too large, ' + "please select a smaller transect)" + ) self.header.value = header_title_text + polyarea_text + warning_text self.gdf_drawn = None @@ -150,7 +142,7 @@ def update_geojson(target, action, geo_json): ######################################### # Create drawing tools - desired_drawtools = ['polyline'] + desired_drawtools = ["polyline"] draw_control = deawidgets.create_drawcontrol(desired_drawtools) # Load DEACoastLines WMS @@ -159,23 +151,24 @@ def update_geojson(target, action, geo_json): deacoastlines = WMSLayer( url=deacl_url, layers=deacl_layer, - format='image/png', + format="image/png", transparent=True, - attribution='DE Africa Coastlines © 2022 Digital Earth Africa') + attribution=f"DE Africa Coastlines © {datetime.now().year} Digital Earth Africa", + ) # 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' + self.map_layers.name = "Map Overlays" # Create map widget - self.m = deawidgets.create_map(map_center=(0.5273, 25.1367), - zoom_level=3, - basemap=basemaps.Esri.WorldImagery) + self.m = deawidgets.create_map( + map_center=(0.5273, 25.1367), zoom_level=3, 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) + self.m.add(draw_control) + self.m.add(self.map_layers) # Store current basemap for future use self.basemap = self.m.basemap @@ -185,18 +178,13 @@ def update_geojson(target, action, geo_json): ############################ # 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]) + 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) + fileupload_transects = widgets.FileUpload(accept="", multiple=True) #################################### # UPDATE FUNCTIONS FOR EACH WIDGET # @@ -216,24 +204,30 @@ def update_geojson(target, action, geo_json): # COLLECTION OF ALL APP CONTROLS # ################################## - parameter_selection = VBox([ - HTML("Output name:"), text_output_name, - HTML( - 'Transect extraction mode:
' - ), - mode_dropdown, - HTML("
Output files:"), - checkbox_plot, - checkbox_csv, - HTML( - "
Advanced
Upload a GeoJSON or ESRI " - "Shapefile (<5 mb) containing one or more transect lines."), - fileupload_transects - ]) - map_selection = VBox([ - HTML("
Map overlay:"), - deaoverlay_dropdown, - ]) + parameter_selection = VBox( + [ + HTML("Output name:"), + text_output_name, + HTML( + 'Transect extraction mode:
' + ), + mode_dropdown, + HTML("
Output files:"), + checkbox_plot, + checkbox_csv, + HTML( + "
Advanced
Upload a GeoJSON or ESRI " + "Shapefile (<5 mb) containing one or more transect lines." + ), + fileupload_transects, + ] + ) + map_selection = VBox( + [ + HTML("
Map overlay:"), + deaoverlay_dropdown, + ] + ) parameter_selection.layout = make_box_layout() map_selection.layout = make_box_layout() @@ -290,27 +284,31 @@ def update_fileupload_transects(self, change): # Save to file for uploaded_filename in change.new.keys(): with open(uploaded_filename, "wb") as output_file: - content = change.new[uploaded_filename]['content'] + content = change.new[uploaded_filename]["content"] output_file.write(content) with self.status_info: - try: + try: - print('Loading vector data...', end='\r') + print("Loading vector data...", end="\r") valid_files = [ - file for file in change.new.keys() - if file.lower().endswith(('.shp', '.geojson')) + file + for file in change.new.keys() + if file.lower().endswith((".shp", ".geojson")) ] valid_file = valid_files[0] - transect_gdf = (gpd.read_file(valid_file).to_crs( - "EPSG:4326").explode().reset_index(drop=True)) + transect_gdf = ( + gpd.read_file(valid_file).to_crs("EPSG:4326").explode().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'.") + 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, " @@ -318,16 +316,12 @@ def update_fileupload_transects(self, change): ) # Create a geodata - geodata = GeoData(geo_dataframe=transect_gdf, - style={ - 'color': 'black', - 'weight': 3 - }) + 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) + self.m.add(geodata) # If completed, add to attribute self.gdf_uploaded = transect_gdf @@ -336,14 +330,16 @@ def update_fileupload_transects(self, change): print( "Cannot read uploaded files. Please ensure that data is " "in either GeoJSON or ESRI Shapefile format.", - end='\r') + 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') + end="\r", + ) self.gdf_uploaded = None # Set output name @@ -375,17 +371,18 @@ def update_deaoverlay(self, change): layers=deacl_layer, format="image/png", transparent=True, - attribution="DE Africa Coastlines © 2022 Digital Earth Africa") + attribution="DE Africa Coastlines © 2024 Digital Earth Africa", + ) if self.product == "none": - self.map_layers.clear_layers() - self.map_layers.add_layer(deacoastlines) + self.map_layers.clear() + self.map_layers.add(deacoastlines) elif self.product == "open_street_map": - self.map_layers.clear_layers() + self.map_layers.clear() layer = basemap_to_tiles(basemaps.OpenStreetMap.Mapnik) - self.map_layers.add_layer(layer) - self.map_layers.add_layer(deacoastlines) + self.map_layers.add(layer) + self.map_layers.add(deacoastlines) def run_app(self, change): @@ -400,62 +397,66 @@ def run_app(self, change): # Load transects from either map or uploaded files if self.gdf_uploaded is not None: transect_gdf = self.gdf_uploaded - run_text = 'uploaded file' + 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' + 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') + print( + "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') + 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 + 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') + print(f'Analysing transect distances using "{self.mode}" mode...', end="\r") dist_df = transect_distances( transect_gdf.to_crs("EPSG:6933"), deacl_gdf.to_crs("EPSG:6933"), - mode=self.mode) + mode=self.mode, + ) # If valid data was produced: if dist_df.any(axis=None): # Successful output - print(f'DE Africa Coastlines data successfully extracted for {run_text}.') + print(f"DE Africa 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) + 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}".') @@ -463,33 +464,34 @@ def run_app(self, change): # Generate plot with self.output_plot: - fig, ax = plt.subplots(constrained_layout=True, - figsize=(15, 5.5)) + 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 Africa Coastlines transect extraction - {self.output_name}") + ax.legend(frameon=False, ncol=3, title="Transect") + ax.set_title( + f"Digital Earth Africa 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) + 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') + 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) + 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}".') @@ -502,4 +504,5 @@ def run_app(self, change): " - 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') \ No newline at end of file + end="\r", + ) diff --git a/Tools/deafrica_tools/app/forestmonitoring.py b/Tools/deafrica_tools/app/forestmonitoring.py index 4e99243c..f5e64628 100644 --- a/Tools/deafrica_tools/app/forestmonitoring.py +++ b/Tools/deafrica_tools/app/forestmonitoring.py @@ -1,19 +1,21 @@ -''' -Functions for loading and interacting with Global Forest Change data in the forest monitoring notebook, inside the Real_world_examples folder. -''' +""" +Functions for loading and interacting with +Global Forest Change data in the forest monitoring notebook, inside the Real_world_examples folder. +""" # 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' + +os.environ["USE_PYGEOS"] = "0" import json import warnings from io import BytesIO -import deafrica_tools.app.widgetconstructors as deawidgets +import fiona import geopandas as gpd import ipywidgets as widgets import matplotlib.colors as mcolors @@ -22,27 +24,27 @@ import pandas as pd import rioxarray import xarray as xr -from deafrica_tools.dask import create_local_dask_cluster -from deafrica_tools.spatial import xr_rasterize -from ipyleaflet import ( - DrawControl, - GeoData, - LayerGroup, - LayersControl, - Map, - WidgetControl, - WMSLayer, - basemap_to_tiles, - basemaps, -) -from ipywidgets import HTML, Button, GridspecLayout, HBox, Layout, Output, VBox +from ipyleaflet import GeoData, LayerGroup, basemap_to_tiles, basemaps +from ipywidgets import Button, GridspecLayout, HBox, Layout, Output, VBox from matplotlib.patches import Patch -from traitlets import Unicode +from odc.geo.geom import Geometry +from shapely.geometry import Polygon + +import deafrica_tools.app.widgetconstructors as deawidgets +from deafrica_tools.dask import create_local_dask_cluster +from deafrica_tools.spatial import add_geobox # Turn off all warnings. warnings.filterwarnings("ignore") warnings.simplefilter("ignore") + +BASE_URL = ( + "https://storage.googleapis.com/earthenginepartners-hansen/GFC-2023-v1.11/Hansen_GFC-2023-v1.11" +) +END_YEAR = 2023 + + def make_box_layout(): """ Defines a number of CSS properties that impact how a widget is laid out. @@ -54,6 +56,7 @@ def make_box_layout(): height="100%", ) + def create_expanded_button(description, button_style): """ Defines a number of CSS properties to create a button to handle mouse clicks. @@ -64,69 +67,116 @@ def create_expanded_button(description, button_style): layout=Layout(width="auto", height="auto"), ) + +def get_lon_label(lon: float) -> str: + """ + Get the longitude label given the longitude value. + Parameters + ---------- + lon: float + Longitude to get the label for. + + Returns + ------- + str: + Longitude label where negative values are East and positive values are West. + """ + if lon >= 0: + label = f"{abs(lon):03d}E" + else: + label = f"{abs(lon):03d}W" + + return label + + +def get_lat_label(lat: float) -> str: + """ + Get the latitude label given the latitude value. + Parameters + ---------- + lat: float + Latitude to get the label for. + + Returns + ------- + str: + Latitude label where negative values are South and positive values are North. + """ + if lat >= 0: + label = f"{abs(lat):02d}N" + else: + label = f"{abs(lat):02d}S" + + return label + + +def get_product_tiles() -> gpd.GeoDataFrame: + """ + Get the individual 10x10 degree granules that comprise the Global Forest Change 2000-2023 + dataset. + + Returns + ------- + gpd.GeoDataFrame + The individual 10x10 degree granules that comprise the Global Forest Change 2000-2023 + dataset + """ + granules_lon_range = np.arange(-180, 180, 10) + granules_lat_range = np.arange(80, -60, -10) + + granules = [] + labels = [] + + for lat in granules_lat_range: + for lon in granules_lon_range: + top_left = (lon, lat) + top_right = (lon + 10, lat) + bottom_right = (lon + 10, lat - 10) + bottom_left = (lon, lat - 10) + + square = Polygon([top_left, top_right, bottom_right, bottom_left, top_left]) + + granules.append(square) + + label = f"{get_lat_label(lat)}_{get_lon_label(lon)}" + labels.append(label) + # Create a GeoDataFrame from the grid polygons + granules_gdf = gpd.GeoDataFrame({"geometry": granules, "label": labels}, crs="EPSG:4326") + return granules_gdf + + def load_gfclayer(gdf_drawn, gfclayer): """ Loads the selected Global Forest Change layer for the area drawn on the map widget. """ + gdf_drawn_geopolygon = Geometry(geom=gdf_drawn.geometry.iloc[0], crs=gdf_drawn.crs) + # Configure local dask cluster. client = create_local_dask_cluster(return_client=True, display_client=True) - # Get the coordinates of the top-left corner for each Global Forest Change tile, - # covering the area of interest. - min_lat, max_lat = ( - gdf_drawn.bounds.miny.item(), - gdf_drawn.bounds.maxy.item(), - ) - min_lon, max_lon = ( - gdf_drawn.bounds.minx.item(), - gdf_drawn.bounds.maxx.item(), - ) + product_tiles = get_product_tiles() - lats = np.arange( - np.floor(min_lat / 10) * 10, np.ceil(max_lat / 10) * 10, 10 - ).astype(int) - lons = np.arange( - np.floor(min_lon / 10) * 10, np.ceil(max_lon / 10) * 10, 10 - ).astype(int) - - coord_list = [] - for lat in lats: - lat = lat + 10 - if lat >= 0: - lat_str = f"{lat:02d}N" - else: - lat_str = f"{abs(lat):02d}S" - for lon in lons: - if lon >= 0: - lon_str = f"{lon:03d}E" - else: - lon_str = f"{abs(lon):03d}W" - coord_str = f"{lat_str}_{lon_str}" - coord_list.append(coord_str) + assert product_tiles.crs == gdf_drawn_geopolygon.crs + + # Find the tiles that intersect with the area of interest + intersecting_tiles = product_tiles[product_tiles.intersects(gdf_drawn_geopolygon.geom)] + tile_labels = intersecting_tiles["label"].to_list() + + # Get the urls for the layers to load + urls = [f"{BASE_URL}_{gfclayer}_{tile_label}.tif" for tile_label in tile_labels] - # Load each Global Forest Change tile covering the area of interest. - base_url = f"https://storage.googleapis.com/earthenginepartners-hansen/GFC-2021-v1.9/Hansen_GFC-2021-v1.9_{gfclayer}_" dask_chunks = dict(x=2048, y=2048) tile_list = [] - for coord in coord_list: - tile_url = f"{base_url}{coord}.tif" - # Load the tile as an xarray.DataArray. - tile = rioxarray.open_rasterio(tile_url, chunks=dask_chunks).squeeze() + for url in urls: + tile = rioxarray.open_rasterio(url, chunks=dask_chunks).squeeze() + tile = add_geobox(tile) + tile = tile.odc.crop(gdf_drawn_geopolygon) tile_list.append(tile) # Merge the tiles into a single xarray.DataArray. ds = xr.combine_by_coords(tile_list) - # Clip the dataset using the bounds of the area of interest. - ds = ds.rio.clip_box( - minx=min_lon - 0.00025, - miny=min_lat - 0.00025, - maxx=max_lon + 0.00025, - maxy=max_lat + 0.00025, - ) - # Rename the y and x variables for DEA convention on xarray.DataArrays where crs="EPSG:4326". - ds = ds.rename({"y": "latitude", "x": "longitude"}) # Mask pixels representing no loss (encoded as 0) in the "lossyear" layer. if gfclayer == "lossyear": @@ -138,22 +188,16 @@ def load_gfclayer(gdf_drawn, gfclayer): elif gfclayer == "treecover2000": ds = ds.where(ds != 0) - # Create a mask from the area of interest GeoDataFrame. - mask = xr_rasterize(gdf_drawn, ds) - # Mask the dataset. - ds = ds.where(mask) # Convert the xarray.DataArray to a dataset. ds = ds.to_dataset(name=gfclayer) # Compute. ds = ds.compute() - # Assign the "EPSG:4326" CRS to the dataset. - ds.rio.write_crs(4326, inplace=True) - ds = ds.transpose("latitude", "longitude") # Close down the dask client. client.close() return ds + def load_all_gfclayers(gdf_drawn): gfclayers = ["treecover2000", "gain", "lossyear"] @@ -165,6 +209,7 @@ def load_all_gfclayers(gdf_drawn): dataset = xr.merge(dataset_list) return dataset + def get_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): """ Preprocess the Global Forest change "treecover2020" layer. @@ -185,16 +230,15 @@ def get_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): counts = np.unique(ds_masked, return_counts=True) # Remove the counts for pixels with the value np.nan. index = np.argwhere(np.isnan(counts[0])) - counts_dict = dict( - zip(np.delete(counts[0], index), np.delete(counts[1], index)) - ) + counts_dict = dict(zip(np.delete(counts[0], index), np.delete(counts[1], index))) # Reproject the dataset to EPSG:6933 which uses metres - ds_reprojected = ds_masked.rio.reproject("EPSG:6933") + ds_reprojected = ds.odc.reproject("EPSG:6933") # Get the area per pixel. - pixel_length = ds_reprojected.geobox.resolution[1] - m_per_km = 1000 - per_pixel_area = (pixel_length / m_per_km) ** 2 + per_pixel_area = ( + abs(ds_reprojected.odc.geobox.resolution.x * ds_reprojected.odc.geobox.resolution.y) + / 1000000 + ) # Save the results as a pandas DataFrame. df = pd.DataFrame( @@ -206,13 +250,17 @@ def get_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): ) # Get the total area. - print_statement = f'Total Forest Cover in {df["Year"].item()}: {round(df["Tree Cover in km$^2$"].item(), 4)} km2' + print_statement = ( + f'Total Forest Cover in {df["Year"].item()}: ' + f'{round(df["Tree Cover in km$^2$"].item(), 4)} km2' + ) # File name to use when exporting results. - file_name = f"forest_cover_in_2000" + file_name = "forest_cover_in_2000" return ds, df, print_statement, file_name + def get_gfclayer_gain(gfclayer_ds, gfclayer="gain"): """ Preprocess the Global Forest Change "gain" layer. @@ -229,36 +277,37 @@ def get_gfclayer_gain(gfclayer_ds, gfclayer="gain"): counts = np.unique(ds, return_counts=True) # Remove the counts for pixels with the value np.nan. index = np.argwhere(np.isnan(counts[0])) - counts_dict = dict( - zip(np.delete(counts[0], index), np.delete(counts[1], index)) - ) + counts_dict = dict(zip(np.delete(counts[0], index), np.delete(counts[1], index))) # Reproject the dataset to EPSG:6933 which uses metres. - ds_reprojected = ds.rio.reproject("EPSG:6933") + ds_reprojected = ds.odc.reproject("EPSG:6933") # Get the area per pixel. - pixel_length = ds_reprojected.geobox.resolution[1] - m_per_km = 1000 - per_pixel_area = (pixel_length / m_per_km) ** 2 + per_pixel_area = ( + abs(ds_reprojected.odc.geobox.resolution.x * ds_reprojected.odc.geobox.resolution.y) + / 1000000 + ) # Save the results as a pandas DataFrame. df = pd.DataFrame( data={ "Year": ["2000-2012"], - "Forest Cover Gain in km$^2$": np.fromiter( - counts_dict.values(), dtype=float - ) + "Forest Cover Gain in km$^2$": np.fromiter(counts_dict.values(), dtype=float) * per_pixel_area, } ) # Get the total area. - print_statement = f'Total Forest Cover Gain {df["Year"].item()}: {round(df["Forest Cover Gain in km$^2$"].item(), 4)} km2' + print_statement = ( + f'Total Forest Cover Gain {df["Year"].item()}: ' + f'{round(df["Forest Cover Gain in km$^2$"].item(), 4)} km2' + ) # File name to use when exporting results. - file_name = f"forest_cover_gain_from_2000_to_2012" + file_name = "forest_cover_gain_from_2000_to_2012" return ds, df, print_statement, file_name + def get_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear"): """ Preprocess the Global Forest Change "lossyear" layer. @@ -267,7 +316,7 @@ def get_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear" ds = gfclayer_ds[gfclayer] # Mask the dataset to the selected time range. - selected_years = list(range(start_year, end_year + 1)) + selected_years = list(range(start_year - 2000, end_year - 2000 + 1)) mask = ds.isin(selected_years) ds = ds.where(mask) @@ -281,37 +330,38 @@ def get_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear" counts = np.unique(ds, return_counts=True) # Remove the counts for pixels with the value np.nan. index = np.argwhere(np.isnan(counts[0])) - counts_dict = dict( - zip(np.delete(counts[0], index), np.delete(counts[1], index)) - ) + counts_dict = dict(zip(np.delete(counts[0], index), np.delete(counts[1], index))) # Reproject the dataset to EPSG:6933 which uses metres - ds_reprojected = ds.rio.reproject("EPSG:6933") + ds_reprojected = ds.odc.reproject("EPSG:6933") # Get the area per pixel. - pixel_length = ds_reprojected.geobox.resolution[1] - m_per_km = 1000 - per_pixel_area = (pixel_length / m_per_km) ** 2 + per_pixel_area = ( + abs(ds_reprojected.odc.geobox.resolution.x * ds_reprojected.odc.geobox.resolution.y) + / 1000000 + ) # For each year get the area of loss. # Save the results as a pandas DataFrame. df = pd.DataFrame( { "Year": 2000 + np.fromiter(counts_dict.keys(), dtype=int), - "Forest Cover Loss in km$^2$": np.fromiter( - counts_dict.values(), dtype=float - ) + "Forest Cover Loss in km$^2$": np.fromiter(counts_dict.values(), dtype=float) * per_pixel_area, } ) # Get the total area. - print_statement = f'Total Forest Cover Loss from {start_year + 2000} to {end_year + 2000}: {round(df["Forest Cover Loss in km$^2$"].sum(), 4)} km2' + print_statement = ( + f"Total Forest Cover Loss from {start_year} to {end_year}: " + f'{round(df["Forest Cover Loss in km$^2$"].sum(), 4)} km2' + ) # File name to use when exporting results. - file_name = f"forest_cover_loss_from_{start_year + 2000}_to_{end_year + 2000}" + file_name = f"forest_cover_loss_from_{start_year}_to_{end_year}" return ds, df, print_statement, file_name + def plot_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): """ Plot the Global Forest Change "treecover2000" layer. @@ -319,7 +369,8 @@ def plot_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): if get_gfclayer_treecover2000(gfclayer_ds) is None: print( - f"No Global Forest Change {gfclayer} layer data found in the selected area. Please select a new polygon over an area with data." + f"No Global Forest Change {gfclayer} layer data found in the " + "selected area. Please select a new polygon over an area with data." ) else: ds, df, print_statement, file_name = get_gfclayer_treecover2000(gfclayer_ds) @@ -331,16 +382,14 @@ def plot_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): # Define the plotting parameters. figure_width = 10 figure_length = 10 - title = f"Tree Canopy Cover for the Year 2000" + title = "Tree Canopy Cover for the Year 2000" # Plot the dataset. fig, ax = plt.subplots(figsize=(figure_width, figure_length)) im = ds.plot(cmap="Greens", add_colorbar=False, ax=ax) # Add a colorbar to the plot. cbar = plt.colorbar(mappable=im) - cbar.set_label( - "Percentage tree canopy cover for year 2000", labelpad=-65, y=0.25 - ) + cbar.set_label("Percentage tree canopy cover for year 2000", labelpad=-65, y=0.25) # Add a title to the plot. plt.title(title) # Save the plot. @@ -350,6 +399,7 @@ def plot_gfclayer_treecover2000(gfclayer_ds, gfclayer="treecover2000"): print(print_statement) + def plot_gfclayer_gain(gfclayer_ds, gfclayer="gain"): """ Plot the Global Forest Change "gain" layer. @@ -357,7 +407,8 @@ def plot_gfclayer_gain(gfclayer_ds, gfclayer="gain"): if get_gfclayer_gain(gfclayer_ds) is None: print( - f"No Global Forest Change {gfclayer} layer data found in the selected area. Please select a new polygon over an area with data." + f"No Global Forest Change {gfclayer} layer data found in the selected area. " + "Please select a new polygon over an area with data." ) else: ds, df, print_statement, file_name = get_gfclayer_gain(gfclayer_ds) @@ -370,7 +421,7 @@ def plot_gfclayer_gain(gfclayer_ds, gfclayer="gain"): color = "#6CAE75" figure_width = 10 figure_length = 10 - title = f"Forest Cover Gain from 2000 to 2012" + title = "Forest Cover Gain from 2000 to 2012" # Plot the dataset. fig, ax = plt.subplots(figsize=(figure_width, figure_length)) @@ -392,17 +443,16 @@ def plot_gfclayer_gain(gfclayer_ds, gfclayer="gain"): print(print_statement) + def plot_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear"): """ Plot the Global Forest change "lossyear" layer. """ - if ( - get_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear") - is None - ): + if get_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear") is None: print( - f"No Global Forest Change {gfclayer} layer data found in the selected area. Please select a new polygon over an area with data." + f"No Global Forest Change {gfclayer} layer data found in the selected area. " + "Please select a new polygon over an area with data." ) else: ds, df, print_statement, file_name = get_gfclayer_lossyear( @@ -418,45 +468,27 @@ def plot_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear figure_length = 15 nrows = 2 ncols = 1 - title = f"Forest Cover Loss from {start_year + 2000} to {end_year + 2000}" + title = f"Forest Cover Loss from {start_year} to {end_year}" + + num_colors = (END_YEAR - 2001) + 1 # Location of transition from one color to the next on the colormap. - color_levels = list(np.arange(1 - 0.5, 22, 1)) + color_levels = list(np.arange(1 - 0.5, num_colors + 1, 1)) + # Ticks to be displayed. - ticks = list(np.arange(1, 22)) - tick_labels = list(2000 + np.arange(1, 22)) - - # Define the color map to use when plotting. - color_list = [ - "#e6194b", - "#3cb44b", - "#ffe119", - "#4363d8", - "#f58231", - "#911eb4", - "#46f0f0", - "#f032e6", - "#bcf60c", - "#fabebe", - "#008080", - "#e6beff", - "#9a6324", - "#fffac8", - "#800000", - "#aaffc3", - "#808000", - "#ffd8b1", - "#000075", - "#808080", - "#7A306C", - ] - cmap = mcolors.ListedColormap(colors=color_list, N=21) + ticks = list(np.arange(1, num_colors + 1)) + tick_labels = list(2000 + np.arange(1, num_colors + 1)) + + # Generate evenly spaced values for the colormap + values = np.linspace(0, 1, num_colors) + # Generate evenly spaced colors in the viridis color space + colors = plt.cm.viridis(values) + # Create a custom colormap + cmap = mcolors.LinearSegmentedColormap.from_list("custom_colormap", colors, N=num_colors) norm = mcolors.BoundaryNorm(boundaries=color_levels, ncolors=cmap.N) # Plot the dataset. - fig, (ax1, ax2) = plt.subplots( - nrows, ncols, figsize=(figure_width, figure_length) - ) + fig, (ax1, ax2) = plt.subplots(nrows, ncols, figsize=(figure_width, figure_length)) im = ds.plot(ax=ax1, cmap=cmap, norm=norm, add_colorbar=False) # Add a title to the subplot. ax1.set_title(title) @@ -479,6 +511,7 @@ def plot_gfclayer_lossyear(gfclayer_ds, start_year, end_year, gfclayer="lossyear print(print_statement) + def plot_gfclayer_all(gfclayer_ds, start_year, end_year): """ Plot all the Global Forest Change Layers loaded. @@ -498,14 +531,10 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): # Define the figure. fig, ax = plt.subplots(figsize=(figure_width, figure_length)) - if ( - get_gfclayer_treecover2000( - gfclayer_ds[["treecover2000"]], gfclayer="treecover2000" - ) - is None - ): + if get_gfclayer_treecover2000(gfclayer_ds[["treecover2000"]], gfclayer="treecover2000") is None: print( - f"No Global Forest Change 'treecover2000' layer data found in the selected area. Please select a new polygon over an area with data." + "No Global Forest Change 'treecover2000' layer data found in the selected area. " + "Please select a new polygon over an area with data." ) else: ( @@ -513,18 +542,12 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): df_treecover2000, print_statement_treecover2000, file_name_treecover2000, - ) = get_gfclayer_treecover2000( - gfclayer_ds[["treecover2000"]], gfclayer="treecover2000" - ) + ) = get_gfclayer_treecover2000(gfclayer_ds[["treecover2000"]], gfclayer="treecover2000") # Plot the treecover2000 layer as the background layer. - background = ds_treecover2000.plot( - cmap=treecover_color, add_colorbar=False, ax=ax - ) + background = ds_treecover2000.plot(cmap=treecover_color, add_colorbar=False, ax=ax) # Add a colorbar to the treecover2000 plot. cbar = plt.colorbar(mappable=background) - cbar.set_label( - "Percentage tree canopy cover for year 2000", labelpad=-65, y=0.25 - ) + cbar.set_label("Percentage tree canopy cover for year 2000", labelpad=-65, y=0.25) # Export the dataframe as a csv. df_treecover2000.to_csv(f"{file_name_treecover2000}.csv", index=False) # Add the print statement to the list. @@ -534,16 +557,15 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): if get_gfclayer_gain(gfclayer_ds[["gain"]], gfclayer="gain") is None: print( - f"No Global Forest Change 'gain' layer data found in the selected area. Please select a new polygon over an area with data." + "No Global Forest Change 'gain' layer data found in the selected area. " + "Please select a new polygon over an area with data." ) else: ds_gain, df_gain, print_statement_gain, file_name_gain = get_gfclayer_gain( gfclayer_ds[["gain"]], gfclayer="gain" ) # Plot the gain layer. - ds_gain.plot( - ax=ax, cmap=mcolors.ListedColormap([gain_color]), add_colorbar=False - ) + ds_gain.plot(ax=ax, cmap=mcolors.ListedColormap([gain_color]), add_colorbar=False) # Export the dataframe as a csv. df_gain.to_csv(f"{file_name_gain}.csv", index=False) # Add the print statement to the list. @@ -552,13 +574,12 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): filename_list.append(f'"{file_name_gain}.csv"') if ( - get_gfclayer_lossyear( - gfclayer_ds[["lossyear"]], start_year, end_year, gfclayer="lossyear" - ) + get_gfclayer_lossyear(gfclayer_ds[["lossyear"]], start_year, end_year, gfclayer="lossyear") is None ): print( - f"No Global Forest Change 'lossyear' layer data found in the selected area. Please select a new polygon over an area with data." + "No Global Forest Change 'lossyear' layer data found in the selected area. " + "Please select a new polygon over an area with data." ) else: ( @@ -570,9 +591,7 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): gfclayer_ds[["lossyear"]], start_year, end_year, gfclayer="lossyear" ) # Plot the lossyear layer. - ds_lossyear.plot( - ax=ax, cmap=mcolors.ListedColormap([lossyear_color]), add_colorbar=False - ) + ds_lossyear.plot(ax=ax, cmap=mcolors.ListedColormap([lossyear_color]), add_colorbar=False) # Export the dataframe as a csv. df_lossyear.to_csv(f"{file_name_lossyear}.csv", index=False) # Add the print statement to the list. @@ -585,7 +604,7 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): [Patch(facecolor=gain_color), Patch(facecolor=lossyear_color)], [ "Global forest cover \n gain 2000–2012", - f"Global forest cover \n loss {str(2000+start_year)}-{str(2000+end_year)}", + f"Global forest cover \n loss {start_year}-{end_year}", ], loc="lower right", bbox_to_anchor=(-0.1, 0.75), @@ -597,7 +616,8 @@ def plot_gfclayer_all(gfclayer_ds, start_year, end_year): plt.show() print(*print_statement_list, sep="\n") print(*filename_list, sep="\n\t") - print(f'\nFigure saved as "{figure_fn}"'); + print(f'\nFigure saved as "{figure_fn}"') + def plot_gfclayer(gfclayer_ds, start_year, end_year, gfclayer): if gfclayer == "treecover2000": @@ -609,6 +629,7 @@ def plot_gfclayer(gfclayer_ds, start_year, end_year, gfclayer): elif gfclayer == "alllayers": plot_gfclayer_all(gfclayer_ds, start_year, end_year) + def update_map_layers(self): """ Updates map widget to add new basemap when selected @@ -622,6 +643,7 @@ def update_map_layers(self): # Add the selected basemap to the layer Group. self.map_layers.add_layer(self.basemap) + class forest_monitoring_app(HBox): def __init__(self): super().__init__() @@ -634,9 +656,7 @@ def __init__(self): header_title_text = "

Digital Earth Africa Forest Change

" instruction_text = """

Select the desired Global Forest Change layer, then zoom in and draw a polygon to select an area for which to plot the selected Global Forest Change layer. Alternatively, upload a vector file of the area of interest.

""" - self.header = deawidgets.create_html( - value=f"{header_title_text}{instruction_text}" - ) + self.header = deawidgets.create_html(value=f"{header_title_text}{instruction_text}") self.header.layout = make_box_layout() ############################ @@ -653,16 +673,12 @@ def __init__(self): # Set the default basemap to be used for the map widget / initial value for the widget. self.basemap = self.basemap_list[0][1] # Dropdown selection widget. - dropdown_basemap = deawidgets.create_dropdown( - options=self.basemap_list, value=self.basemap - ) + dropdown_basemap = deawidgets.create_dropdown(options=self.basemap_list, value=self.basemap) # Register the update function to run when a new value is selected # on the dropdown_basemap widget. dropdown_basemap.observe(self.update_basemap, "value") # Text to accompany the dropdown selection widget. - basemap_selection_html = deawidgets.create_html( - value=f"
Map overlay:" - ) + basemap_selection_html = deawidgets.create_html(value="
Map overlay:") # Combine the basemap_selection_html text and the dropdown_basemap widget in a single container. basemap_selection = VBox([basemap_selection_html, dropdown_basemap]) @@ -679,12 +695,12 @@ def __init__(self): ## Selection widget for the data time range. # Set the default time range for which to load data for. - self.start_year = 1 - self.end_year = 21 + self.start_year = 2001 + self.end_year = END_YEAR # Create the time range selector. time_range = list(range(self.start_year, self.end_year + 1)) - time_range_str = [str(2000 + i) for i in time_range] + time_range_str = [str(i) for i in time_range] timerange_options = tuple(zip(time_range_str, time_range)) timerange_selection_slide = widgets.SelectionRangeSlider( @@ -697,12 +713,10 @@ def __init__(self): timerange_selection_slide.observe(self.update_timerange, "value") # Text to accompany the timerange_selection widget. timerange_selection_html = deawidgets.create_html( - value=f"
Forest Cover Loss Time Range:" + value="
Forest Cover Loss Time Range:" ) # Combine the timerange_selection_text and the timerange_selection_slide in a single container. - timerange_selection = VBox( - [timerange_selection_html, timerange_selection_slide] - ) + timerange_selection = VBox([timerange_selection_html, timerange_selection_slide]) # Set the initial parameter for the GFC layer dataset. self.gfclayer_ds = None @@ -715,7 +729,7 @@ def __init__(self): dropdown_gfclayer.observe(self.update_gfclayer, "value") # Text to accompany the dropdown selection widget. gfclayer_selection_html = deawidgets.create_html( - value=f"
Global Forest Change Layer:" + value="
Global Forest Change Layer:" ) # Combine the gfclayer_selection_html text and the dropdown_gfclayer widget in a single container. gfclayer_selection = VBox([gfclayer_selection_html, dropdown_gfclayer]) @@ -730,24 +744,24 @@ def __init__(self): ) # Text to accompany the CheckBox widget. checkbox_max_size_html = deawidgets.create_html( - value=f"""
Override maximum size limit: - (use with caution; may cause memory issues/crashes)""" + value="""
Override maximum size limit: (use with caution; may cause memory issues/crashes)""" ) # Register the update function to run when the checkbox is ticked. # on the checkbox_max_size CheckBox checkbox_max_size.observe(self.update_checkbox_max_size, "value") # # Combine the checkbox_max_size_html text and the checkbox_max_size widget in a single container. enable_max_size = VBox([checkbox_max_size_html, checkbox_max_size]) - - # Add widget to enable uploading a geojson or ESRI shapefile. + + # Add widget to enable uploading a geojson or ESRI shapefile. self.gdf_uploaded = None fileupload_aoi = widgets.FileUpload(accept="", multiple=True) - # Register the update function to be called for the file upload. + # Register the update function to be callegit d for the file upload. fileupload_aoi.observe(self.update_fileupload_aoi, "value") - fileupload_html = deawidgets.create_html(value=f"""
Advanced
Upload a GeoJSON or ESRI Shapefile (<5 mb) containing a single area of interest.""") + fileupload_html = deawidgets.create_html( + value="""
Advanced
Upload a GeoJSON or ESRI Shapefile (<5 mb) containing a single area of interest.""" + ) fileupload = VBox([fileupload_html, fileupload_aoi]) - - + ## Put the app controls widgets into a single container. parameter_selection = VBox( [ @@ -755,20 +769,16 @@ def __init__(self): gfclayer_selection, timerange_selection, enable_max_size, - fileupload + fileupload, ] ) parameter_selection.layout = make_box_layout() ## Button to click to run the app. - run_button = create_expanded_button( - description="Generate plot", button_style="info" - ) + run_button = create_expanded_button(description="Generate plot", button_style="info") # Register the update function to be called when the run_button button # is clicked. run_button.on_click(self.run_app) - - ########################### # WIDGETS FOR APP OUTPUTS # @@ -793,13 +803,13 @@ def __init__(self): # Name of the Layer Group layer. self.map_layers.name = "Map Overlays" # Add the empty Layer Group as a single layer to the map widget. - self.m.add_layer(self.map_layers) + self.m.add(self.map_layers) # Create the desired drawing tools. desired_drawtools = ["rectangle", "polygon"] draw_control = deawidgets.create_drawcontrol(desired_drawtools) # Add drawing tools to the map widget. - self.m.add_control(draw_control) + self.m.add(draw_control) # Set the initial parameters for the drawing tools. self.target = None self.action = None @@ -810,7 +820,6 @@ def __init__(self): ##################################### def handle_draw(target, action, geo_json): - """ Defines the action to take once something is drawn on the map widget. @@ -818,7 +827,7 @@ def handle_draw(target, action, geo_json): # Remove previously uploaded data if present self.gdf_uploaded = None fileupload_aoi._counter = 0 - + self.target = target self.action = action @@ -838,20 +847,14 @@ def handle_draw(target, action, geo_json): m2_per_ha = 10000 area = gdf_drawn_nsidc.area.values[0] / m2_per_ha - polyarea_label = ( - f"Total area of Global Forest Change {self.gfclayer} layer to load" - ) + polyarea_label = f"Total area of Global Forest Change {self.gfclayer} layer to load" polyarea_text = f"{polyarea_label}: {area:.2f} ha" # Test the size of the polygon drawn. if self.max_size: - confirmation_text = """ - (Overriding maximum size limit; use with caution as may lead to memory issues)""" + 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 + header_title_text + instruction_text + polyarea_text + confirmation_text ) self.gdf_drawn = gdf elif area <= 50000: @@ -859,10 +862,7 @@ def handle_draw(target, action, geo_json): (Area to extract falls within recommended 50000 ha limit)""" self.header.value = ( - header_title_text - + instruction_text - + polyarea_text - + confirmation_text + header_title_text + instruction_text + polyarea_text + confirmation_text ) self.gdf_drawn = gdf else: @@ -886,9 +886,7 @@ def handle_draw(target, action, geo_json): grid_columns = 11 grid_height = "1500px" grid_width = "auto" - grid = GridspecLayout( - grid_rows, grid_columns, height=grid_height, width=grid_width - ) + grid = GridspecLayout(grid_rows, grid_columns, height=grid_height, width=grid_width) # Place app widgets and components in app layout. # [rows, columns] @@ -911,7 +909,6 @@ def update_basemap(self, change): selected value of the dropdown_basemap widget. """ self.basemap = change.new - self.output_plot_basemap = get_basemap(self.basemap.url) update_map_layers(self) def update_gfclayer(self, change): @@ -932,42 +929,40 @@ def update_checkbox_max_size(self, change): checkbox_max_size CheckBox is checked. """ self.max_size = change.new - + def update_fileupload_aoi(self, change): # Clear any drawn data if present self.gdf_drawn = None - + # Save to file for uploaded_filename in change.new.keys(): with open(uploaded_filename, "wb") as output_file: - content = change.new[uploaded_filename]['content'] + content = change.new[uploaded_filename]["content"] output_file.write(content) with self.status_info: - try: + try: - print('Loading vector data...', end='\r') + print("Loading vector data...", end="\r") valid_files = [ - file for file in change.new.keys() - if file.lower().endswith(('.shp', '.geojson')) + file + for file in change.new.keys() + if file.lower().endswith((".shp", ".geojson")) ] valid_file = valid_files[0] - aoi_gdf = (gpd.read_file(valid_file).to_crs( - "EPSG:4326").explode().reset_index(drop=True)) + aoi_gdf = ( + gpd.read_file(valid_file).to_crs("EPSG:4326").explode().reset_index(drop=True) + ) # Create a geodata - geodata = GeoData(geo_dataframe=aoi_gdf, - style={ - 'color': 'black', - 'weight': 3 - }) + geodata = GeoData(geo_dataframe=aoi_gdf, style={"color": "black", "weight": 3}) # Add to map xmin, ymin, xmax, ymax = aoi_gdf.total_bounds self.m.fit_bounds([[ymin, xmin], [ymax, xmax]]) - self.m.add_layer(geodata) + self.m.add(geodata) # If completed, add to attribute self.gdf_uploaded = aoi_gdf @@ -976,14 +971,16 @@ def update_fileupload_aoi(self, change): print( "Cannot read uploaded files. Please ensure that data is " "in either GeoJSON or ESRI Shapefile format.", - end='\r') + 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') + end="\r", + ) self.gdf_uploaded = None def run_app(self, change): @@ -991,7 +988,7 @@ def run_app(self, change): # Clear progress bar and output areas before running. self.status_info.clear_output() self.output_plot.clear_output() - + with self.status_info: # Load the area of interest from the map or uploaded files. if self.gdf_uploaded is not None: @@ -999,13 +996,15 @@ def run_app(self, change): elif self.gdf_drawn is not None: aoi_gdf = self.gdf_drawn else: - print(f'No valid polygon drawn on the map or uploaded. Please draw a valid a transect on the map, or upload a GeoJSON or ESRI Shapefile.', - end='\r') + print( + "No valid polygon drawn on the map or uploaded. Please draw a valid a transect on the map, or upload a GeoJSON or ESRI Shapefile.", + end="\r", + ) aoi_gdf = None # If valid area of interest data returned. Load the selected Global Forest Change data. if aoi_gdf is not None: - + if self.gfclayer_ds is None: if self.gfclayer != "alllayers": self.gfclayer_ds = load_gfclayer(gdf_drawn=aoi_gdf, gfclayer=self.gfclayer) @@ -1013,14 +1012,18 @@ def run_app(self, change): self.gfclayer_ds = load_all_gfclayers(gdf_drawn=aoi_gdf) else: print("Using previously loaded data") - + # Plot the selected Global Forest Change layer. if self.gfclayer_ds is not None: with self.output_plot: - plot_gfclayer(gfclayer_ds=self.gfclayer_ds, - start_year=self.start_year, - end_year=self.end_year, - gfclayer=self.gfclayer) + plot_gfclayer( + gfclayer_ds=self.gfclayer_ds, + start_year=self.start_year, + end_year=self.end_year, + gfclayer=self.gfclayer, + ) else: with self.status_info: - print(f"No Global Forest Change {self.gfclayer} layer data found in the selected area. Please select a new polygon over an area with data.") \ No newline at end of file + print( + f"No Global Forest Change {self.gfclayer} layer data found in the selected area. Please select a new polygon over an area with data." + ) diff --git a/Tools/deafrica_tools/app/geomedian.py b/Tools/deafrica_tools/app/geomedian.py index 4c7f45e1..29653284 100644 --- a/Tools/deafrica_tools/app/geomedian.py +++ b/Tools/deafrica_tools/app/geomedian.py @@ -1,126 +1,421 @@ """ Geomedian widget: generates an interactive visualisation of -the geomedian summary statistic. +the geomedian summary statistic. """ -# Load modules +# Load modules import ipywidgets as widgets import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D import numpy as np import xarray as xr from odc.algo import xr_geomedian + def run_app(): - """ - An interactive app that allows users to visualise the difference between the median and geomedian time-series summary statistics. By modifying the red-green-blue values of three timesteps for a given pixel, the user changes the output summary statistics. - - This allows a visual representation of the difference through the output values, RGB colour, as well as showing values plotted as a vector on a 3-dimensional space. - + An interactive app that allows users to visualise the difference between the median and geomedian + time-series summary statistics. By modifying the red-green-blue values of three timesteps for a given pixel, + the user changes the output summary statistics. + + This allows a visual representation of the difference through the output values, + RGB colour, as well as showing values plotted as a vector on a 3-dimensional space. + Last modified: December 2021 """ - + # Define the red-green-blue sliders for timestep 1 - p1r = widgets.IntSlider(description='Red', max=255, value=58) - p1g = widgets.IntSlider(description='Green', max=255, value=153) - p1b = widgets.IntSlider(description='Blue', max=255, value=68) + p1r = widgets.IntSlider(description="Red", max=255, value=58) + p1g = widgets.IntSlider(description="Green", max=255, value=153) + p1b = widgets.IntSlider(description="Blue", max=255, value=68) # Define the red-green-blue sliders for timestep 2 - p2r = widgets.IntSlider(description='Red', max=255, value=208) - p2g = widgets.IntSlider(description='Green', max=255, value=221) - p2b = widgets.IntSlider(description='Blue', max=255, value=203) + p2r = widgets.IntSlider(description="Red", max=255, value=208) + p2g = widgets.IntSlider(description="Green", max=255, value=221) + p2b = widgets.IntSlider(description="Blue", max=255, value=203) # Define the red-green-blue sliders for timestep 3 - p3r = widgets.IntSlider(description='Red', max=255, value=202) - p3g = widgets.IntSlider(description='Green', max=255, value=82) - p3b = widgets.IntSlider(description='Blue', max=255, value=33) + p3r = widgets.IntSlider(description="Red", max=255, value=202) + p3g = widgets.IntSlider(description="Green", max=255, value=82) + p3b = widgets.IntSlider(description="Blue", max=255, value=33) # Define the median calculation for the timesteps def f(p1r, p1g, p1b, p2r, p2g, p2b, p3r, p3g, p3b): - print('Red Median = {}'.format(np.median([p1r, p2r, p3r]))) - print('Green Median = {}'.format(np.median([p1g, p2g, p3g]))) - print('Blue Median = {}'.format(np.median([p1b, p2b, p3b]))) + print("Red Median = {}".format(np.median([p1r, p2r, p3r]))) + print("Green Median = {}".format(np.median([p1g, p2g, p3g]))) + print("Blue Median = {}".format(np.median([p1b, p2b, p3b]))) # Define the geomedian calculation for the timesteps def g(p1r, p1g, p1b, p2r, p2g, p2b, p3r, p3g, p3b): - print('Red Geomedian = {:.2f}'.format(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).red.values.ravel()[0])) - print('Green Geomedian = {:.2f}'.format(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).green.values.ravel()[0])) - print('Blue Geomedian = {:.2f}'.format(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).blue.values.ravel()[0])) + print( + "Red Geomedian = {:.2f}".format( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).red.values.ravel()[0] + ) + ) + print( + "Green Geomedian = {:.2f}".format( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).green.values.ravel()[0] + ) + ) + print( + "Blue Geomedian = {:.2f}".format( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).blue.values.ravel()[0] + ) + ) # Define the Timestep 1 box colour def h(p1r, p1g, p1b): - fig1, axes1 = plt.subplots(figsize=(2,2)) + fig1, axes1 = plt.subplots(figsize=(2, 2)) fig1 = plt.imshow([[(p1r, p1g, p1b)]]) - axes1.set_title('Timestep 1') - axes1.axis('off') + axes1.set_title("Timestep 1") + axes1.axis("off") plt.show(fig1) # Define the Timestep 2 box colour - def hh(p2r, p2g, p2b): - fig2, axes2 = plt.subplots(figsize=(2,2)) + def hh(p2r, p2g, p2b): + fig2, axes2 = plt.subplots(figsize=(2, 2)) fig2 = plt.imshow([[(p2r, p2g, p2b)]]) - axes2.set_title('Timestep 2') - axes2.axis('off') + axes2.set_title("Timestep 2") + axes2.axis("off") plt.show(fig2) # Define the Timestep 3 box colour - def hhh(p3r, p3g, p3b): - fig3, axes3 = plt.subplots(figsize=(2,2)) + def hhh(p3r, p3g, p3b): + fig3, axes3 = plt.subplots(figsize=(2, 2)) fig3 = plt.imshow([[(p3r, p3g, p3b)]]) - axes3.set_title('Timestep 3') - axes3.axis('off') + axes3.set_title("Timestep 3") + axes3.axis("off") plt.show(fig3) # Define the Median RGB colour box def i(p1r, p1g, p1b, p2r, p2g, p2b, p3r, p3g, p3b): - fig4, axes4 = plt.subplots(figsize=(3,3)) - fig4 = plt.imshow([[(int(np.median([p1r, p2r, p3r])), int(np.median([p1g, p2g, p3g])), int(np.median([p1b, p2b, p3b])))]]) - axes4.set_title('Median RGB - All timesteps') - axes4.axis('off') + fig4, axes4 = plt.subplots(figsize=(3, 3)) + fig4 = plt.imshow( + [ + [ + ( + int(np.median([p1r, p2r, p3r])), + int(np.median([p1g, p2g, p3g])), + int(np.median([p1b, p2b, p3b])), + ) + ] + ] + ) + axes4.set_title("Median RGB - All timesteps") + axes4.axis("off") plt.show(fig4) # Define the Geomedian RGB colour box def ii(p1r, p1g, p1b, p2r, p2g, p2b, p3r, p3g, p3b): - fig5, axes5 = plt.subplots(figsize=(3,3)) - fig5 = plt.imshow([[(int(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).red.values.ravel()[0]), int(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).green.values.ravel()[0]), int(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).blue.values.ravel()[0]))]]) - axes5.set_title('Geomedian RGB - All timesteps') - axes5.axis('off') + fig5, axes5 = plt.subplots(figsize=(3, 3)) + fig5 = plt.imshow( + [ + [ + ( + int( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).red.values.ravel()[0] + ), + int( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).green.values.ravel()[0] + ), + int( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).blue.values.ravel()[0] + ), + ) + ] + ] + ) + axes5.set_title("Geomedian RGB - All timesteps") + axes5.axis("off") plt.show(fig5) - # Define 3-D axis to display vectors on + # Define 3-D axis to display vectors on def j(p1r, p1g, p1b, p2r, p2g, p2b, p3r, p3g, p3b): fig6 = plt.figure() - axes6 = fig6.add_subplot(111, projection='3d') - x = [p1r, p2r, p3r, int(np.median([p1r, p2r, p3r])), int(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).red.values.ravel()[0])] - y = [p1g, p2g, p3g, int(np.median([p1g, p2g, p3g])), int(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).green.values.ravel()[0])] - z = [p1b, p2b, p3b, int(np.median([p1b, p2b, p3b])), int(xr_geomedian(xr.Dataset({"red": (("x", "y", "time"), [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]]), "green": (("x", "y", "time"), [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]]), "blue": (("x", "y", "time"), [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]])})).blue.values.ravel()[0])] - labels = [' 1', ' 2', ' 3', ' median', ' geomedian'] - axes6.scatter(x, y, z, c=['black','black','black','r', 'blue'], marker='o') - axes6.set_xlabel('Red') - axes6.set_ylabel('Green') - axes6.set_zlabel('Blue') + axes6 = fig6.add_subplot(111, projection="3d") + x = [ + p1r, + p2r, + p3r, + int(np.median([p1r, p2r, p3r])), + int( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).red.values.ravel()[0] + ), + ] + y = [ + p1g, + p2g, + p3g, + int(np.median([p1g, p2g, p3g])), + int( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).green.values.ravel()[0] + ), + ] + z = [ + p1b, + p2b, + p3b, + int(np.median([p1b, p2b, p3b])), + int( + xr_geomedian( + xr.Dataset( + { + "red": ( + ("x", "y", "time"), + [[[np.float32(p1r), np.float32(p2r), np.float32(p3r)]]], + ), + "green": ( + ("x", "y", "time"), + [[[np.float32(p1g), np.float32(p2g), np.float32(p3g)]]], + ), + "blue": ( + ("x", "y", "time"), + [[[np.float32(p1b), np.float32(p2b), np.float32(p3b)]]], + ), + } + ) + ).blue.values.ravel()[0] + ), + ] + labels = [" 1", " 2", " 3", " median", " geomedian"] + axes6.scatter(x, y, z, c=["black", "black", "black", "r", "blue"], marker="o") + axes6.set_xlabel("Red") + axes6.set_ylabel("Green") + axes6.set_zlabel("Blue") axes6.set_xlim3d(0, 255) axes6.set_ylim3d(0, 255) axes6.set_zlim3d(0, 255) for ax, ay, az, label in zip(x, y, z, labels): axes6.text(ax, ay, az, label) - plt.title('Each band represents a dimension.') + plt.title("Each band represents a dimension.") plt.show() # Define outputs - outf = widgets.interactive_output(f, {'p1r': p1r, 'p2r': p2r,'p3r': p3r, 'p1g': p1g, 'p2g': p2g,'p3g': p3g, 'p1b': p1b, 'p2b': p2b,'p3b': p3b}) - outg = widgets.interactive_output(g, {'p1r': p1r, 'p2r': p2r,'p3r': p3r, 'p1g': p1g, 'p2g': p2g,'p3g': p3g, 'p1b': p1b, 'p2b': p2b,'p3b': p3b}) + outf = widgets.interactive_output( + f, + { + "p1r": p1r, + "p2r": p2r, + "p3r": p3r, + "p1g": p1g, + "p2g": p2g, + "p3g": p3g, + "p1b": p1b, + "p2b": p2b, + "p3b": p3b, + }, + ) + outg = widgets.interactive_output( + g, + { + "p1r": p1r, + "p2r": p2r, + "p3r": p3r, + "p1g": p1g, + "p2g": p2g, + "p3g": p3g, + "p1b": p1b, + "p2b": p2b, + "p3b": p3b, + }, + ) + + outh = widgets.interactive_output(h, {"p1r": p1r, "p1g": p1g, "p1b": p1b}) + outhh = widgets.interactive_output(hh, {"p2r": p2r, "p2g": p2g, "p2b": p2b}) + outhhh = widgets.interactive_output(hhh, {"p3r": p3r, "p3g": p3g, "p3b": p3b}) - outh = widgets.interactive_output(h, {'p1r': p1r, 'p1g': p1g, 'p1b': p1b}) - outhh = widgets.interactive_output(hh, {'p2r': p2r, 'p2g': p2g, 'p2b': p2b}) - outhhh = widgets.interactive_output(hhh, {'p3r': p3r, 'p3g': p3g, 'p3b': p3b}) + outi = widgets.interactive_output( + i, + { + "p1r": p1r, + "p2r": p2r, + "p3r": p3r, + "p1g": p1g, + "p2g": p2g, + "p3g": p3g, + "p1b": p1b, + "p2b": p2b, + "p3b": p3b, + }, + ) + outii = widgets.interactive_output( + ii, + { + "p1r": p1r, + "p2r": p2r, + "p3r": p3r, + "p1g": p1g, + "p2g": p2g, + "p3g": p3g, + "p1b": p1b, + "p2b": p2b, + "p3b": p3b, + }, + ) - outi = widgets.interactive_output(i, {'p1r': p1r, 'p2r': p2r,'p3r': p3r, 'p1g': p1g, 'p2g': p2g,'p3g': p3g, 'p1b': p1b, 'p2b': p2b,'p3b': p3b}) - outii = widgets.interactive_output(ii, {'p1r': p1r, 'p2r': p2r,'p3r': p3r, 'p1g': p1g, 'p2g': p2g,'p3g': p3g, 'p1b': p1b, 'p2b': p2b,'p3b': p3b}) + outj = widgets.interactive_output( + j, + { + "p1r": p1r, + "p2r": p2r, + "p3r": p3r, + "p1g": p1g, + "p2g": p2g, + "p3g": p3g, + "p1b": p1b, + "p2b": p2b, + "p3b": p3b, + }, + ) - outj = widgets.interactive_output(j, {'p1r': p1r, 'p2r': p2r,'p3r': p3r, 'p1g': p1g, 'p2g': p2g,'p3g': p3g, 'p1b': p1b, 'p2b': p2b,'p3b': p3b}) + app_output = widgets.HBox( + [ + widgets.VBox( + [ + widgets.HBox([outh, widgets.VBox([p1r, p1g, p1b])]), + widgets.HBox([outhh, widgets.VBox([p2r, p2g, p2b])]), + widgets.HBox([outhhh, widgets.VBox([p3r, p3g, p3b])]), + ] + ), + widgets.VBox( + [widgets.HBox([widgets.VBox([outf, outi]), widgets.VBox([outg, outii])]), outj] + ), + ] + ) - app_output = widgets.HBox([widgets.VBox([widgets.HBox([outh, widgets.VBox([ p1r, p1g, p1b])]), widgets.HBox([outhh, widgets.VBox([p2r, p2g, p2b])]), widgets.HBox([outhhh, widgets.VBox([ p3r, p3g, p3b])])]), widgets.VBox([widgets.HBox([widgets.VBox([outf, outi]), widgets.VBox([outg, outii])]), outj])]) - - return app_output \ No newline at end of file + return app_output diff --git a/Tools/deafrica_tools/app/imageexport.py b/Tools/deafrica_tools/app/imageexport.py index 8b0cf1c3..376c70ba 100644 --- a/Tools/deafrica_tools/app/imageexport.py +++ b/Tools/deafrica_tools/app/imageexport.py @@ -2,25 +2,24 @@ Create an interactive map for selecting satellite imagery and exporting image files. """ +import itertools + # Load modules import datacube -import itertools -import numpy as np import matplotlib.pyplot as plt -from odc.ui import select_on_a_map -from datacube.utils.geometry import CRS +import numpy as np from datacube.utils import masking +from datacube.utils.geometry import CRS +from ipyleaflet import WMSLayer, basemap_to_tiles, basemaps +from odc.ui import select_on_a_map from skimage import exposure -from ipyleaflet import (WMSLayer, basemaps, basemap_to_tiles) from traitlets import Unicode -from deafrica_tools.spatial import reverse_geocode from deafrica_tools.dask import create_local_dask_cluster +from deafrica_tools.spatial import reverse_geocode -def select_region_app(date, - satellites, - size_limit=10000): +def select_region_app(date, satellites, size_limit=10000): """ An interactive app that allows the user to select a region from a map using imagery from Sentinel-2 and Landsat. The output of this @@ -35,7 +34,7 @@ def select_region_app(date, The exact date used to plot imagery on the interactive map (e.g. ``date='1988-01-01'``). satellites : str - The satellite data to plot on the interactive map. The + The satellite data to plot on the interactive map. The following options are supported: ``'Landsat-9'``: data from the Landsat 9 satellite @@ -66,61 +65,67 @@ def select_region_app(date, # Load DEA WMS class TimeWMSLayer(WMSLayer): - time = Unicode('').tag(sync=True, o=True) + time = Unicode("").tag(sync=True, o=True) # WMS layers wms_params = { - 'Landsat-9': 'ls9_sr', - 'Landsat-8': 'ls8_sr', - 'Landsat-7': 'ls7_sr', - 'Landsat-5': 'ls5_sr', - 'Sentinel-2': 's2_l2a', - 'Sentinel-2 geomedian': 'gm_s2_annual' + "Landsat-9": "ls9_sr", + "Landsat-8": "ls8_sr", + "Landsat-7": "ls7_sr", + "Landsat-5": "ls5_sr", + "Sentinel-2": "s2_l2a", + "Sentinel-2 geomedian": "gm_s2_annual", } - time_wms = TimeWMSLayer(url='https://ows.digitalearth.africa/', - layers=wms_params[satellites], - time=date, - format='image/png', - transparent=True, - attribution='Digital Earth Africa') + time_wms = TimeWMSLayer( + url="https://ows.digitalearth.africa/", + layers=wms_params[satellites], + time=date, + format="image/png", + transparent=True, + attribution="Digital Earth Africa", + ) # Plot interactive map to select area basemap = basemap_to_tiles(basemaps.OpenStreetMap.Mapnik) - geopolygon = select_on_a_map(height='1000px', - layers=( - basemap, - time_wms, - ), - center=(4, 20), - zoom=4) + geopolygon = select_on_a_map( + height="1000px", + layers=( + basemap, + time_wms, + ), + center=(4, 20), + zoom=4, + ) # Test size of selected area - area = geopolygon.to_crs(crs=CRS('epsg:6933')).area / 1000000 + area = geopolygon.to_crs(crs=CRS("epsg:6933")).area / 1000000 if area > size_limit: - print(f'Warning: Your selected area is {area:.00f} sq km. ' - f'Please select an area of less than {size_limit} sq km.' - f'\nTo select a smaller area, re-run the cell ' - f'above and draw a new polygon.') + print( + f"Warning: Your selected area is {area:.00f} sq km. " + f"Please select an area of less than {size_limit} sq km." + f"\nTo select a smaller area, re-run the cell " + f"above and draw a new polygon." + ) else: - return {'geopolygon': geopolygon, - 'date': date, - 'satellites': satellites} - - -def export_image_app(geopolygon, - date, - satellites, - style='True colour', - resolution=None, - vmin=0, - vmax=2000, - percentile_stretch=None, - power=None, - image_proc_funcs=None, - output_format="jpg", - standardise_name=False): + return {"geopolygon": geopolygon, "date": date, "satellites": satellites} + + +def export_image_app( + geopolygon, + date, + satellites, + style="True colour", + resolution=None, + vmin=0, + vmax=2000, + percentile_stretch=None, + power=None, + image_proc_funcs=None, + output_format="jpg", + standardise_name=False, +): """ Exports Digital Earth Africa satellite data as an image file based on the extent and time period selected using @@ -146,7 +151,7 @@ def export_image_app(geopolygon, The exact date used to extract imagery (e.g. `date='1988-01-01'`). satellites : str - The satellite data to be used to extract imagery. The + The satellite data to be used to extract imagery. The following options are supported: ``'Landsat-9'``: data from the Landsat 9 satellite @@ -206,53 +211,53 @@ def export_image_app(geopolygon, ########################### sat_params = { - 'Landsat-9': { - 'products': ['ls9_sr'], - 'resolution': [-30, 30], - 'styles': { - 'True colour': ['red', 'green', 'blue'], - 'False colour': ['swir_1', 'nir', 'green'] - } + "Landsat-9": { + "products": ["ls9_sr"], + "resolution": [-30, 30], + "styles": { + "True colour": ["red", "green", "blue"], + "False colour": ["swir_1", "nir", "green"], + }, }, - 'Landsat-8': { - 'products': ['ls8_sr'], - 'resolution': [-30, 30], - 'styles': { - 'True colour': ['red', 'green', 'blue'], - 'False colour': ['swir_1', 'nir', 'green'] - } + "Landsat-8": { + "products": ["ls8_sr"], + "resolution": [-30, 30], + "styles": { + "True colour": ["red", "green", "blue"], + "False colour": ["swir_1", "nir", "green"], + }, }, - 'Landsat-7': { - 'products': ['ls7_sr'], - 'resolution': [-30, 30], - 'styles': { - 'True colour': ['red', 'green', 'blue'], - 'False colour': ['swir_1', 'nir', 'green'] - } + "Landsat-7": { + "products": ["ls7_sr"], + "resolution": [-30, 30], + "styles": { + "True colour": ["red", "green", "blue"], + "False colour": ["swir_1", "nir", "green"], + }, }, - 'Landsat-5': { - 'products': ['ls5_sr'], - 'resolution': [-30, 30], - 'styles': { - 'True colour': ['red', 'green', 'blue'], - 'False colour': ['swir_1', 'nir', 'green'] - } + "Landsat-5": { + "products": ["ls5_sr"], + "resolution": [-30, 30], + "styles": { + "True colour": ["red", "green", "blue"], + "False colour": ["swir_1", "nir", "green"], + }, }, - 'Sentinel-2': { - 'products': ['s2_l2a'], - 'resolution': [-10, 10], - 'styles': { - 'True colour': ['red', 'green', 'blue'], - 'False colour': ['swir_2', 'nir_1', 'green'] - } + "Sentinel-2": { + "products": ["s2_l2a"], + "resolution": [-10, 10], + "styles": { + "True colour": ["red", "green", "blue"], + "False colour": ["swir_2", "nir_1", "green"], + }, }, - 'Sentinel-2 geomedian': { - 'products': ['gm_s2_annual'], - 'resolution': [-10, 10], - 'styles': { - 'True colour': ['red', 'green', 'blue'], - 'False colour': ['swir_2', 'nir_1', 'green'] - } + "Sentinel-2 geomedian": { + "products": ["gm_s2_annual"], + "resolution": [-10, 10], + "styles": { + "True colour": ["red", "green", "blue"], + "False colour": ["swir_2", "nir_1", "green"], + }, }, } @@ -261,7 +266,7 @@ def export_image_app(geopolygon, ############# # Connect to datacube database - dc = datacube.Datacube(app='Exporting_satellite_images') + dc = datacube.Datacube(app="Exporting_satellite_images") # Configure local dask cluster client = create_local_dask_cluster(return_client=True) @@ -269,48 +274,40 @@ def export_image_app(geopolygon, # Create query after adjusting interval time to UTC by # adding a UTC offset of -10 hours. start_date = np.datetime64(date) - query_params = { - 'time': (str(start_date)), - 'geopolygon': geopolygon - } + query_params = {"time": (str(start_date)), "geopolygon": geopolygon} # Find matching datasets - dss = [ - dc.find_datasets(product=i, **query_params) - for i in sat_params[satellites]['products'] - ] + dss = [dc.find_datasets(product=i, **query_params) for i in sat_params[satellites]["products"]] dss = list(itertools.chain.from_iterable(dss)) # Get CRS and sensor crs = str(dss[0].crs) - if satellites == 'Sentinel-2 geomedian': + if satellites == "Sentinel-2 geomedian": sensor = satellites else: - sensor = dss[0].metadata_doc['properties']['eo:platform'].capitalize() - sensor = sensor[0:-1].replace('_', '-') + sensor[-1].capitalize() + sensor = dss[0].metadata_doc["properties"]["eo:platform"].capitalize() + sensor = sensor[0:-1].replace("_", "-") + sensor[-1].capitalize() # Use resolution if provided, otherwise use default if resolution: - sat_params[satellites]['resolution'] = resolution + sat_params[satellites]["resolution"] = resolution load_params = { - 'output_crs': crs, - 'resolution': sat_params[satellites]['resolution'], - 'resampling': 'bilinear' + "output_crs": crs, + "resolution": sat_params[satellites]["resolution"], + "resampling": "bilinear", } # Load data from datasets - ds = dc.load(datasets=dss, - measurements=sat_params[satellites]['styles'][style], - group_by='solar_day', - dask_chunks={ - 'time': 1, - 'x': 3000, - 'y': 3000 - }, - **load_params, - **query_params) + ds = dc.load( + datasets=dss, + measurements=sat_params[satellites]["styles"][style], + group_by="solar_day", + dask_chunks={"time": 1, "x": 3000, "y": 3000}, + **load_params, + **query_params, + ) ds = masking.mask_invalid_data(ds) rgb_array = ds.isel(time=0).to_array().values @@ -322,18 +319,16 @@ def export_image_app(geopolygon, # Create unique file name centre_coords = geopolygon.centroid.coords[0][::-1] site = reverse_geocode(coords=centre_coords) - fname = (f"{sensor} - {date} - {site} - {style}, " - f"{load_params['resolution'][1]} m resolution.{output_format}") + fname = ( + f"{sensor} - {date} - {site} - {style}, " + f"{load_params['resolution'][1]} m resolution.{output_format}" + ) # Remove spaces and commas if requested if standardise_name: - fname = fname.replace(' - ', '_').replace(', ', - '-').replace(' ', - '-').lower() + fname = fname.replace(" - ", "_").replace(", ", "-").replace(" ", "-").lower() - print( - f'\nExporting image to {fname}.\nThis may take several minutes to complete...' - ) + print(f"\nExporting image to {fname}.\nThis may take several minutes to complete...") # Convert to numpy array rgb_array = np.transpose(rgb_array, axes=[1, 2, 0]) @@ -350,14 +345,14 @@ def export_image_app(geopolygon, vmin, vmax = vmin**power, vmax**power # Rescale/stretch imagery between vmin and vmax - rgb_rescaled = exposure.rescale_intensity(rgb_array.astype(float), - in_range=(vmin, vmax), - out_range=(0.0, 1.0)) + rgb_rescaled = exposure.rescale_intensity( + rgb_array.astype(float), in_range=(vmin, vmax), out_range=(0.0, 1.0) + ) # Apply image processing funcs if image_proc_funcs: for i, func in enumerate(image_proc_funcs): - print(f'Applying custom function {i + 1}') + print(f"Applying custom function {i + 1}") rgb_rescaled = func(rgb_rescaled) # Plot RGB @@ -369,4 +364,4 @@ def export_image_app(geopolygon, # Close dask client client.shutdown() - print('Finished exporting image.') + print("Finished exporting image.") diff --git a/Tools/deafrica_tools/app/wetlandsinsighttool.py b/Tools/deafrica_tools/app/wetlandsinsighttool.py index e8b6ce52..a3be100c 100644 --- a/Tools/deafrica_tools/app/wetlandsinsighttool.py +++ b/Tools/deafrica_tools/app/wetlandsinsighttool.py @@ -8,52 +8,33 @@ # 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 datacube -import warnings -import seaborn as sns -import matplotlib.pyplot as plt -from datacube.utils.geometry import CRS -from ipyleaflet import ( - WMSLayer, - basemaps, - basemap_to_tiles, - Map, - DrawControl, - WidgetControl, - LayerGroup, - LayersControl, -) -from traitlets import Unicode -from ipywidgets import ( - GridspecLayout, - Button, - Layout, - HBox, - VBox, - HTML, - Output, -) +os.environ["USE_PYGEOS"] = "0" + import json -import geopandas as gpd +import warnings from io import BytesIO -from dask.diagnostics import ProgressBar + +import geopandas as gpd +import matplotlib.pyplot as plt +import seaborn as sns +from ipyleaflet import LayerGroup, basemap_to_tiles, basemaps +from ipywidgets import HTML, Button, GridspecLayout, HBox, Layout, Output, VBox import deafrica_tools +import deafrica_tools.app.widgetconstructors as deawidgets from deafrica_tools.dask import create_local_dask_cluster from deafrica_tools.wetlands import WIT_drill -import deafrica_tools.app.widgetconstructors as deawidgets def make_box_layout(): - return Layout( - #border='solid 1px black', - margin='0px 10px 10px 0px', - padding='5px 5px 5px 5px', - width='100%', - height='100%', - ) + 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): @@ -67,7 +48,7 @@ def create_expanded_button(description, button_style): class wit_app(HBox): def __init__(self, lang=None): super().__init__() - + deafrica_tools.set_lang(lang) ########################################################## @@ -80,30 +61,31 @@ def __init__(self, lang=None): self.out_csv = "example_WIT.csv" self.out_plot = "example_WIT.png" self.product_list = [ - (_("None"), "none"), - (_("ESRI World Imagery"), "esri_world_imagery"), - (_("Sentinel-2 Geomedian"), "gm_s2_annual"), - (_("Water Observations from Space"), "wofs_ls_summary_annual"), - + (("None"), "none"), + (("ESRI World Imagery"), "esri_world_imagery"), + (("Sentinel-2 Geomedian"), "gm_s2_annual"), + (("Water Observations from Space"), "wofs_ls_summary_annual"), ] self.product = self.product_list[0][1] self.product_year = "2020-01-01" self.target = None self.action = None self.gdf_drawn = None - + ########################################################## # HEADER FOR APP # - + # Create the Header widget - header_title_text = _("Wetlands Insight Tool") - instruction_text = _("Select parameters and AOI") - self.header = deawidgets.create_html(f"

{header_title_text}

{instruction_text}

") + header_title_text = "Wetlands Insight Tool" + instruction_text = "Select parameters and AOI" + 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): @@ -119,16 +101,22 @@ def update_geojson(target, action, geo_json): self.gdf_drawn = gdf gdf_drawn_epsg6933 = gdf.copy().to_crs("EPSG:6933") - m2_per_km2 = 10 ** 6 + m2_per_km2 = 10**6 area = gdf_drawn_epsg6933.area.values[0] / m2_per_km2 - polyarea_label = _('Total polygon area') + polyarea_label = "Total polygon area" polyarea_text = f"

{polyarea_label}: {area:.2f} km2

" if area <= 3000: - confirmation_text = '

' + _('Area falls within recommended limit') + '

' + confirmation_text = ( + '

' + ("Area falls within recommended limit") + "

" + ) self.header.value = header_title_text + polyarea_text + confirmation_text else: - warning_text = '

' + _('Area is too large, please update your polygon') + '

' + warning_text = ( + '

' + + ("Area is too large, please update your polygon") + + "

" + ) self.header.value = header_title_text + polyarea_text + warning_text ########################################################## @@ -143,22 +131,22 @@ def update_geojson(target, action, geo_json): # MAP WIDGET, DRAWING TOOLS, WMS LAYERS # # Create drawing tools - desired_drawtools = ['rectangle', 'polygon'] + desired_drawtools = ["rectangle", "polygon"] draw_control = deawidgets.create_drawcontrol(desired_drawtools) - + # Begin by displaying an empty layer group, and update the group with desired WMS on interaction. self.deafrica_layers = LayerGroup(layers=()) - self.deafrica_layers.name = _('Map Overlays') + self.deafrica_layers.name = "Map Overlays" # Create map widget self.m = deawidgets.create_map() - + self.m.layout = make_box_layout() - + # Add tools to map widget self.m.add_control(draw_control) self.m.add_layer(self.deafrica_layers) - + # Store current basemap for future use self.basemap = self.m.basemap @@ -173,26 +161,26 @@ def update_geojson(target, action, geo_json): output_csv = deawidgets.create_inputtext(self.out_csv, self.out_csv) output_plot = deawidgets.create_inputtext(self.out_plot, self.out_plot) deaoverlay_dropdown = deawidgets.create_dropdown(self.product_list, self.product_list[0][1]) - run_button = create_expanded_button(_("Run"), "info") + run_button = create_expanded_button(("Run"), "info") ########################################################## # COLLECTION OF ALL APP CONTROLS # - + parameter_selection = VBox( [ - HTML("" + _("Map Overlay:") + ""), + HTML("" + ("Map Overlay:") + ""), deaoverlay_dropdown, - HTML("" + _("Start Date:") + ""), + HTML("" + ("Start Date:") + ""), startdate_picker, - HTML("" + _("End Date:") + ""), + HTML("" + ("End Date:") + ""), enddate_picker, - HTML("" + _("Minimum Good Data:") + ""), + HTML("" + ("Minimum Good Data:") + ""), min_good_data, - HTML("" + _("Resampling Frequency:") + ""), + HTML("" + ("Resampling Frequency:") + ""), resampling_freq, - HTML("" + _("Output CSV:") + ""), + HTML("" + ("Output CSV:") + ""), output_csv, - HTML("" + _("Output Plot:") + ""), + HTML("" + ("Output Plot:") + ""), output_plot, ] ) @@ -208,7 +196,7 @@ def update_geojson(target, action, geo_json): grid[0, :] = self.header grid[1:6, 0:2] = parameter_selection grid[6, 0:2] = run_button - + # Dask and Progress info grid[1, 7:] = self.dask_client grid[2:7, 7:] = self.progress_bar @@ -280,37 +268,32 @@ def update_deaoverlay(self, change): self.deafrica_layers.add_layer(layer) def run_app(self, change): - + # Clear progress bar and output areas before running self.dask_client.clear_output() self.progress_bar.clear_output() self.wit_plot.clear_output() - # Connect to datacube database - dc = datacube.Datacube(app="wetland_app") - # Configure local dask cluster with self.dask_client: - client = create_local_dask_cluster( - return_client=True, display_client=True - ) + client = create_local_dask_cluster(return_client=True, display_client=True) # Set any defaults TCW_threshold = -0.035 dask_chunks = dict(x=1000, y=1000, time=1) - - #check resampling freq - if self.resamplingfreq == 'None': + + # check resampling freq + if self.resamplingfreq == "None": rsf = None else: rsf = self.resamplingfreq - - self.progress_header.value = f"

"+_("Progress")+"

" - + + self.progress_header.value = "

" + ("Progress") + "

" + # run wetlands polygon drill with self.progress_bar: -# with ProgressBar(): - warnings.filterwarnings("ignore") + # with ProgressBar(): + warnings.filterwarnings("ignore") try: df = WIT_drill( gdf=self.gdf_drawn, @@ -323,10 +306,10 @@ def run_app(self, change): verbose=False, verbose_progress=True, ) - print(_("WIT complete")) + print(("WIT complete")) except AttributeError: - print(_("No polygon selected")) - + print(("No polygon selected")) + # close down the dask client client.shutdown() @@ -335,7 +318,7 @@ def run_app(self, change): df.to_csv(self.out_csv, index_label="Datetime") # ---Plotting------------------------------ - + with self.wit_plot: fontsize = 17 @@ -362,11 +345,11 @@ def run_app(self, change): df.dry_veg_percent, df.bare_soil_percent, labels=[ - _("open water"), - _("wet"), - _("green veg"), - _("dry veg"), - _("bare soil"), + ("open water"), + ("wet"), + ("green veg"), + ("dry veg"), + ("bare soil"), ], colors=pal, alpha=0.6, @@ -379,7 +362,7 @@ def run_app(self, change): # add a legend and a tight plot box ax.legend(loc="lower left", framealpha=0.6) - ax.set_title(_("Percentage Fractional Cover, Wetness, and Water")) + ax.set_title(("Percentage Fractional Cover, Wetness, and Water")) # plt.tight_layout() plt.show() diff --git a/Tools/pyproject.toml b/Tools/pyproject.toml index 37323028..44d25d38 100644 --- a/Tools/pyproject.toml +++ b/Tools/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "setuptools.build_meta" [project] name = "deafrica-tools" # reflect version changes in deafrica_tools/__init__.py -version = "2.4.8" +version = "2.4.9" description = "Functions and algorithms for analysing Digital Earth Africa data." authors = [{name = "Digital Earth Africa", email = "systems@digitalearthafrica.org"}] maintainers = [{name = "Digital Earth Africa", email = "systems@digitalearthafrica.org"}] @@ -51,7 +51,7 @@ dependencies= [ "numexpr", "numpy", "odc-algo", - "odc-geo", + "odc-geo>=0.4.2", "odc-ui", "OWSLib", "packaging", @@ -65,6 +65,7 @@ dependencies= [ "rasterio", "rasterstats", "requests", + "rioxarray", "scikit-image", "scikit-learn", "scipy",