nghiadang
338 lines
10 KiB
Python
338 lines
10 KiB
Python
# 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
|
|
)
|
|
import datetime as dt
|
|
import datacube
|
|
from osgeo import ogr
|
|
import matplotlib as mpl
|
|
import matplotlib.pyplot as plt
|
|
import rasterio
|
|
from rasterio.features import geometry_mask
|
|
import xarray as xr
|
|
from IPython.display import display
|
|
import warnings
|
|
import ipywidgets as widgets
|
|
import json
|
|
import geopandas as gpd
|
|
from io import BytesIO
|
|
|
|
# Load utility functions
|
|
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):
|
|
"""
|
|
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
|
|
The central latitude to analyse
|
|
lon: float
|
|
The central longitude to analyse
|
|
buffer:
|
|
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.
|
|
The app will automatically load all data available for the two years prior to this date.
|
|
|
|
Returns
|
|
----------
|
|
ds: xarray.Dataset
|
|
data set containing combined, masked data
|
|
Masked values are set to 'nan'
|
|
"""
|
|
|
|
# Suppress warnings
|
|
warnings.filterwarnings('ignore')
|
|
|
|
# Initialise the data cube. 'app' argument is used to identify this app
|
|
dc = datacube.Datacube(app='Crophealth-app')
|
|
|
|
# Define area to load
|
|
latitude = (lat - buffer, lat + buffer)
|
|
longitude = (lon - buffer, lon + buffer)
|
|
|
|
# Specify the date range
|
|
# Calculated as today's date, subtract 730 days to collect two years of data
|
|
# Dates are converted to strings as required by loading function below
|
|
end_date = dt.datetime.strptime(date, "%Y-%m-%d")
|
|
start_date = end_date - dt.timedelta(days=730)
|
|
|
|
time = (start_date.strftime("%Y-%m-%d"), end_date.strftime("%Y-%m-%d"))
|
|
|
|
# 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)
|
|
}
|
|
|
|
# Load the data and mask out bad quality pixels
|
|
ds = load_ard(dc, products=products, min_gooddata=0.5, **query)
|
|
|
|
# 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')
|
|
|
|
# Return the data
|
|
return(ds)
|
|
|
|
|
|
def run_crophealth_app(ds, lat, lon, buffer):
|
|
"""
|
|
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
|
|
data set containing combined, masked data
|
|
Masked values are set to 'nan'
|
|
lat: float
|
|
The central latitude corresponding to the area of loaded ds
|
|
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.
|
|
For reasonable loading times, set this as `0.1` or lower.
|
|
"""
|
|
|
|
# Suppress warnings
|
|
warnings.filterwarnings('ignore')
|
|
|
|
# Update plotting functionality through rcParams
|
|
mpl.rcParams.update({'figure.autolayout': True})
|
|
|
|
# Define polygon bounds
|
|
latitude = (lat - buffer, lat + buffer)
|
|
longitude = (lon - buffer, lon + buffer)
|
|
|
|
# Define the bounding box that will be overlayed on the interactive map
|
|
# The bounds are hard-coded to match those from the loaded data
|
|
geom_obj = {
|
|
"type": "Feature",
|
|
"properties": {
|
|
"style": {
|
|
"stroke": True,
|
|
"color": 'red',
|
|
"weight": 4,
|
|
"opacity": 0.8,
|
|
"fill": True,
|
|
"fillColor": False,
|
|
"fillOpacity": 0,
|
|
"showArea": 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]
|
|
]
|
|
]
|
|
]
|
|
}
|
|
}
|
|
|
|
# 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_center = [
|
|
loadeddata_geometry.Centroid().GetY(),
|
|
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
|
|
)
|
|
|
|
# define the drawing controls
|
|
studyarea_drawctrl = DrawControl(
|
|
polygon={"shapeOptions": {"fillOpacity": 0}},
|
|
marker={},
|
|
circle={},
|
|
circlemarker={},
|
|
polyline={},
|
|
)
|
|
|
|
# add drawing controls and data bound geometry to the map
|
|
studyarea_map.add_control(studyarea_drawctrl)
|
|
studyarea_map.add_layer(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'})
|
|
with instruction:
|
|
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'})
|
|
with info:
|
|
print("Plot status:")
|
|
|
|
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']
|
|
|
|
# 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':
|
|
|
|
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')
|
|
|
|
# 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)
|
|
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(
|
|
GeoJSON(
|
|
data=geo_json,
|
|
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
|
|
)
|
|
|
|
# reset titles back to custom
|
|
ax.set_title("Average NDVI from Sentinel-2")
|
|
ax.set_xlabel("Date")
|
|
ax.set_ylabel("NDVI")
|
|
|
|
# refresh display
|
|
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.")
|
|
|
|
# Iterate the polygon number before drawing another polygon
|
|
polygon_number = polygon_number + 1
|
|
|
|
else:
|
|
info.clear_output(wait=True)
|
|
with info:
|
|
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)
|
|
|
|
with fig_display:
|
|
# TODO: update with user friendly something
|
|
display(widgets.HTML(""))
|
|
|
|
# Construct UI:
|
|
# +-----------------------+
|
|
# | instruction |
|
|
# +-----------+-----------+
|
|
# | map | plot |
|
|
# | | |
|
|
# +-----------+-----------+
|
|
# | info |
|
|
# +-----------------------+
|
|
ui = widgets.VBox([instruction,
|
|
widgets.HBox([studyarea_map, fig_display]),
|
|
info])
|
|
display(ui)
|