diff --git a/intertidal/exposure.py b/intertidal/exposure.py index de87bbd..e9e5b98 100644 --- a/intertidal/exposure.py +++ b/intertidal/exposure.py @@ -1,19 +1,39 @@ import xarray as xr import numpy as np import geopandas as gpd +import pandas as pd +from shapely.geometry import Point +from shapely.ops import unary_union +import sunriset +from math import ceil +import datetime +from datetime import timedelta +import re +import pytz +from pyproj import CRS, Transformer +from scipy.signal import argrelmax, argrelmin +from numpy import interp + +from dea_tools.coastal import pixel_tides, model_tides from intertidal.tide_modelling import pixel_tides_ensemble -from intertidal.utils import configure_logging +from intertidal.utils import configure_logging, round_date_strings def exposure( - dem, - times, - tide_model="FES2014", - tide_model_dir="/var/share/tide_models", - run_id=None, - log=None, -): + dem, + times, + start_date, + end_date, + modelled_freq = "30min", + tide_model="FES2014", + tide_model_dir="/var/share/tide_models", + filters = ['unfiltered'], + filters_combined = None, + run_id=None, + log=None, + ): + """ Calculate intertidal exposure for each pixel, indicating the proportion of time that each pixel was "exposed" from tidal @@ -21,15 +41,24 @@ def exposure( The exposure calculation is based on tide-height differences between the elevation value and modelled tide height percentiles. + + Calculate exposure percentage for each pixel based on tide-height + differences between the elevation value and percentile values of the + tide model for a given time range. Parameters ---------- dem : xarray.DataArray - xarray.DataArray containing Digital Elevation Model (DEM) data. - times : pandas.DatetimeIndex or list of pandas.Timestamps - High frequency times used to model tides across the period of - interest. These are used to evaluate how frequently each pixel - in the DEM was exposed or inundated by the tide. + xarray.DataArray containing Digital Elevation Model (DEM) data + and coordinates and attributes metadata. + time_range : tuple + Tuple containing start and end time of time range to be used for + tide model in the format of "YYYY-MM-DD". + modelled_freq : str + A pandas time offset alias for the frequency with which to + calculate the tide model during exposure calculations. Examples + include '30min' for 30 minute cadence or '1h' for a one-hourly + cadence. Defaults to '30min'. tide_model : str, optional The tide model or a list of models used to model tides, as supported by the `pyTMD` Python package. Options include: @@ -44,7 +73,12 @@ def exposure( The directory containing tide model data files. Defaults to "/var/share/tide_models"; for more information about the directory structure, refer to `dea_tools.coastal.model_tides`. - run_id : string, optional + filters : list of strings, optional + A list of customisation options to input into the tidal + modelling to calculate exposure. Defaults to ['unfiltered'] + filters_combined : list of two-object tuples, optional + Defaults to None. + run_id : string, optional An optional string giving the name of the analysis; used to prefix log entries. log : logging.Logger, optional @@ -53,12 +87,13 @@ def exposure( Returns ------- exposure : xarray.DataArray - An array containing the percentage time 'exposure' of - each pixel from tidal inundation for the duration of the modelling - period. - tide_cq : xarray.DataArray - An array containing the quantiled high temporal resolution tide - modelling for each pixel. Dimensions should be 'quantile', 'x' and 'y'. + An xarray.Dataset containing an array for each filter of + the percentage time exposure of each pixel from seawater for + the duration of the modelling period `timerange`. + tide_cq : xarray.DataArray ##Revise to dataset + An xarray.DataArray containing the quantiled high temporal + resolution tide modelling for each pixel. Dimesions should be + 'quantile', 'x' and 'y'. Notes ----- @@ -69,42 +104,756 @@ def exposure( - The 'idxmin' variable is the index of the smallest tide-height difference (i.e., maximum similarity) per pixel and is equivalent to the exposure percent. + - filters = 'unfiltered' produces exposure for the full input time + period. + - temporal filters include any of: 'dry', 'wet', 'summer', 'autumn', + 'winter', 'spring', 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', + 'Aug', 'Sep', 'Oct', 'Nov', 'Dec', 'Daylight', 'Night' + - spatial filters include any of: 'Spring_high', 'Spring_low', + 'Neap_high', 'Neap_low', 'Hightide', 'Lowtide' + - filters_combined can be any combination of one temporal and one + spatial filter + - if filters is set to `None`, no exposure will be calculated and + the program will fail unless a tuple is nominated in `filters_combined` + """ - # Set up logs if no log is passed in if log is None: log = configure_logging() # Use run ID name for logs if it exists run_id = "Processing" if run_id is None else run_id - # Create the tide-height percentiles from which to calculate # exposure statistics - tide_percentiles = np.linspace(0, 1, 101) + calculate_quantiles = np.linspace(0, 1, 101) #nb formerly 'pc_range' + + # Generate range of times covering entire period of satellite record for exposure and bias/offset calculation + time_range = pd.date_range( + start=round_date_strings(start_date, round_type="start"), + end=round_date_strings(end_date, round_type="end"), + freq=modelled_freq, + ) + # Separate 'filters' into spatial and temporal categories to define + # which exposure workflow to use + temporal_filters = ['dry', 'wet', 'summer', 'autumn', 'winter', 'spring', 'Jan', 'Feb', 'Mar', 'Apr', + 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec', 'Daylight', 'Night'] + spatial_filters = ['Spring_high', 'Spring_low', 'Neap_high', 'Neap_low', 'Hightide', 'Lowtide'] + + ## Set the required range of tide-height percentiles for exposure calculation + # calculate_quantiles = np.linspace(0, 1, 101) + calculate_quantiles = np.linspace(0, 1, 1001)#Temporary to separate exposure values + + ## Create empty datasets to store outputs into + exposure = xr.Dataset(coords=dict(y=(['y'], dem.y.values), + x=(['x'], dem.x.values))) + tide_cq_dict = xr.Dataset(coords=dict(y=(['y'], dem.y.values), + x=(['x'], dem.x.values))) + + ## Create an empty dict to store temporal `time_range` variables into + timeranges = {} + + ## If filter combinations are desired, make sure each filter is calculated individually for later combination + if filters_combined is not None: + for x in filters_combined: + filters.append(str(x[0])) if str(x[0]) not in filters else next + filters.append(str(x[1])) if x[1] not in filters else next - # Run the `pixel_tides_ensemble` function with the `calculate_quantiles` - # option. For each pixel, an array of tide heights is returned, - # corresponding to the percentiles of the tide range specified by - # `tide_percentiles`. - log.info(f"{run_id}: Modelling tide heights for each pixel") - tide_cq, _ = pixel_tides_ensemble( - ds=dem, - ancillary_points="data/raw/tide_correlations_2017-2019.geojson", - model=tide_model, - directory=tide_model_dir, - calculate_quantiles=tide_percentiles, - times=times, - ) - - # Calculate the tide-height difference between the elevation value and - # each percentile value per pixel - diff = abs(tide_cq - dem) - - # Take the percentile of the smallest tide-height difference as the - # exposure % per pixel - idxmin = diff.idxmin(dim="quantile") - - # Convert to percentage - exposure = idxmin * 100 - - return exposure, tide_cq + if 'unfiltered' in filters: + print('Calculating unfiltered exposure') + + if (tide_model[0] == "ensemble") or (tide_model == "ensemble"): + # Use ensemble model combining multiple input ocean tide models + tide_cq, _ = pixel_tides_ensemble( + dem, + calculate_quantiles=calculate_quantiles, + times=time_range, + directory=tide_model_dir, + ancillary_points="data/raw/tide_correlations_2017-2019.geojson", + top_n=3, + reduce_method='mean', + resolution=3000, + ) + + else: + # Use single input ocean tide model + tide_cq, _ = pixel_tides( + dem, + resample=True, + calculate_quantiles=calculate_quantiles, + times=time_range, + model=tide_model, + directory=tide_model_dir, + ) + + # Add tide_cq to output dict + tide_cq_dict['unfiltered']=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['unfiltered'] = idxmin * 100 + + # return exposure + + if any (x in spatial_filters for x in filters): + if (tide_model[0] == "ensemble") or (tide_model == "ensemble"): + # Use ensemble model combining multiple input ocean tide models + ModelledTides, _ = pixel_tides_ensemble( + dem, + calculate_quantiles=calculate_quantiles, + times=time_range, + directory=tide_model_dir, + ancillary_points="data/raw/tide_correlations_2017-2019.geojson", + top_n=3, + reduce_method='mean', + resolution=3000, + ) + + else: + # Use single input ocean tide model + ModelledTides, _ = pixel_tides( + dem, + calculate_quantiles=calculate_quantiles, + times=time_range, + resample=True, + model=tide_model, + directory=tide_model_dir, + ) + + ## For use with spatial filter options + ## stack the y and x dimensions + stacked_everything = ModelledTides.stack(z=['y','x']).groupby('z') + + # Filter the input timerange to include only dates or tide ranges of interest + # if filters is not None: + for x in filters: + if x in temporal_filters: + print(f'Calculating temporal filter: {x}') + + if x == 'dry': + timeranges['dry'] = time_range.drop(time_range[(time_range.month == 10) ## Wet season: Oct-Mar + |(time_range.month == 11) + |(time_range.month == 12) + |(time_range.month == 1) + |(time_range.month == 2) + |(time_range.month == 3) + ]) + elif x == 'wet': + timeranges['wet'] = time_range.drop(time_range[(time_range.month == 4) ## Dry season: Apr-Sep + |(time_range.month == 5) + |(time_range.month == 6) + |(time_range.month == 7) + |(time_range.month == 8) + |(time_range.month == 9) + ]) + elif x == 'summer': + timeranges['summer'] = time_range.drop(time_range[time_range.quarter != 1]) + elif x == 'autumn': + timeranges['autumn'] = time_range.drop(time_range[time_range.quarter != 2]) + elif x == 'winter': + timeranges['winter'] = time_range.drop(time_range[time_range.quarter != 3]) + elif x == 'spring': + timeranges['spring'] = time_range.drop(time_range[time_range.quarter != 4]) + elif x == 'Jan': + timeranges['Jan'] = time_range.drop(time_range[time_range.month != 1]) + elif x == 'Feb': + timeranges['Feb'] = time_range.drop(time_range[time_range.month != 2]) + elif x == 'Mar': + timeranges['Mar'] = time_range.drop(time_range[time_range.month != 3]) + elif x == 'Apr': + timeranges['Apr'] = time_range.drop(time_range[time_range.month != 4]) + elif x == 'May': + timeranges['May'] = time_range.drop(time_range[time_range.month != 5]) + elif x == 'Jun': + timeranges['Jun'] = time_range.drop(time_range[time_range.month != 6]) + elif x == 'Jul': + timeranges['Jul'] = time_range.drop(time_range[time_range.month != 7]) + elif x == 'Aug': + timeranges['Aug'] = time_range.drop(time_range[time_range.month != 8]) + elif x == 'Sep': + timeranges['Sep'] = time_range.drop(time_range[time_range.month != 9]) + elif x == 'Oct': + timeranges['Oct'] = time_range.drop(time_range[time_range.month != 10]) + elif x == 'Nov': + timeranges['Nov'] = time_range.drop(time_range[time_range.month != 11]) + elif x == 'Dec': + timeranges['Dec'] = time_range.drop(time_range[time_range.month != 12]) + elif x == 'Daylight' or 'Night': + + timezones = {'wa':'../../gdata1/data/boundaries/GEODATA_COAST_100K/western_australia/cstwacd_r.shp', + 'nt':'../../gdata1/data/boundaries/GEODATA_COAST_100K/northern_territory/cstntcd_r.shp', + 'sa':'../../gdata1/data/boundaries/GEODATA_COAST_100K/south_australia/cstsacd_r.shp', + 'qld':'../../gdata1/data/boundaries/GEODATA_COAST_100K/queensland/cstqldmd_r.shp', + 'nsw':'../../gdata1/data/boundaries/GEODATA_COAST_100K/new_south_wales/cstnswcd_r.shp', + 'vic':'../../gdata1/data/boundaries/GEODATA_COAST_100K/victoria/cstviccd_r.shp', + 'tas':'../../gdata1/data/boundaries/GEODATA_COAST_100K/tasmania/csttascd_r.shp' + } + + ## Determine the timezone of the pixels + ## Bring in the state polygons (note: native crs = epsg:4283) + wa = gpd.read_file(timezones['wa']) + nt = gpd.read_file(timezones['nt']) + sa = gpd.read_file(timezones['sa']) + qld = gpd.read_file(timezones['qld']) + nsw = gpd.read_file(timezones['nsw']) + vic = gpd.read_file(timezones['vic']) + tas = gpd.read_file(timezones['tas']) + + # Merge the polygons to create single state/territory boundaries + wa = gpd.GeoSeries(unary_union(wa.geometry)) + nt = gpd.GeoSeries(unary_union(nt.geometry)) + sa = gpd.GeoSeries(unary_union(sa.geometry)) + qld = gpd.GeoSeries(unary_union(qld.geometry)) + nsw = gpd.GeoSeries(unary_union(nsw.geometry)) + vic = gpd.GeoSeries(unary_union(vic.geometry)) + tas = gpd.GeoSeries(unary_union(tas.geometry)) + + ## Note: day and night times will be calculated once per area-of-interest(ds) + ## for the median latitude and longitude position of the aoi + tidepost_lat_3577 = dem.x.median(dim='x').values + tidepost_lon_3577 = dem.y.median(dim='y').values + + ## Set the Datacube native crs (GDA/Aus Albers (meters)) + crs_3577 = CRS.from_epsg(3577) + + ## Translate the crs of the tidepost to determine (1) local timezone + ## and (2) the local sunrise and sunset times: + + ## (1) Create a transform to convert default epsg3577 coords to epsg4283 to compare + ## against state/territory boundary polygons and assign a timezone + + ## GDA94 CRS (degrees) + crs_4283 = CRS.from_epsg(4283) + ## Transfer coords from/to + transformer_4283 = Transformer.from_crs(crs_3577, crs_4283) + ## Translate tidepost coords + tidepost_lat_4283, tidepost_lon_4283 = transformer_4283.transform(tidepost_lat_3577, + tidepost_lon_3577) + ## Coordinate point to test for timezone + point_4283 = Point(tidepost_lon_4283, tidepost_lat_4283) + + ## (2) Create a transform to convert default epsg3577 coords to epsg4326 for use in + ## sunise/sunset library + + ## World WGS84 (degrees) + crs_4326 = CRS.from_epsg(4326) + ## Transfer coords from/to + transformer_4326 = Transformer.from_crs(crs_3577, crs_4326) + ## Translate the tidepost coords + tidepost_lat_4326, tidepost_lon_4326 = transformer_4326.transform(tidepost_lat_3577, + tidepost_lon_3577) + ## Coordinate point to locate the sunriset calculation + point_4326 = Point(tidepost_lon_4326, tidepost_lat_4326) + + ## Set the local timezone for the analysis area of interest + if wa.contains(point_4283)[0] == True: + timezone = 'Australia/West' + local_tz = 8 + + elif nt.contains(point_4283)[0] == True: + timezone = 'Australia/North' + local_tz = 9.5 + + elif sa.contains(point_4283)[0] == True: + timezone = 'Australia/South' + local_tz = 9.5 + + elif qld.contains(point_4283)[0] == True: + timezone = 'Australia/Queensland' + local_tz = 10 + + elif nsw.contains(point_4283)[0] == True: + timezone = 'Australia/NSW' + local_tz = 10 + + elif vic.contains(point_4283)[0] == True: + timezone = 'Australia/Victoria' + local_tz = 10 + + elif tas.contains(point_4283)[0] == True: + timezone = 'Australia/Tasmania' + local_tz = 10 + + ## Calculate the local sunrise and sunset times + # Place start and end dates in correct format + start = time_range[0] + end = time_range[-1] + startdate = datetime.date(pd.to_datetime(start).year, + pd.to_datetime(start).month, + pd.to_datetime(start).day) + + # Make 'all_timerange' time-zone aware + localtides = time_range.tz_localize(tz=pytz.UTC).tz_convert(timezone) + + # Replace the UTC datetimes from all_timerange with local times + modelledtides = pd.DataFrame(index = localtides) + + # Return the difference in years for the time-period. + # Round up to ensure all modelledtide datetimes are captured in the solar model + diff = pd.to_datetime(end) - pd.to_datetime(start) + diff = int(ceil(diff.days/365)) + + ## Model sunrise and sunset + sun_df = sunriset.to_pandas(startdate, tidepost_lat_4326, tidepost_lon_4326, local_tz, diff) + + ## Set the index as a datetimeindex to match the modelledtide df + sun_df = sun_df.set_index(pd.DatetimeIndex(sun_df.index)) + + ## Append the date to each Sunrise and Sunset time + sun_df['Sunrise dt'] = sun_df.index + sun_df['Sunrise'] + sun_df['Sunset dt'] = sun_df.index + (sun_df['Sunset']) + + ## Create new dataframes where daytime and nightime datetimes are recorded, then merged + ## on a new `Sunlight` column + daytime=pd.DataFrame(data = 'Sunrise', index=sun_df['Sunrise dt'], columns=['Sunlight']) + nighttime=pd.DataFrame(data = 'Sunset', index=sun_df['Sunset dt'], columns=['Sunlight']) + DayNight = pd.concat([daytime, nighttime], join='outer') + DayNight.sort_index(inplace=True) + DayNight.index.rename('Datetime', inplace=True) + + ## Create an xarray object from the merged day/night dataframe + day_night = xr.Dataset.from_dataframe(DayNight) + + ## Remove local timezone timestamp column in modelledtides dataframe. Xarray doesn't handle + ## timezone aware datetimeindexes 'from_dataframe' very well. + modelledtides.index = modelledtides.index.tz_localize(tz=None) + + ## Create an xr Dataset from the modelledtides pd.dataframe + mt = modelledtides.to_xarray() + + ## Filter the modelledtides (mt) by the daytime, nighttime datetimes from the sunriset module + ## Modelled tides are designated as either day or night by propogation of the last valid index + ## value forward + Solar=day_night.sel(Datetime=mt.index, method='ffill') + + ## Assign the day and night tideheight datasets + SolarDayTides = mt.where(Solar.Sunlight=='Sunrise', drop=True) + SolarNightTides = mt.where(Solar.Sunlight=='Sunset', drop=True) + + ## Extract DatetimeIndexes to use in exposure calculations + all_timerange_day = pd.DatetimeIndex(SolarDayTides.index) + all_timerange_night = pd.DatetimeIndex(SolarNightTides.index) + + if x == 'Daylight': + timeranges['Daylight'] = all_timerange_day + if x == 'Night': + timeranges['Night'] = all_timerange_night + + elif x in spatial_filters: + print(f'Calculating statial filter: {x}') + + # # Extract the modelling freq units + # Split the number and text characters in modelled_freq + freq_time = int(re.findall(r'(\d+)(\w+)', modelled_freq)[0][0]) + freq_unit = str(re.findall(r'(\d+)(\w+)', modelled_freq)[0][-1]) + + # Extract the number of modelled timesteps per 14 days (half lunar cycle) for neap/spring calcs + mod_timesteps = pd.Timedelta(14,"d")/pd.Timedelta(freq_time, freq_unit) + + ## Identify kwargs for peak detection algorithm + order=(int(mod_timesteps/2)) + + ## Calculate the spring highest and spring lowest tides per 14 day half lunar cycle + if x == 'Spring_high': + + print ('Calculating Spring_high') + + ## apply the peak detection routine + stacked_everything_high = stacked_everything.apply( + lambda x: xr.DataArray(argrelmax(x.values, order=order)[0]) + ) + ## Unstack + springhighs_all = stacked_everything_high.unstack('z') + ##Reorder the y axis. Uncertain why it gets reversed during the stack/unstack. + springhighs_all = springhighs_all.reindex(y=springhighs_all.y[::-1]) + ## Rename the time axis + springhighs_all = springhighs_all.rename({'dim_0':'time'}) + ## Convert to dataset + springhighs_all = springhighs_all.to_dataset(name = 'time') + ## Reorder the dims + springhighs_all = springhighs_all[['time','y','x']] + + # Select dates associated with detected peaks + ## removed reference below to ModelledTides[0]. Possibly an artefact of new + ## pixel_tides_ensemble func. If using pixel_tides, may need to revert to ModelledTides[0]. + + # springhighs_all = ModelledTides[0].to_dataset().isel(time=springhighs_all.time) + springhighs_all = ModelledTides.to_dataset().isel(time=springhighs_all.time) + + ## Save datetimes for calculation of combined filter exposure + timeranges['Spring_high'] = pd.to_datetime(springhighs_all.isel(x=1,y=1).time) + + ## Extract the peak height dates + tide_cq = springhighs_all.tide_m.quantile(q=calculate_quantiles,dim='time') + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['Spring_high'] = idxmin * 100 + + + ## Calculate the spring highest and spring lowest tides per 14 day half lunar cycle + if x == 'Spring_low': + print ('Calculating Spring_low') + + ## apply the peak detection routine + stacked_everything_low = stacked_everything.apply(lambda x: xr.DataArray(argrelmin(x.values, order=order)[0])) + ## Unstack + springlows_all = stacked_everything_low.unstack('z') + ##Reorder the y axis. Uncertain why it gets reversed during the stack/unstack. + springlows_all = springlows_all.reindex(y=springlows_all.y[::-1]) + ## Rename the time axis + springlows_all = springlows_all.rename({'dim_0':'time'}) + ## Convert to dataset + springlows_all = springlows_all.to_dataset(name = 'time') + ## Reorder the dims + springlows_all = springlows_all[['time','y','x']] + ## Select dates associated with detected peaks + # springlows_all = ModelledTides[0].to_dataset().isel(time=springlows_all.time) + springlows_all = ModelledTides.to_dataset().isel(time=springlows_all.time)## removed reference to ModelledTides[0]. Possibly an artefact of new pixel_tides_ensemble func. If using pixel_tides, may need to revert to ModelledTides[0]. + + ## Save datetimes for calculation of combined filter exposure + timeranges['Spring_low'] = pd.to_datetime(springlows_all.isel(x=1,y=1).time) + + tide_cq = springlows_all.tide_m.quantile(q=calculate_quantiles,dim='time') + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['Spring_low'] = idxmin * 100 + + if x == 'Neap_high': + print ('Calculating Neap_high') + ## Calculate the number of spring high tides to support calculation of neap highs + ## apply the peak detection routine + stacked_everything_high = stacked_everything.apply(lambda x: xr.DataArray(argrelmax(x.values, order=order)[0])) + ## Unstack + springhighs_all = stacked_everything_high.unstack('z') + + ## apply the peak detection routine to calculate all the high tide maxima + Max_testarray = stacked_everything.apply(lambda x: xr.DataArray(argrelmax(x.values)[0])) + ## extract the corresponding dates from the peaks + Max_testarray = (Max_testarray.unstack('z')) + Max_testarray = (Max_testarray.reindex(y=Max_testarray.y[::-1]) + .rename({'dim_0':'time'}) + .to_dataset(name = 'time') + [['time','y','x']] + ) + ## extract all hightide peaks + # Max_testarray = ModelledTides[0].to_dataset().isel(time=Max_testarray.time) + Max_testarray = ModelledTides.to_dataset().isel(time=Max_testarray.time)## removed reference to ModelledTides[0]. Possibly an artefact of new pixel_tides_ensemble func. If using pixel_tides, may need to revert to ModelledTides[0]. + + ## repeat the peak detection to identify neap high tides (minima in the high tide maxima) + stacked_everything2 = Max_testarray.tide_m.stack(z=['y','x']).groupby('z') + ## extract neap high tides based on 14 day half lunar cycle - determined as the fraction of all high tide points + ## relative to the number of spring high tide values + order_nh = int(ceil((len(Max_testarray.time)/(len(springhighs_all))/2))) + ## apply the peak detection routine to calculate all the neap high tide minima within the high tide peaks + neaphighs_all = stacked_everything2.apply(lambda x: xr.DataArray(argrelmin(x.values, order=order_nh)[0])) + ## unstack and format as above + neaphighs_all = neaphighs_all.unstack('z') + neaphighs_all = ( + neaphighs_all + .reindex(y=neaphighs_all.y[::-1]) + .rename({'dim_0':'time'}) + .to_dataset(name = 'time') + [['time','y','x']] + ) + ## extract neap high tides + neaphighs_all = Max_testarray.isel(time=neaphighs_all.time) + + ## Save datetimes for calculation of combined filter exposure + timeranges['Neap_high'] = pd.to_datetime(neaphighs_all.isel(x=1,y=1).time) + + tide_cq = neaphighs_all.tide_m.quantile(q=calculate_quantiles,dim='time') + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['Neap_high'] = idxmin * 100 + + if x == 'Neap_low': + print ('Calculating Neap_low') + ## Calculate the number of spring low tides to support calculation of neap low tides + ## apply the peak detection routine + stacked_everything_low = stacked_everything.apply(lambda x: xr.DataArray(argrelmin(x.values, order=order)[0])) + + ## Unstack + springlows_all = stacked_everything_low.unstack('z') + + ## apply the peak detection routine to calculate all the low tide maxima + Min_testarray = stacked_everything.apply(lambda x: xr.DataArray(argrelmin(x.values)[0])) + + ## extract the corresponding dates from the peaks + Min_testarray = (Min_testarray.unstack('z')) + Min_testarray = (Min_testarray.reindex(y=Min_testarray.y[::-1]) + .rename({'dim_0':'time'}) + .to_dataset(name = 'time') + [['time','y','x']] + ) + + ## extract all lowtide peaks + # Min_testarray = ModelledTides[0].to_dataset().isel(time=Min_testarray.time) + Max_testarray = ModelledTides.to_dataset().isel(time=Max_testarray.time)## removed reference to ModelledTides[0]. Possibly an artefact of new pixel_tides_ensemble func. If using pixel_tides, may need to revert to ModelledTides[0]. + + ## repeat the peak detection to identify neap low tides (maxima in the low tide maxima) + stacked_everything2 = Min_testarray.tide_m.stack(z=['y','x']).groupby('z') + + ## extract neap high tides based on 14 day half lunar cycle - determined as the fraction of all high tide points + ## relative to the number of spring high tide values + order_nl = int(ceil((len(Min_testarray.time)/(len(springlows_all))/2))) + + ## apply the peak detection routine to calculate all the neap high tide minima within the high tide peaks + neaplows_all = stacked_everything2.apply(lambda x: xr.DataArray(argrelmax(x.values, order=order_nl)[0])) + + ## unstack and format as above + neaplows_all = neaplows_all.unstack('z') + neaplows_all = ( + neaplows_all + .reindex(y=neaplows_all.y[::-1]) + .rename({'dim_0':'time'}) + .to_dataset(name = 'time') + [['time','y','x']] + ) + + ## extract neap high tides + neaplows_all = Min_testarray.isel(time=neaplows_all.time) + + ## Save datetimes for calculation of combined filter exposure + timeranges['Neap_low'] = pd.to_datetime(neaplows_all.isel(x=1,y=1).time) + + tide_cq = neaplows_all.tide_m.quantile(q=calculate_quantiles,dim='time') + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['Neap_low'] = idxmin * 100 + + + if x == 'Hightide': + print ('Calculating Hightide') + def lowesthightides(x): + ''' + x is a grouping of x and y pixels from the peaks_array (labelled as 'z') + ''' + + ## apply the peak detection routine to calculate all the high tide maxima + high_peaks = np.array(argrelmax(x.values)[0]) + + ## extract all hightide peaks + Max_testarray = x.isel(time=high_peaks) + + ## Identify all lower hightide peaks + lowhigh_peaks = np.array(argrelmin(Max_testarray.values)[0]) + + ## Interpolate the lower hightide curve + neap_high_linear = interp( + ## Create an array to interpolate into + np.arange(0,len(x.time)), + ## low high peaks as a subset of all high tide peaks + high_peaks[lowhigh_peaks], + ## Corresponding tide heights + Max_testarray.isel(time=lowhigh_peaks).squeeze(['z']).values, + ) + + # # Extract hightides as all tides higher than/equal to the extrapolated lowest high tide line + hightide = x.squeeze(['z']).where(x.squeeze(['z']) >= neap_high_linear, drop=True) + + return hightide + + ## Vectorise the hightide calculation + lowhighs_all = stacked_everything.apply(lambda x: xr.DataArray(lowesthightides(x))) + + # ## Unstack and re-format the array + lowhighs_all = lowhighs_all.unstack('z') + lowhighs_all_unstacked = ( + lowhighs_all + .reindex(y=lowhighs_all.y[::-1]) + .to_dataset() + [['tide_m','time','y','x']] + ) + + ## Save datetimes for calculation of combined filter exposure + timeranges['Hightide'] = pd.to_datetime(lowhighs_all_unstacked.isel(x=1,y=1).time) + + tide_cq = lowhighs_all_unstacked.tide_m.quantile(q=calculate_quantiles,dim='time') + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['Hightide'] = idxmin * 100 + + if x == 'Lowtide': + print ('Calculating Lowtide') + def highestlowtides(x): + ''' + x is a grouping of x and y pixels from the peaks_array (labelled as 'z') + ''' + + ## apply the peak detection routine to calculate all the high tide maxima + low_peaks = np.array(argrelmin(x.values)[0]) + + ## extract all hightide peaks + Min_testarray = x.isel(time=low_peaks) + + ## Identify all lower hightide peaks + highlow_peaks = np.array(argrelmax(Min_testarray.values)[0]) + + ## Interpolate the lower hightide curve + neap_low_linear = interp( + ## Create an array to interpolate into + np.arange(0,len(x.time)), + ## low high peaks as a subset of all high tide peaks + low_peaks[highlow_peaks], + ## Corresponding tide heights + Min_testarray.isel(time=highlow_peaks).squeeze(['z']).values, + ) + + # # Extract hightides as all tides higher than/equal to the extrapolated lowest high tide line + lowtide = x.squeeze(['z']).where(x.squeeze(['z']) <= neap_low_linear, drop=True) + + return lowtide + + ## Vectorise the lowtide calculation + highlows_all = stacked_everything.apply(lambda x: xr.DataArray(highestlowtides(x))) + + # ## Unstack and re-format the array + highlows_all = highlows_all.unstack('z') + highlows_all_unstacked = ( + highlows_all + .reindex(y=highlows_all.y[::-1]) + .to_dataset() + [['tide_m','time','y','x']] + ) + + ## Save datetimes for calculation of combined filter exposure + timeranges['Lowtide'] = pd.to_datetime(highlows_all.isel(x=1,y=1).time) + + tide_cq = highlows_all_unstacked.tide_m.quantile(q=calculate_quantiles,dim='time') + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure['Lowtide'] = idxmin * 100 + + ## Intersect the filters of interest to extract the common datetimes for calc of combined filters + if filters_combined is not None: + for x in filters_combined: + y=x[0] + z=x[1] + timeranges[str(y+"_"+z)] = timeranges[y].intersection(timeranges[z]) + + ## Generator expression to calculate exposure for each nominated filter in temporal_filters + # Don't calculate exposure for spatial filters. This has already been calculated. + gen = (x for x in timeranges if x not in spatial_filters) + + for x in gen: + # Run the pixel_tides function with the calculate_quantiles option. + # For each pixel, an array of tideheights is returned, corresponding + # to the percentiles from `calculate_quantiles` of the timerange-tide model that + # each tideheight appears in the model. + + if (tide_model[0] == "ensemble") or (tide_model == "ensemble"): + # Use ensemble model combining multiple input ocean tide models + tide_cq, _ = pixel_tides_ensemble( + dem, + calculate_quantiles=calculate_quantiles, + times=timeranges[str(x)], + directory=tide_model_dir, + ancillary_points="data/raw/tide_correlations_2017-2019.geojson", + top_n=3, + reduce_method='mean', + resolution=3000, + ) + + else: + # Use single input ocean tide model + tide_cq, _ = pixel_tides_ensemble( + dem, + resample=True, + calculate_quantiles=calculate_quantiles, + times=timeranges[str(x)], + model=tide_model, + directory=tide_model_dir, + ) + + # Add tide_cq to output dict + tide_cq_dict[str(x)]=tide_cq + + # Calculate the tide-height difference between the elevation value and + # each percentile value per pixel + diff = abs(tide_cq - dem) + + # Take the percentile of the smallest tide-height difference as the + # exposure % per pixel + idxmin = diff.idxmin(dim="quantile") + + # Convert to percentage + exposure[str(x)] = idxmin * 100 + + + return exposure, tide_cq_dict + + + + diff --git a/notebooks/Intertidal_workflow.ipynb b/notebooks/Intertidal_workflow.ipynb index bf57ab3..26d92a7 100644 --- a/notebooks/Intertidal_workflow.ipynb +++ b/notebooks/Intertidal_workflow.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "44932e9f", + "id": "a1eed1f8", "metadata": { "tags": [] }, @@ -16,7 +16,7 @@ { "cell_type": "code", "execution_count": 1, - "id": "b18166da", + "id": "7f782832", "metadata": { "tags": [] }, @@ -25,7 +25,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "/home/jovyan/Robbi/dea-intertidal\n" + "/home/jovyan/dea_intertidal/dea-intertidal\n" ] } ], @@ -35,7 +35,7 @@ }, { "cell_type": "markdown", - "id": "51d593a1", + "id": "cdb87709", "metadata": {}, "source": [ "Install additional packages directly from the requirements file" @@ -43,19 +43,27 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "a9e5e54a", + "execution_count": 2, + "id": "adce236f", "metadata": { "tags": [] }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Note: you may need to restart the kernel to use updated packages.\n" + ] + } + ], "source": [ "pip install -r requirements.in --quiet" ] }, { "cell_type": "markdown", - "id": "cff69d46", + "id": "fb6f5bac", "metadata": {}, "source": [ "### Load packages" @@ -63,8 +71,8 @@ }, { "cell_type": "code", - "execution_count": 2, - "id": "06f536eb", + "execution_count": 3, + "id": "a723fcdc", "metadata": { "tags": [] }, @@ -107,7 +115,7 @@ }, { "cell_type": "markdown", - "id": "34f2a478", + "id": "55347d87", "metadata": { "tags": [] }, @@ -117,7 +125,7 @@ }, { "cell_type": "markdown", - "id": "a83b7160", + "id": "c3043ff7", "metadata": { "tags": [] }, @@ -128,7 +136,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "f00107c2", + "id": "1aa6c76d", "metadata": { "tags": [] }, @@ -144,35 +152,33 @@ "include_ls = True # Include Landsat data in the analysis?\n", "filter_gqa = True # Filter to remove poorly georeferenced scenes?\n", "tide_model = \"ensemble\" # Tide model to use in analysis\n", - "# tide_model_dir = \"/var/share/tide_models\" # Directory containing tide model files\n", - "# tide_model = [\"FES2014\", \"FES2012\", \"TPXO9-atlas-v5\"]\n", "tide_model_dir = \"/gdata1/data/tide_models_clipped\"\n", "\n", "# Exposure variables\n", "modelled_freq = \"3h\" # Frequency to run tidal model e.g '30min' or '1h'\n", + "# filters=[None] # Exposure filters.\n", + "# filters_combined = None ## Must be a list of tuples containing one temporal and spatial filter each\n", "\n", - "# Generate range of times covering entire period of satellite record for exposure and bias/offset calculation\n", - "all_times = pd.date_range(\n", - " start=round_date_strings(start_date, round_type=\"start\"),\n", - " end=round_date_strings(end_date, round_type=\"end\"),\n", - " freq=modelled_freq,\n", - ")" + "# # Generate range of times covering entire period of satellite record for exposure and bias/offset calculation\n", + "# all_timerange = pd.date_range(\n", + "# start=round_date_strings(start_date, round_type=\"start\"),\n", + "# end=round_date_strings(end_date, round_type=\"end\"),\n", + "# freq=modelled_freq,\n", + "# )" ] }, { "cell_type": "markdown", - "id": "65b1e0b2", + "id": "a7321065", "metadata": {}, "source": [ - "#### Set study area\n", - "\n", "##### Option 1: load study area from 32 km tile GridSpec" ] }, { "cell_type": "code", "execution_count": null, - "id": "a54c5653", + "id": "4caa3bd3", "metadata": { "tags": [] }, @@ -185,7 +191,7 @@ }, { "cell_type": "markdown", - "id": "cbba20c2", + "id": "5d7ccd92", "metadata": { "tags": [] }, @@ -196,7 +202,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e35efaea", + "id": "49b4231b", "metadata": { "tags": [] }, @@ -215,7 +221,7 @@ }, { "cell_type": "markdown", - "id": "cc7c76f2", + "id": "0a484060", "metadata": {}, "source": [ "##### Option 3: load study area using interactive map" @@ -224,11 +230,52 @@ { "cell_type": "code", "execution_count": 5, - "id": "9895b5c5", + "id": "79d8aba4", "metadata": { "tags": [] }, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<<<<<<< local\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3cb6c9c264bb46ca9861d181a255638c", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Map(center=[-26, 135], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_out_t…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "" + ], + "text/plain": [ + "Geometry(POLYGON ((116.547775 -20.585843, 116.547775 -20.566558, 116.566057 -20.566558, 116.566057 -20.585843, 116.547775 -20.585843)), EPSG:4326)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n" + ] + }, { "data": { "application/vnd.jupyter.widget-view+json": { @@ -255,6 +302,26 @@ "execution_count": 5, "metadata": {}, "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ">>>>>>> remote\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "" + ], + "text/plain": [ + "Geometry(POLYGON ((131.853962 -12.233161, 131.853962 -12.194405, 131.90546 -12.194405, 131.90546 -12.233161, 131.853962 -12.233161)), EPSG:4326)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -269,7 +336,7 @@ }, { "cell_type": "markdown", - "id": "3a77a221", + "id": "41a23d82", "metadata": {}, "source": [ "## Intertidal workflow\n", @@ -280,11 +347,205 @@ { "cell_type": "code", "execution_count": 6, - "id": "7b142abf", + "id": "611834aa", "metadata": { "tags": [] }, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<<<<<<< local \n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "
\n", + "
\n", + "

