analysis/webservice/algorithms/LongitudeLatitudeMap.py

# 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 math
from datetime import datetime

from pytz import timezone
from shapely.geometry import box

from webservice.NexusHandler import nexus_handler
from webservice.algorithms.NexusCalcHandler import NexusCalcHandler
from webservice.webmodel import NexusResults, NexusProcessingException

SENTINEL = 'STOP'
EPOCH = timezone('UTC').localize(datetime(1970, 1, 1))
tile_service = None


@nexus_handler
class LongitudeLatitudeMapCalcHandlerImpl(NexusCalcHandler):
    name = "Longitude/Latitude Time Average Map"
    path = "/longitudeLatitudeMap"
    description = "Computes a Latitude/Longitude Time Average plot given an arbitrary geographical area and time range"
    params = {
        "ds": {
            "name": "Dataset",
            "type": "string",
            "description": "One or more comma-separated dataset shortnames"
        },
        "minLat": {
            "name": "Minimum Latitude",
            "type": "float",
            "description": "Minimum (Southern) bounding box Latitude"
        },
        "maxLat": {
            "name": "Maximum Latitude",
            "type": "float",
            "description": "Maximum (Northern) bounding box Latitude"
        },
        "minLon": {
            "name": "Minimum Longitude",
            "type": "float",
            "description": "Minimum (Western) bounding box Longitude"
        },
        "maxLon": {
            "name": "Maximum Longitude",
            "type": "float",
            "description": "Maximum (Eastern) bounding box Longitude"
        },
        "startTime": {
            "name": "Start Time",
            "type": "string",
            "description": "Starting time in format YYYY-MM-DDTHH:mm:ssZ or seconds since epoch (Jan 1st, 1970)"
        },
        "endTime": {
            "name": "End Time",
            "type": "string",
            "description": "Ending time in format YYYY-MM-DDTHH:mm:ssZ or seconds since epoch (Jan 1st, 1970)"
        }
    }
    singleton = True

    def parse_arguments(self, request):
        # Parse input arguments
        try:
            ds = request.get_dataset()[0]
        except:
            raise NexusProcessingException(reason="'ds' argument is required", code=400)

        try:
            bounding_polygon = box(request.get_min_lon(), request.get_min_lat(), request.get_max_lon(),
                                   request.get_max_lat())
        except:
            raise NexusProcessingException(
                reason="'minLon', 'minLat', 'maxLon', and 'maxLat' arguments are required.",
                code=400)

        try:
            start_time = request.get_start_datetime()
        except:
            raise NexusProcessingException(
                reason="'startTime' argument is required. Can be int value milliseconds from epoch or string format YYYY-MM-DDTHH:mm:ssZ",
                code=400)
        try:
            end_time = request.get_end_datetime()
        except:
            raise NexusProcessingException(
                reason="'endTime' argument is required. Can be int value milliseconds from epoch or string format YYYY-MM-DDTHH:mm:ssZ",
                code=400)

        start_seconds_from_epoch = int((start_time - EPOCH).total_seconds())
        end_seconds_from_epoch = int((end_time - EPOCH).total_seconds())

        return ds, bounding_polygon, start_seconds_from_epoch, end_seconds_from_epoch

    def calc(self, request, **args):

        ds, bounding_polygon, start_seconds_from_epoch, end_seconds_from_epoch = self.parse_arguments(request)

        boxes = self._get_tile_service().get_distinct_bounding_boxes_in_polygon(bounding_polygon, ds,
                                                                          start_seconds_from_epoch,
                                                                          end_seconds_from_epoch)
        point_avg_over_time = lat_lon_map_driver(bounding_polygon, start_seconds_from_epoch, end_seconds_from_epoch, ds,
                                                 [a_box.bounds for a_box in boxes])

        kwargs = {
            "minLon": bounding_polygon.bounds[0],
            "minLat": bounding_polygon.bounds[1],
            "maxLon": bounding_polygon.bounds[2],
            "maxLat": bounding_polygon.bounds[3],
            "ds": ds,
            "startTime": datetime.utcfromtimestamp(start_seconds_from_epoch).strftime('%Y-%m-%dT%H:%M:%SZ'),
            "endTime": datetime.utcfromtimestamp(end_seconds_from_epoch).strftime('%Y-%m-%dT%H:%M:%SZ')
        }
        return LongitudeLatitudeMapResults(
            results=LongitudeLatitudeMapCalcHandlerImpl.results_to_dicts(point_avg_over_time), meta=None,
            **kwargs)

    @staticmethod
    def results_to_dicts(results):

        # ((lon, lat), (slope, intercept, r_value, p_value, std_err, mean, pmax, pmin, pstd, pcnt))
        return [{
            'lon': result[0][0],
            'lat': result[0][1],
            'slope': result[1][0] if not math.isnan(result[1][0]) else 'NaN',
            'intercept': result[1][1] if not math.isnan(result[1][1]) else 'NaN',
            'r': result[1][2],
            'p': result[1][3],
            'stderr': result[1][4] if not math.isinf(result[1][4]) else 'Inf',
            'avg': result[1][5],
            'max': result[1][6],
            'min': result[1][7],
            'std': result[1][8],
            'cnt': result[1][9],
        } for result in results]


