analysis/webservice/algorithms/doms/MatchupQuery.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
import uuid
from datetime import datetime

import numpy as np
import utm
from nexustiles.model.nexusmodel import get_approximate_value_for_lat_lon
from scipy import spatial

from . import BaseDomsHandler
from . import ResultsStorage
from . import datafetch
from . import fetchedgeimpl
from . import geo
from . import workerthread
from webservice.NexusHandler import nexus_handler


@nexus_handler
class CombinedDomsMatchupQueryHandler(BaseDomsHandler.BaseDomsQueryCalcHandler):
    name = "Experimental Combined DOMS In-Situ Matchup"
    path = "/domsmatchup"
    description = ""
    params = {}
    singleton = True

    def __init__(self):
        BaseDomsHandler.BaseDomsQueryCalcHandler.__init__(self)

    def fetchData(self, endpoints, startTime, endTime, bbox, depth_min, depth_max, platforms):

        boundsConstrainer = geo.BoundsConstrainer(asString=bbox)
        threads = []
        for endpoint in endpoints:
            thread = workerthread.WorkerThread(datafetch.fetchData,
                                               params=(endpoint, startTime, endTime, bbox, depth_min, depth_max))
            threads.append(thread)
        workerthread.wait(threads, startFirst=True, poll=0.01)

        data2 = []
        for thread in threads:
            data, bounds = thread.results
            data2 += data
            boundsConstrainer.testOtherConstrainer(bounds)

        return data2, boundsConstrainer

    def __parseDatetime(self, dtString):
        dt = datetime.strptime(dtString, "%Y-%m-%dT%H:%M:%SZ")
        epoch = datetime.utcfromtimestamp(0)
        time = (dt - epoch).total_seconds() * 1000.0
        return time

    def calc(self, computeOptions, **args):
        primary = computeOptions.get_argument("primary", None)
        matchup = computeOptions.get_argument("matchup", None)
        startTime = computeOptions.get_argument("s", None)
        endTime = computeOptions.get_argument("e", None)
        bbox = computeOptions.get_argument("b", None)
        timeTolerance = computeOptions.get_float_arg("tt")
        depth_min = computeOptions.get_float_arg("depthMin", default=None)
        depth_max = computeOptions.get_float_arg("depthMax", default=None)
        radiusTolerance = computeOptions.get_float_arg("rt")
        platforms = computeOptions.get_argument("platforms", None)

        if primary is None or len(primary) == 0:
            raise Exception("No primary dataset specified")

        if matchup is None or len(matchup) == 0:
            raise Exception("No matchup datasets specified")

        start = self._now()

        primarySpec = self.getDataSourceByName(primary)
        if primarySpec is None:
            raise Exception("Specified primary dataset not found using identifier '%s'" % primary)

        primaryData, bounds = self.fetchData([primarySpec], startTime, endTime, bbox, depth_min, depth_max, platforms)

        primaryContext = MatchupContext(primaryData)

        matchupIds = matchup.split(",")

        for matchupId in matchupIds:
            matchupSpec = self.getDataSourceByName(matchupId)

            if matchupSpec is not None:  # Then it's in the in-situ configuration
                proc = InsituDatasetProcessor(primaryContext, matchupSpec, startTime, endTime, bbox, depth_min,
                                              depth_max,
                                              platforms, timeTolerance, radiusTolerance)
                proc.start()
            else:  # We assume it to be a Nexus tiled dataset

                '''
                Single Threaded at the moment...
                '''
                daysinrange = self._get_tile_service().find_days_in_range_asc(bounds.south, bounds.north, bounds.west,
                                                                        bounds.east, matchupId,
                                                                        self.__parseDatetime(startTime) / 1000,
                                                                        self.__parseDatetime(endTime) / 1000)

                tilesByDay = {}
                for dayTimestamp in daysinrange:
                    ds1_nexus_tiles = self._get_tile_service().get_tiles_bounded_by_box_at_time(bounds.south, bounds.north,
                                                                                          bounds.west, bounds.east,
                                                                                          matchupId, dayTimestamp)

                    # print "***", type(ds1_nexus_tiles)
                    # print ds1_nexus_tiles[0].__dict__
                    tilesByDay[dayTimestamp] = ds1_nexus_tiles

                primaryContext.processGridded(tilesByDay, matchupId, radiusTolerance, timeTolerance)

        matches, numMatches = primaryContext.getFinal(len(matchupIds))

        end = self._now()

        args = {
            "primary": primary,
            "matchup": matchupIds,
            "startTime": startTime,
            "endTime": endTime,
            "bbox": bbox,
            "timeTolerance": timeTolerance,
            "depthMin": depth_min,
            "depthMax": depth_max,
            "radiusTolerance": radiusTolerance,
            "platforms": platforms
        }

        details = {
            "timeToComplete": (end - start),
            "numInSituRecords": primaryContext.insituCount,
            "numInSituMatched": primaryContext.insituMatches,
            "numGriddedChecked": primaryContext.griddedCount,
            "numGriddedMatched": primaryContext.griddedMatched
        }

        with ResultsStorage.ResultsStorage() as resultsStorage:
            execution_id = resultsStorage.insertResults(results=matches, params=args, stats=details, startTime=start,
                                                        completeTime=end, userEmail="")

        return BaseDomsHandler.DomsQueryResults(results=matches, args=args, details=details, bounds=None, count=None,
                                                computeOptions=None, executionId=execution_id)