Client

\n", + "

Client-00d37937-c0b9-11ee-9b20-8a70e1ff1440

\n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "
Connection method: Cluster objectCluster type: distributed.LocalCluster
\n", + " Dashboard: /user/claire.phillips@ga.gov.au/proxy/8787/status\n", + "
\n", + "\n", + " \n", + " \n", + " \n", + "\n", + " \n", + "
\n", + "

Cluster Info

\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

LocalCluster

\n", + "

7a5c7c97

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + "
\n", + " Dashboard: /user/claire.phillips@ga.gov.au/proxy/8787/status\n", + " \n", + " Workers: 1\n", + "
\n", + " Total threads: 15\n", + " \n", + " Total memory: 117.21 GiB\n", + "
Status: runningUsing processes: True
\n", + "\n", + "
\n", + " \n", + "

Scheduler Info

\n", + " \n", + "\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

Scheduler

\n", + "

Scheduler-5405bd62-7764-4420-9818-871c8ffb897a

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Comm: tcp://127.0.0.1:37581\n", + " \n", + " Workers: 1\n", + "
\n", + " Dashboard: /user/claire.phillips@ga.gov.au/proxy/8787/status\n", + " \n", + " Total threads: 15\n", + "
\n", + " Started: Just now\n", + " \n", + " Total memory: 117.21 GiB\n", + "
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "

