# Licensed to the Apache Software Foundation (ASF) under one or more # contributor license agreements. See the NOTICE file distributed with # this work for additional information regarding copyright ownership. # The ASF licenses this file to You under the Apache License, Version 2.0 # (the "License"); you may not use this file except in compliance with # the License. You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. import logging import traceback from io import StringIO from datetime import datetime from multiprocessing.dummy import Pool, Manager import matplotlib.dates as mdates import matplotlib.pyplot as plt import numpy as np from nexustiles.nexustiles import NexusTileService from scipy import stats from webservice import Filtering as filt from webservice.NexusHandler import nexus_handler, DEFAULT_PARAMETERS_SPEC from webservice.algorithms.NexusCalcHandler import NexusCalcHandler from webservice.webmodel import NexusResults, NexusProcessingException, NoDataException SENTINEL = 'STOP' logger = logging.getLogger(__name__) @nexus_handler class TimeSeriesCalcHandlerImpl(NexusCalcHandler): name = "Time Series Solr" path = "/statsSolr" description = "Computes a time series plot between one or more datasets given an arbitrary geographical area and time range" params = DEFAULT_PARAMETERS_SPEC singleton = True def calc(self, computeOptions, **args): """ :param computeOptions: StatsComputeOptions :param args: dict :return: """ ds = computeOptions.get_dataset() if type(ds) != list and type(ds) != tuple: ds = (ds,) resultsRaw = [] for shortName in ds: results, meta = self.getTimeSeriesStatsForBoxSingleDataSet(computeOptions.get_min_lat(), computeOptions.get_max_lat(), computeOptions.get_min_lon(), computeOptions.get_max_lon(), shortName, computeOptions.get_start_time(), computeOptions.get_end_time(), computeOptions.get_apply_seasonal_cycle_filter(), computeOptions.get_apply_low_pass_filter()) resultsRaw.append([results, meta]) results = self._mergeResults(resultsRaw) if len(ds) == 2: stats = self.calculateComparisonStats(results, suffix="") if computeOptions.get_apply_seasonal_cycle_filter(): s = self.calculateComparisonStats(results, suffix="Seasonal") stats = self._mergeDicts(stats, s) if computeOptions.get_apply_low_pass_filter(): s = self.calculateComparisonStats(results, suffix="LowPass") stats = self._mergeDicts(stats, s) if computeOptions.get_apply_seasonal_cycle_filter() and computeOptions.get_apply_low_pass_filter(): s = self.calculateComparisonStats(results, suffix="SeasonalLowPass") stats = self._mergeDicts(stats, s) else: stats = {} meta = [] for singleRes in resultsRaw: meta.append(singleRes[1]) res = TimeSeriesResults(results=results, meta=meta, stats=stats, computeOptions=computeOptions) return res def getTimeSeriesStatsForBoxSingleDataSet(self, min_lat, max_lat, min_lon, max_lon, ds, start_time=0, end_time=-1, applySeasonalFilter=True, applyLowPass=True): daysinrange = self._get_tile_service().find_days_in_range_asc(min_lat, max_lat, min_lon, max_lon, ds, start_time, end_time) if len(daysinrange) == 0: raise NoDataException(reason="No data found for selected timeframe") maxprocesses = int(self.algorithm_config.get("multiprocessing", "maxprocesses")) results = [] if maxprocesses == 1: calculator = TimeSeriesCalculator() for dayinseconds in daysinrange: result = calculator.calc_average_on_day(min_lat, max_lat, min_lon, max_lon, ds, dayinseconds) results.append(result) else: # Create a task to calc average difference for each day manager = Manager() work_queue = manager.Queue() done_queue = manager.Queue() for dayinseconds in daysinrange: work_queue.put( ('calc_average_on_day', min_lat, max_lat, min_lon, max_lon, ds, dayinseconds)) [work_queue.put(SENTINEL) for _ in range(0, maxprocesses)] # Start new processes to handle the work pool = Pool(maxprocesses) [pool.apply_async(pool_worker, (work_queue, done_queue)) for _ in range(0, maxprocesses)] pool.close() # Collect the results as [(day (in ms), average difference for that day)] for i in range(0, len(daysinrange)): result = done_queue.get() try: error_str = result['error'] logger.error(error_str) raise NexusProcessingException(reason="Error calculating average by day.") except KeyError: pass results.append(result) pool.terminate() manager.shutdown() results = sorted(results, key=lambda entry: entry["time"]) filt.applyAllFiltersOnField(results, 'mean', applySeasonal=applySeasonalFilter, applyLowPass=applyLowPass) filt.applyAllFiltersOnField(results, 'max', applySeasonal=applySeasonalFilter, applyLowPass=applyLowPass) filt.applyAllFiltersOnField(results, 'min', applySeasonal=applySeasonalFilter, applyLowPass=applyLowPass) return results, {} def calculateComparisonStats(self, results, suffix=""): xy = [[], []] for item in results: if len(item) == 2: xy[item[0]["ds"]].append(item[0]["mean%s" % suffix]) xy[item[1]["ds"]].append(item[1]["mean%s" % suffix]) slope, intercept, r_value, p_value, std_err = stats.linregress(xy[0], xy[1]) comparisonStats = { "slope%s" % suffix: slope, "intercept%s" % suffix: intercept, "r%s" % suffix: r_value, "p%s" % suffix: p_value, "err%s" % suffix: std_err } return comparisonStats class TimeSeriesResults(NexusResults): LINE_PLOT = "line" SCATTER_PLOT = "scatter" __SERIES_COLORS = ['red', 'blue'] def __init__(self, results=None, meta=None, stats=None, computeOptions=None): NexusResults.__init__(self, results=results, meta=meta, stats=stats, computeOptions=computeOptions) def toImage(self): type = self.computeOptions().get_plot_type() if type == TimeSeriesResults.LINE_PLOT or type == "default": return self.createLinePlot() elif type == TimeSeriesResults.SCATTER_PLOT: return self.createScatterPlot() else: raise Exception("Invalid or unsupported time series plot specified") def createScatterPlot(self): timeSeries = [] series0 = [] series1 = [] res = self.results() meta = self.meta() plotSeries = self.computeOptions().get_plot_series() if self.computeOptions is not None else None if plotSeries is None: plotSeries = "mean" for m in res: if len(m) == 2: timeSeries.append(datetime.fromtimestamp(m[0]["time"] / 1000)) series0.append(m[0][plotSeries]) series1.append(m[1][plotSeries]) title = ', '.join(set([m['title'] for m in meta])) sources = ', '.join(set([m['source'] for m in meta])) dateRange = "%s - %s" % (timeSeries[0].strftime('%b %Y'), timeSeries[-1].strftime('%b %Y')) fig, ax = plt.subplots() fig.set_size_inches(11.0, 8.5) ax.scatter(series0, series1, alpha=0.5) ax.set_xlabel(meta[0]['units']) ax.set_ylabel(meta[1]['units']) ax.set_title("%s\n%s\n%s" % (title, sources, dateRange)) par = np.polyfit(series0, series1, 1, full=True) slope = par[0][0] intercept = par[0][1] xl = [min(series0), max(series0)] yl = [slope * xx + intercept for xx in xl] plt.plot(xl, yl, '-r') # r = self.stats()["r"] # plt.text(0.5, 0.5, "r = foo") ax.grid(True) fig.tight_layout() sio = StringIO() plt.savefig(sio, format='png') return sio.getvalue() def createLinePlot(self): nseries = len(self.meta()) res = self.results() meta = self.meta() timeSeries = [datetime.fromtimestamp(m[0]["time"] / 1000) for m in res] means = [[np.nan] * len(res) for n in range(0, nseries)] plotSeries = self.computeOptions().get_plot_series() if self.computeOptions is not None else None if plotSeries is None: plotSeries = "mean" for n in range(0, len(res)): timeSlot = res[n] for seriesValues in timeSlot: means[seriesValues['ds']][n] = seriesValues[plotSeries] x = timeSeries fig, axMain = plt.subplots() fig.set_size_inches(11.0, 8.5) fig.autofmt_xdate() title = ', '.join(set([m['title'] for m in meta])) sources = ', '.join(set([m['source'] for m in meta])) dateRange = "%s - %s" % (timeSeries[0].strftime('%b %Y'), timeSeries[-1].strftime('%b %Y')) axMain.set_title("%s\n%s\n%s" % (title, sources, dateRange)) axMain.set_xlabel('Date') axMain.grid(True) axMain.xaxis.set_major_locator(mdates.YearLocator()) axMain.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y')) axMain.xaxis.set_minor_locator(mdates.MonthLocator()) axMain.format_xdata = mdates.DateFormatter('%Y-%m-%d') plots = [] for n in range(0, nseries): if n == 0: ax = axMain else: ax = ax.twinx() plots += ax.plot(x, means[n], color=self.__SERIES_COLORS[n], zorder=10, linewidth=3, label=meta[n]['title']) ax.set_ylabel(meta[n]['units']) labs = [l.get_label() for l in plots] axMain.legend(plots, labs, loc=0) sio = StringIO() plt.savefig(sio, format='png') return sio.getvalue() class TimeSeriesCalculator(object): def __init__(self): self.__tile_service = NexusTileService() def calc_average_on_day(self, min_lat, max_lat, min_lon, max_lon, dataset, timeinseconds): # Get stats using solr only ds1_nexus_tiles_stats = self.__tile_service.get_stats_within_box_at_time(min_lat, max_lat, min_lon, max_lon, dataset, timeinseconds) data_min_within = min([tile["tile_min_val_d"] for tile in ds1_nexus_tiles_stats]) data_max_within = max([tile["tile_max_val_d"] for tile in ds1_nexus_tiles_stats]) data_sum_within = sum([tile["product(tile_avg_val_d, tile_count_i)"] for tile in ds1_nexus_tiles_stats]) data_count_within = sum([tile["tile_count_i"] for tile in ds1_nexus_tiles_stats]) # Get boundary tiles and calculate stats ds1_nexus_tiles = self.__tile_service.get_boundary_tiles_at_time(min_lat, max_lat, min_lon, max_lon, dataset, timeinseconds) tile_data_agg = np.ma.array([tile.data for tile in ds1_nexus_tiles]) data_min_boundary = np.ma.min(tile_data_agg) data_max_boundary = np.ma.max(tile_data_agg) # daily_mean = np.ma.mean(tile_data_agg).item() data_sum_boundary = np.ma.sum(tile_data_agg) data_count_boundary = np.ma.count(tile_data_agg).item() # data_std = np.ma.std(tile_data_agg) # Combine stats data_min = min(data_min_within, data_min_boundary) data_max = max(data_max_within, data_max_boundary) data_count = data_count_within + data_count_boundary daily_mean = (data_sum_within + data_sum_boundary) / data_count data_std = 0 # Return Stats by day stat = { 'min': data_min, 'max': data_max, 'mean': daily_mean, 'cnt': data_count, 'std': data_std, 'time': int(timeinseconds) } return stat def pool_worker(work_queue, done_queue): try: calculator = TimeSeriesCalculator() for work in iter(work_queue.get, SENTINEL): scifunction = work[0] args = work[1:] result = calculator.__getattribute__(scifunction)(*args) done_queue.put(result) except Exception as e: e_str = traceback.format_exc(e) done_queue.put({'error': e_str})