class MatchupContextMap:
    def __init__(self):
        pass

    def add(self, context):
        pass

    def delete(self, context):
        pass


class MatchupContext:
    def __init__(self, primaryData):
        self.id = str(uuid.uuid4())

        self.griddedCount = 0
        self.griddedMatched = 0

        self.insituCount = len(primaryData)
        self.insituMatches = 0

        self.primary = primaryData
        for r in self.primary:
            r["matches"] = []

        self.data = []
        for s in primaryData:
            u = utm.from_latlon(s["y"], s["x"])
            v = (u[0], u[1], 0.0)
            self.data.append(v)

        if len(self.data) > 0:
            self.tree = spatial.KDTree(self.data)
        else:
            self.tree = None

    def getFinal(self, minMatchesToInclude):

        matched = []
        ttlMatches = 0
        for m in self.primary:
            if len(m["matches"]) >= minMatchesToInclude:
                matched.append(m)
                ttlMatches += len(m["matches"])

        return matched, ttlMatches

    def processGridded(self, tilesByDay, source, xyTolerance, timeTolerance):
        for r in self.primary:
            foundSatNodes = self.__getSatNodeForLatLonAndTime(tilesByDay, source, r["y"], r["x"], r["time"],
                                                              xyTolerance)
            self.griddedCount += 1
            self.griddedMatched += len(foundSatNodes)
            r["matches"].extend(foundSatNodes)

    def processInSitu(self, records, xyTolerance, timeTolerance):
        if self.tree is not None:
            for s in records:
                self.insituCount += 1
                u = utm.from_latlon(s["y"], s["x"])
                coords = np.array([u[0], u[1], 0])
                ball = self.tree.query_ball_point(coords, xyTolerance)

                self.insituMatches += len(ball)

                for i in ball:
                    match = self.primary[i]
                    if abs(match["time"] - s["time"]) <= (timeTolerance * 1000.0):
                        match["matches"].append(s)

    def __getValueForLatLon(self, chunks, lat, lon, arrayName="data"):
        value = get_approximate_value_for_lat_lon(chunks, lat, lon, arrayName)
        return value

    def __checkNumber(self, value):
        if isinstance(value, float) and (math.isnan(value) or value == np.nan):
            value = None
        elif value is not None:
            value = float(value)
        return value

    def __buildSwathIndexes(self, chunk):
        latlons = []
        utms = []
        indexes = []
        for i in range(0, len(chunk.latitudes)):
            _lat = chunk.latitudes[i]
            if isinstance(_lat, np.ma.core.MaskedConstant):
                continue
            for j in range(0, len(chunk.longitudes)):
                _lon = chunk.longitudes[j]
                if isinstance(_lon, np.ma.core.MaskedConstant):
                    continue

                value = self.__getChunkValueAtIndex(chunk, (i, j))
                if isinstance(value, float) and (math.isnan(value) or value == np.nan):
                    continue

                u = utm.from_latlon(_lat, _lon)
                v = (u[0], u[1], 0.0)
                latlons.append((_lat, _lon))
                utms.append(v)
                indexes.append((i, j))

        tree = None
        if len(latlons) > 0:
            tree = spatial.KDTree(utms)

        chunk.swathIndexing = {
            "tree": tree,
            "latlons": latlons,
            "indexes": indexes
        }

    def __getChunkIndexesForLatLon(self, chunk, lat, lon, xyTolerance):
        foundIndexes = []
        foundLatLons = []

        if "swathIndexing" not in chunk.__dict__:
            self.__buildSwathIndexes(chunk)

        tree = chunk.swathIndexing["tree"]
        if tree is not None:
            indexes = chunk.swathIndexing["indexes"]
            latlons = chunk.swathIndexing["latlons"]
            u = utm.from_latlon(lat, lon)
            coords = np.array([u[0], u[1], 0])
            ball = tree.query_ball_point(coords, xyTolerance)
            for i in ball:
                foundIndexes.append(indexes[i])
                foundLatLons.append(latlons[i])
        return foundIndexes, foundLatLons

    def __getChunkValueAtIndex(self, chunk, index, arrayName=None):

        if arrayName is None or arrayName == "data":
            data_val = chunk.data[0][index[0]][index[1]]
        else:
            data_val = chunk.meta_data[arrayName][0][index[0]][index[1]]
        return data_val.item() if (data_val is not np.ma.masked) and data_val.size == 1 else float('Nan')

    def __getSatNodeForLatLonAndTime(self, chunksByDay, source, lat, lon, searchTime, xyTolerance):
        timeDiff = 86400 * 365 * 1000
        foundNodes = []

        for ts in chunksByDay:
            chunks = chunksByDay[ts]
            if abs((ts * 1000) - searchTime) < timeDiff:
                for chunk in chunks:
                    indexes, latlons = self.__getChunkIndexesForLatLon(chunk, lat, lon, xyTolerance)

                    # for index in indexes:
                    for i in range(0, len(indexes)):
                        index = indexes[i]
                        latlon = latlons[i]
                        sst = None
                        sss = None
                        windSpeed = None
                        windDirection = None
                        windU = None
                        windV = None

                        value = self.__getChunkValueAtIndex(chunk, index)

                        if isinstance(value, float) and (math.isnan(value) or value == np.nan):
                            continue

                        if "GHRSST" in source:
                            sst = value
                        elif "ASCATB" in source:
                            windU = value
                        elif "SSS" in source:  # SMAP
                            sss = value

                        if len(chunks) > 0 and "wind_dir" in chunks[0].meta_data:
                            windDirection = self.__checkNumber(self.__getChunkValueAtIndex(chunk, index, "wind_dir"))
                        if len(chunks) > 0 and "wind_v" in chunks[0].meta_data:
                            windV = self.__checkNumber(self.__getChunkValueAtIndex(chunk, index, "wind_v"))
                        if len(chunks) > 0 and "wind_speed" in chunks[0].meta_data:
                            windSpeed = self.__checkNumber(self.__getChunkValueAtIndex(chunk, index, "wind_speed"))

                        foundNode = {
                            "sea_water_temperature": sst,
                            "sea_water_salinity": sss,
                            "wind_speed": windSpeed,
                            "wind_direction": windDirection,
                            "wind_u": windU,
                            "wind_v": windV,
                            "time": ts,
                            "x": self.__checkNumber(latlon[1]),
                            "y": self.__checkNumber(latlon[0]),
                            "depth": 0,
                            "sea_water_temperature_depth": 0,
                            "source": source,
                            "id": "%s:%s:%s" % (ts, lat, lon)
                        }

                        foundNodes.append(foundNode)
                timeDiff = abs(ts - searchTime)

        return foundNodes

    def __getSatNodeForLatLonAndTime__(self, chunksByDay, source, lat, lon, searchTime):

        timeDiff = 86400 * 365 * 1000
        foundNodes = []

        for ts in chunksByDay:
            chunks = chunksByDay[ts]
            # print chunks
            # ts = calendar.timegm(chunks.start.utctimetuple()) * 1000
            if abs((ts * 1000) - searchTime) < timeDiff:
                value = self.__getValueForLatLon(chunks, lat, lon, arrayName="data")
                value = self.__checkNumber(value)

                # _Really_ don't like doing it this way...

                sst = None
                sss = None
                windSpeed = None
                windDirection = None
                windU = None
                windV = None

                if "GHRSST" in source:
                    sst = value

                if "ASCATB" in source:
                    windU = value

                if len(chunks) > 0 and "wind_dir" in chunks[0].meta_data:
                    windDirection = self.__checkNumber(self.__getValueForLatLon(chunks, lat, lon, arrayName="wind_dir"))
                if len(chunks) > 0 and "wind_v" in chunks[0].meta_data:
                    windV = self.__checkNumber(self.__getValueForLatLon(chunks, lat, lon, arrayName="wind_v"))
                if len(chunks) > 0 and "wind_speed" in chunks[0].meta_data:
                    windSpeed = self.__checkNumber(self.__getValueForLatLon(chunks, lat, lon, arrayName="wind_speed"))

                foundNode = {
                    "sea_water_temperature": sst,
                    "sea_water_salinity": sss,
                    "wind_speed": windSpeed,
                    "wind_direction": windDirection,
                    "wind_uv": {
                        "u": windU,
                        "v": windV
                    },
                    "time": ts,
                    "x": lon,
                    "y": lat,
                    "depth": 0,
                    "sea_water_temperature_depth": 0,
                    "source": source,
                    "id": "%s:%s:%s" % (ts, lat, lon)
                }

                isValidNode = True
                if "ASCATB" in source and windSpeed is None:
                    isValidNode = None

                if isValidNode:
                    foundNodes.append(foundNode)
                timeDiff = abs(ts - searchTime)

        return foundNodes


class InsituDatasetProcessor:
    def __init__(self, primary, datasource, startTime, endTime, bbox, depth_min, depth_max, platforms, timeTolerance,
                 radiusTolerance):
        self.primary = primary
        self.datasource = datasource
        self.startTime = startTime
        self.endTime = endTime
        self.bbox = bbox
        self.depth_min = depth_min
        self.depth_max = depth_max
        self.platforms = platforms
        self.timeTolerance = timeTolerance
        self.radiusTolerance = radiusTolerance

    def start(self):
        def callback(pageData):
            self.primary.processInSitu(pageData, self.radiusTolerance, self.timeTolerance)

        fetchedgeimpl.fetch(self.datasource, self.startTime, self.endTime, self.bbox, self.depth_min, self.depth_max,
                            self.platforms, pageCallback=callback)


class InsituPageProcessor:
    def __init__(self):
        pass