Workers

\n", + " \n", + "\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 0

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:46491\n", + " \n", + " Total threads: 15\n", + "
\n", + " Dashboard: /user/claire.phillips@ga.gov.au/proxy/42263/status\n", + " \n", + " Memory: 117.21 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:36407\n", + "
\n", + " Local directory: /tmp/dask-worker-space/worker-vs62yglb\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "\n", + "
\n", + "
\n", + "\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "\n", + "
\n", + "
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n" + ] + }, { "data": { "text/html": [ @@ -436,69 +697,596 @@ " \n", " \n", "\n", - " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + "\n", + " \n", + " \n", + "\n", + " \n", + " \n", + "\n", + " \n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ">>>>>>> remote \n", + "<<<<<<< local \n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Dimensions: (time: 325, y: 237, x: 215)\n", + "Coordinates:\n", + " * time (time) datetime64[ns] 2020-01-01T01:55:31.756059 ... 2022-12...\n", + " * y (y) float64 -2.304e+06 -2.304e+06 ... -2.306e+06 -2.306e+06\n", + " * x (x) float64 -1.598e+06 -1.598e+06 ... -1.596e+06 -1.596e+06\n", + " spatial_ref int32 3577\n", + "Data variables:\n", + " ndwi (time, y, x) float32 dask.array\n", + "Attributes:\n", + " crs: EPSG:3577\n", + " grid_mapping: spatial_ref\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Dimensions: (time: 272, y: 278, x: 289)\n", + "Coordinates:\n", + " * time (time) datetime64[ns] 2020-01-02T00:57:27.736492 ... 2022-12...\n", + " * y (y) float64 -1.286e+06 -1.286e+06 ... -1.288e+06 -1.288e+06\n", + " * x (x) float64 -1.408e+04 -1.408e+04 ... -1.122e+04 -1.12e+04\n", + " spatial_ref int32 3577\n", + "Data variables:\n", + " ndwi (time, y, x) float32 dask.array\n", + "Attributes:\n", + " crs: EPSG:3577\n", + " grid_mapping: spatial_ref\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ">>>>>>> remote \n", + "/env/lib/python3.10/site-packages/datacube/drivers/driver_cache.py:54: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html\n", + " from pkg_resources import iter_entry_points\n", + "/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.\n", + " _reproject(\n", + "<<<<<<< local \n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 5.23 s, sys: 504 ms, total: 5.73 s\n", + "Wall time: 37.1 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 9.18 s, sys: 720 ms, total: 9.9 s\n", + "Wall time: 52.2 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ">>>>>>> remote \n", + "<<<<<<< local\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (time: 325, y: 237, x: 215)\n",
+       "Coordinates:\n",
+       "  * time         (time) datetime64[ns] 2020-01-01T01:55:31.756059 ... 2022-12...\n",
+       "  * y            (y) float64 -2.304e+06 -2.304e+06 ... -2.306e+06 -2.306e+06\n",
+       "  * x            (x) float64 -1.598e+06 -1.598e+06 ... -1.596e+06 -1.596e+06\n",
+       "    spatial_ref  int32 3577\n",
+       "Data variables:\n",
+       "    ndwi         (time, y, x) float32 nan nan nan ... -0.1535 -0.1542 -0.1555\n",
+       "Attributes:\n",
+       "    crs:           EPSG:3577\n",
+       "    grid_mapping:  spatial_ref
" ], "text/plain": [ - "" + "\n", + "Dimensions: (time: 325, y: 237, x: 215)\n", + "Coordinates:\n", + " * time (time) datetime64[ns] 2020-01-01T01:55:31.756059 ... 2022-12...\n", + " * y (y) float64 -2.304e+06 -2.304e+06 ... -2.306e+06 -2.306e+06\n", + " * x (x) float64 -1.598e+06 -1.598e+06 ... -1.596e+06 -1.596e+06\n", + " spatial_ref int32 3577\n", + "Data variables:\n", + " ndwi (time, y, x) float32 nan nan nan ... -0.1535 -0.1542 -0.1555\n", + "Attributes:\n", + " crs: EPSG:3577\n", + " grid_mapping: spatial_ref" ] }, + "execution_count": 6, "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Dimensions: (time: 272, y: 278, x: 289)\n", - "Coordinates:\n", - " * time (time) datetime64[ns] 2020-01-02T00:57:27.736492 ... 2022-12...\n", - " * y (y) float64 -1.286e+06 -1.286e+06 ... -1.288e+06 -1.288e+06\n", - " * x (x) float64 -1.408e+04 -1.408e+04 ... -1.122e+04 -1.12e+04\n", - " spatial_ref int32 3577\n", - "Data variables:\n", - " ndwi (time, y, x) float32 dask.array\n", - "Attributes:\n", - " crs: EPSG:3577\n", - " grid_mapping: spatial_ref\n" - ] + "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ - "/env/lib/python3.10/site-packages/datacube/drivers/driver_cache.py:54: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html\n", - " from pkg_resources import iter_entry_points\n", - "/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.\n", - " _reproject(\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 9.18 s, sys: 720 ms, total: 9.9 s\n", - "Wall time: 52.2 s\n" + "=======\n" ] }, { @@ -962,6 +1750,13 @@ "execution_count": 6, "metadata": {}, "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ">>>>>>> remote\n" + ] } ], "source": [ @@ -973,7 +1768,7 @@ "# Create local dask cluster to improve data load time\n", "client = create_local_dask_cluster(return_client=True)\n", "\n", - "satellite_ds, _, _ = load_data(\n", + "satellite_ds = load_data(\n", " dc=dc,\n", " study_area=study_area,\n", " geom=geom,\n", @@ -995,7 +1790,7 @@ }, { "cell_type": "markdown", - "id": "ad0ab3e3", + "id": "1e110fbd", "metadata": {}, "source": [ "### Load optional topobathy mask\n", @@ -1006,7 +1801,7 @@ { "cell_type": "code", "execution_count": 7, - "id": "2adcc77b", + "id": "86e0342f", "metadata": { "tags": [] }, @@ -1018,7 +1813,7 @@ }, { "cell_type": "markdown", - "id": "fd2c1630", + "id": "61de9aa9", "metadata": {}, "source": [ "### Intertidal elevation\n", @@ -1028,7 +1823,7 @@ { "cell_type": "code", "execution_count": 8, - "id": "45de9453-5ca0-4879-85fd-09864229bbb2", + "id": "498a6da5", "metadata": { "tags": [] }, @@ -1037,7 +1832,13 @@ "name": "stderr", "output_type": "stream", "text": [ - "2024-03-08 03:49:42 INFO Processing: Modelling tide heights for each pixel\n" + "<<<<<<< local \n", + "2024-02-01 04:19:28 INFO Modelling tide heights for each pixel\n", + "=======\n", + "2024-03-08 03:49:42 INFO Processing: Modelling tide heights for each pixel\n", + ">>>>>>> remote \n", + "<<<<<<< local\n", + "=======\n" ] }, { @@ -1053,7 +1854,45 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 35/35 [00:19<00:00, 1.83it/s]\n" + ">>>>>>> remote\n", + "<<<<<<< local \n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating reduced resolution 5000 x 5000 metre tide modelling array\n", + "Modelling tides using FES2014, FES2012, TPXO8-atlas-v1, TPXO9-atlas-v5, EOT20, HAMTIDE11, GOT4.10 in parallel\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n", + ">>>>>>> remote \n", + "<<<<<<< local \n", + "100%|██████████| 35/35 [00:19<00:00, 1.77it/s]\n", + "=======\n", + "100%|██████████| 35/35 [00:19<00:00, 1.83it/s]\n", + ">>>>>>> remote \n", + "<<<<<<< local\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Returning low resolution tide array\n", + "Generating ensemble tide model from point inputs\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n" ] }, { @@ -1079,10 +1918,33 @@ "name": "stderr", "output_type": "stream", "text": [ + ">>>>>>> remote\n", + "<<<<<<< local \n", + "2024-02-01 04:19:52 INFO Masking nodata and adding tide heights to satellite data array\n", + "2024-02-01 04:19:53 INFO Flattening satellite data array and filtering to intertidal candidate pixels\n", + "2024-02-01 04:19:53 INFO Applying valid data mask to constrain study area\n", + "2024-02-01 04:19:53 INFO Running per-pixel rolling median\n", + "=======\n", "2024-03-08 03:50:04 INFO Processing: Masking nodata and adding tide heights to satellite data array\n", "2024-03-08 03:50:04 INFO Processing: Flattening satellite data array and filtering to intertidal candidate pixels\n", "2024-03-08 03:50:04 INFO Processing: Applying valid data mask to constrain study area\n", - "2024-03-08 03:50:05 INFO Processing: Running per-pixel rolling median\n" + "2024-03-08 03:50:05 INFO Processing: Running per-pixel rolling median\n", + ">>>>>>> remote \n", + "<<<<<<< local \n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reducing analysed pixels from 50955 to 8141 (15.98%)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "=======\n" ] }, { @@ -1092,6 +1954,37 @@ "Reducing analysed pixels from 80342 to 42943 (53.45%)\n" ] }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ">>>>>>> remote \n", + "<<<<<<< local\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "9e8537d7507f4ff999103adc0dde4ffa", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/105 [00:00>>>>>> remote\n", + "<<<<<<< local \n", + "=======\n", + "2024-03-08 03:50:31 INFO Processing: Modelling intertidal elevation\n", + ">>>>>>> remote \n" ] }, { @@ -1125,18 +2022,24 @@ "name": "stderr", "output_type": "stream", "text": [ + "<<<<<<< local \n", + "2024-02-01 04:20:06 INFO Unflattening data back to its original spatial dimensions\n", + "2024-02-01 04:20:06 INFO Cleaning inaccurate upper intertidal pixels\n", + "2024-02-01 04:20:06 INFO Successfully completed intertidal elevation modelling\n", + "=======\n", "2024-03-08 03:50:32 INFO Processing: Modelling intertidal uncertainty\n", "2024-03-08 03:50:33 INFO Processing: Unflattening data back to its original spatial dimensions\n", "2024-03-08 03:50:33 INFO Processing: Cleaning inaccurate upper intertidal pixels\n", - "2024-03-08 03:50:33 INFO Processing: Successfully completed intertidal elevation modelling\n" + "2024-03-08 03:50:33 INFO Processing: Successfully completed intertidal elevation modelling\n", + ">>>>>>> remote \n" ] } ], "source": [ "# Model elevation for each pixel\n", - "ds, tide_m = elevation(\n", + "ds, ds_aux, tide_m = elevation(\n", " satellite_ds,\n", - " valid_mask=topobathy_mask,\n", + " valid_mask=topobathy_ds.height_depth > -20,\n", " tide_model=tide_model,\n", " tide_model_dir=tide_model_dir,\n", ")" @@ -1144,7 +2047,7 @@ }, { "cell_type": "markdown", - "id": "82962441", + "id": "704fc753", "metadata": { "tags": [] }, @@ -1155,7 +2058,7 @@ { "cell_type": "code", "execution_count": 9, - "id": "f49adb23-1eb4-457b-bffe-c8b90e117c97", + "id": "f7f72117", "metadata": { "tags": [] }, @@ -1176,13 +2079,13 @@ { "cell_type": "code", "execution_count": null, - "id": "08dbad9b-ad1e-46bd-a55c-dd442e6e8446", + "id": "2e92bfa6", "metadata": { "tags": [] }, "outputs": [], "source": [ - "# ## Plot\n", + "# ## Plot with labels\n", "# labels=['0 Dry',\n", "# '1 Intermittent\\n wet\\n inland',\n", "# '2 Wet inland',\n", @@ -1200,19 +2103,85 @@ }, { "cell_type": "markdown", - "id": "b30492d7", + "id": "063c72b9", "metadata": { "tags": [] }, "source": [ "### Intertidal exposure\n", - "Calculate exposure using the script function" + "Calculate exposure using the script function.\n", + "To calculate exposure for the full time period, leave `filters` commented out or set as ['unfiltered'].\n", + "See the function documentation for the full range of available filters and filter_combinations.\n", + "The code accepts lists of multiple filters and filter_combination tuples." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "e1ff8db8", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Calculating unfiltered exposure\n", + "Creating reduced resolution 5000 x 5000 metre tide modelling array\n", + "Modelling tides using FES2014, FES2012, TPXO8-atlas-v1, TPXO9-atlas-v5, EOT20, HAMTIDE11, GOT4.10 in parallel\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 35/35 [00:19<00:00, 1.78it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Computing tide quantiles\n", + "Returning low resolution tide array\n", + "Generating ensemble tide model from point inputs\n" + ] + } + ], + "source": [ + "%reload_ext autoreload\n", + "\n", + "exposure_filters, tide_cq_dict = exposure(\n", + " dem=ds.elevation,\n", + " start_date=start_date,\n", + " end_date=end_date,\n", + " modelled_freq = modelled_freq,\n", + " tide_model=tide_model,\n", + " tide_model_dir=tide_model_dir,\n", + " # filters=filters,\n", + " # filters_combined=filters_combined,\n", + ")\n", + "\n", + "for x in list(exposure_filters.keys()):\n", + " ds['exposure_'+str(x)]=exposure_filters[str(x)]" + ] + }, + { + "cell_type": "markdown", + "id": "a3157252", + "metadata": {}, + "source": [ + "### Spread and offset\n", + "Calculate the spread and high/low tide biases of input observed images as a percentage of the modelled tide heights.\n", + "\n", + "Warning: this code will only work if the exposure function has been run and produced 'unfiltered' exposure results.\n" ] }, { "cell_type": "code", "execution_count": 10, - "id": "ad27ef63-3e62-4167-8369-0d6fa79bab8e", + "id": "e8c240d4", "metadata": { "tags": [] }, @@ -1270,19 +2239,10 @@ ")" ] }, - { - "cell_type": "markdown", - "id": "89b86fb9", - "metadata": {}, - "source": [ - "### Spread and offset\n", - "Calculate the spread and high/low tide biases of input observed images as a percentage of the modelled tide heights" - ] - }, { "cell_type": "code", - "execution_count": 11, - "id": "ccc3123d", + "execution_count": null, + "id": "f7a4cd0a", "metadata": { "tags": [] }, @@ -1297,13 +2257,13 @@ " ds[\"ta_offset_low\"],\n", " ds[\"ta_offset_high\"],\n", ") = bias_offset(\n", - " tide_m=tide_m, tide_cq=tide_cq, extents=ds.extents, lot_hot=True, lat_hat=True\n", + " tide_m=tide_m, tide_cq=tide_cq_dict['unfiltered'], extents=ds.extents, lot_hot=True, lat_hat=True\n", ")" ] }, { "cell_type": "markdown", - "id": "2986d91d", + "id": "01e2aa71", "metadata": {}, "source": [ "## Plot all layers" @@ -1311,8 +2271,8 @@ }, { "cell_type": "code", - "execution_count": 12, - "id": "b80d35e1", + "execution_count": null, + "id": "187c867c", "metadata": { "tags": [] }, @@ -1895,8 +2855,8 @@ }, { "cell_type": "code", - "execution_count": 13, - "id": "0b8f96f8", + "execution_count": null, + "id": "f188da27", "metadata": { "tags": [] }, @@ -1913,14 +2873,15 @@ } ], "source": [ - "fig = plt.figure(figsize=(12, 12), tight_layout=True)\n", + "fig = plt.figure(figsize=(16, 18), tight_layout=True)\n", "ax_dict = fig.subplot_mosaic(\n", " \"\"\"\n", " AAAABBBBCCCC\n", " DDDEEEFFFGGG\n", " HHHIIIJJJKKK\n", + " LLLMMM......\n", " \"\"\",\n", - " height_ratios=[1, 0.8, 0.8],\n", + " height_ratios=[1, 0.8, 0.8, 0.8],\n", ")\n", "\n", "# label_params = dict(add_labels=False, yticks=[], xticks=[])\n", @@ -1947,7 +2908,7 @@ "ax_dict[\"B\"].set_facecolor(\"#2E2E2E\")\n", "\n", "# Plot Exposure\n", - "ds[\"exposure\"].plot.imshow(\n", + "ds[\"exposure_unfiltered\"].plot.imshow(\n", " ax=ax_dict[\"C\"],\n", " cmap=\"RdYlGn\",\n", " levels=np.arange(0, 100, 1),\n", @@ -1961,7 +2922,7 @@ "ax_dict[\"D\"].set_title(\"Extents\")\n", "\n", "# Plot the observation spread\n", - "ds[\"ta_spread\"].plot.imshow(\n", + "ds[\"oa_spread\"].plot.imshow(\n", " ax=ax_dict[\"E\"],\n", " vmin=0,\n", " vmax=100,\n", @@ -1971,7 +2932,7 @@ "ax_dict[\"E\"].set_title(\"Observation Spread (%)\")\n", "\n", "# Plot the high-tide offset\n", - "ds[\"ta_offset_high\"].plot.imshow(\n", + "ds[\"oa_offset_hightide\"].plot.imshow(\n", " ax=ax_dict[\"F\"],\n", " vmin=0,\n", " vmax=40,\n", @@ -1981,7 +2942,7 @@ "ax_dict[\"F\"].set_title(\"High-tide offset (%)\")\n", "\n", "# Plot the low-tide offset\n", - "ds[\"ta_offset_low\"].plot.imshow(\n", + "ds[\"oa_offset_lowtide\"].plot.imshow(\n", " ax=ax_dict[\"G\"],\n", " vmin=0,\n", " vmax=40,\n", @@ -1991,38 +2952,68 @@ "ax_dict[\"G\"].set_title(\"Low-tide offset (%)\")\n", "\n", "# Plot the LAT\n", - "ds[\"ta_lat\"].plot.imshow(\n", + "ds[\"oa_lat\"].plot.imshow(\n", " ax=ax_dict[\"H\"],\n", - " vmin=ds[\"ta_lat\"].min(),\n", - " vmax=ds[\"ta_hat\"].max(),\n", + " vmin=ds[\"oa_lat\"].min(),\n", + " vmax=ds[\"oa_hat\"].max(),\n", " add_labels=False,\n", ")\n", "ax_dict[\"H\"].set_title(\"Lowest Astronomical Tide\")\n", "\n", "# Plot the LOT\n", - "ds[\"ta_lot\"].plot.imshow(\n", - " ax=ax_dict[\"I\"], vmin=ds[\"ta_lat\"].min(), vmax=ds[\"ta_hat\"].max()\n", + "ds[\"oa_lot\"].plot.imshow(\n", + " ax=ax_dict[\"I\"], vmin=ds[\"oa_lat\"].min(), vmax=ds[\"oa_hat\"].max()\n", ")\n", "ax_dict[\"I\"].set_title(\"Lowest Observed Tide\")\n", "\n", "# Plot the HAT\n", - "ds[\"ta_hat\"].plot.imshow(\n", + "ds[\"oa_hat\"].plot.imshow(\n", " ax=ax_dict[\"J\"],\n", - " vmin=ds[\"ta_lat\"].min(),\n", - " vmax=ds[\"ta_hat\"].max(),\n", + " vmin=ds[\"oa_lat\"].min(),\n", + " vmax=ds[\"oa_hat\"].max(),\n", " add_labels=False,\n", " yticks=[],\n", ")\n", "ax_dict[\"J\"].set_title(\"Highest Astronomical Tide\")\n", "\n", "# Plot the HOT\n", - "ds[\"ta_hot\"].plot.imshow(\n", + "ds[\"oa_hot\"].plot.imshow(\n", " ax=ax_dict[\"K\"],\n", - " vmin=ds[\"ta_lat\"].min(),\n", - " vmax=ds[\"ta_hat\"].max(),\n", + " vmin=ds[\"oa_lat\"].min(),\n", + " vmax=ds[\"oa_hat\"].max(),\n", ")\n", "ax_dict[\"K\"].set_title(\"Highest Observed Tide\")\n", "\n", + "# Plot the high and low tidelines with respective offset\n", + "ax_dict[\"L\"].set_title(\"Lowtide line and lowtide offset\")\n", + "# lowtideline.plot(\n", + "# column=\"offset_lowtide\",\n", + "# legend=True,\n", + "# vmin=0,\n", + "# vmax=40,\n", + "# cmap=\"magma\",\n", + "# ax=ax_dict[\"L\"],\n", + "# zorder=2,\n", + "# )\n", + "# tidelines_gdf.loc[[0], \"geometry\"].plot(ax=ax_dict[\"L\"], zorder=1)\n", + "ax_dict[\"L\"].set_xlim(left=ds.elevation.x.min(), right=ds.elevation.x.max())\n", + "ax_dict[\"L\"].set_ylim(bottom=ds.elevation.y.min(), top=ds.elevation.y.max())\n", + "\n", + "ax_dict[\"M\"].set_title(\"Hightide line and hightide offset\")\n", + "# hightideline.plot(\n", + "# column=\"offset_hightide\",\n", + "# legend=True,\n", + "# vmin=0,\n", + "# vmax=40,\n", + "# cmap=\"magma\",\n", + "# ax=ax_dict[\"M\"],\n", + "# zorder=2,\n", + "# )\n", + "# tidelines_gdf.loc[[1], \"geometry\"].plot(ax=ax_dict[\"M\"], zorder=1)\n", + "ax_dict[\"M\"].set_yticks([])\n", + "ax_dict[\"M\"].set_xlim(left=ds.elevation.x.min(), right=ds.elevation.x.max())\n", + "ax_dict[\"M\"].set_ylim(bottom=ds.elevation.y.min(), top=ds.elevation.y.max())\n", + "\n", "# Remove axis labels\n", "for label, ax in ax_dict.items():\n", " ax.set_yticks([])\n", @@ -2033,7 +3024,7 @@ }, { "cell_type": "markdown", - "id": "e91f81d3", + "id": "6f63fbdf", "metadata": { "tags": [] }, @@ -2044,7 +3035,7 @@ { "cell_type": "code", "execution_count": 14, - "id": "1a027ec4", + "id": "69e246d3", "metadata": {}, "outputs": [], "source": [ @@ -2063,7 +3054,7 @@ { "cell_type": "code", "execution_count": 15, - "id": "4cc049a7", + "id": "c3399861", "metadata": {}, "outputs": [], "source": [ @@ -2073,7 +3064,7 @@ }, { "cell_type": "markdown", - "id": "f84ff11f", + "id": "ecbb2c87", "metadata": { "tags": [] }, @@ -2084,7 +3075,7 @@ { "cell_type": "code", "execution_count": 17, - "id": "a59ba636", + "id": "0c6c7149", "metadata": {}, "outputs": [], "source": [