def pool_initializer():
    from nexustiles.nexustiles import NexusTileService
    global tile_service
    tile_service = NexusTileService()
    # TODO This is a hack to make sure each sub-process uses it's own connection to cassandra. data-access needs to be updated
    from cassandra.cqlengine import connection
    from multiprocessing import current_process

    connection.register_connection(current_process().name, [host.address for host in connection.get_session().hosts])
    connection.set_default_connection(current_process().name)


def lat_lon_map_driver(search_bounding_polygon, search_start, search_end, ds, distinct_boxes):
    from functools import partial
    from nexustiles.nexustiles import NexusTileService
    # Start new processes to handle the work
    # pool = Pool(5, pool_initializer)

    func = partial(regression_on_tiles, search_bounding_polygon_wkt=search_bounding_polygon.wkt,
                   search_start=search_start, search_end=search_end, ds=ds)

    global tile_service
    tile_service = NexusTileService()
    map_result = list(map(func, distinct_boxes))
    return [item for sublist in map_result for item in sublist]
    # TODO Use for multiprocessing:
    # result = pool.map_async(func, distinct_boxes)
    #
    # pool.close()
    # while not result.ready():
    #     print "Not ready"
    #     time.sleep(5)
    #
    # print result.successful()
    # print result.get()[0]
    # pool.join()
    # pool.terminate()
    #
    # return [item for sublist in result.get() for item in sublist]


def calc_linear_regression(arry, xarry):
    from scipy.stats import linregress
    slope, intercept, r_value, p_value, std_err = linregress(arry, xarry)
    return slope, intercept, r_value, p_value, std_err


def regression_on_tiles(tile_bounds, search_bounding_polygon_wkt, search_start, search_end, ds):
    if len(tile_bounds) < 1:
        return []

    import numpy as np
    import operator
    from shapely import wkt
    from shapely.geometry import box
    search_bounding_shape = wkt.loads(search_bounding_polygon_wkt)
    tile_bounding_shape = box(*tile_bounds)

    # Load all tiles for given (exact) bounding box across the search time range
    tiles = tile_service.find_tiles_by_exact_bounds(tile_bounds, ds, search_start, search_end)
    if search_bounding_shape.contains(tile_bounding_shape):
        # The tile bounds are totally contained in the search area, we don't need to mask it.
        pass
    else:
        # The tile bounds cross the search area borders, we need to mask the tiles to the search area
        tiles = tile_service.mask_tiles_to_polygon(wkt.loads(search_bounding_polygon_wkt), tiles)
    # If all tiles end up being masked, there is no work to do
    if len(tiles) < 1:
        return []
    tiles.sort(key=operator.attrgetter('min_time'))

    stacked_tile_data = np.stack(tuple([np.squeeze(tile.data, 0) for tile in tiles]))
    stacked_tile_lons = np.stack([tile.longitudes for tile in tiles])
    stacked_tile_lats = np.stack([tile.latitudes for tile in tiles])

    x_array = np.arange(stacked_tile_data.shape[0])
    point_regressions = np.apply_along_axis(calc_linear_regression, 0, stacked_tile_data, x_array)
    point_means = np.nanmean(stacked_tile_data, 0)
    point_maximums = np.nanmax(stacked_tile_data, 0)
    point_minimums = np.nanmin(stacked_tile_data, 0)
    point_st_deviations = np.nanstd(stacked_tile_data, 0)
    point_counts = np.ma.count(np.ma.masked_invalid(stacked_tile_data), 0)

    results = []
    for lat_idx, lon_idx in np.ndindex(point_means.shape):
        lon, lat = np.max(stacked_tile_lons[:, lon_idx]), np.max(stacked_tile_lats[:, lat_idx])
        pcnt = point_counts[lat_idx, lon_idx]

        if pcnt == 0:
            continue

        mean = point_means[lat_idx, lon_idx]
        pmax = point_maximums[lat_idx, lon_idx]
        pmin = point_minimums[lat_idx, lon_idx]
        pstd = point_st_deviations[lat_idx, lon_idx]
        slope, intercept, r_value, p_value, std_err = point_regressions[:, lat_idx, lon_idx]

        results.append(((lon, lat), (slope, intercept, r_value, p_value, std_err, mean, pmax, pmin, pstd, pcnt)))

    return results


class LongitudeLatitudeMapResults(NexusResults):
    def __init__(self, results=None, meta=None, computeOptions=None, **kwargs):
        NexusResults.__init__(self, results=results, meta=meta, stats=None, computeOptions=computeOptions, **kwargs)