analysis/webservice/algorithms/TimeAvgMap.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.


# distutils: include_dirs = /usr/local/lib/python2.7/site-packages/cassandra


import sys
import numpy as np
from time import time
from webservice.NexusHandler import DEFAULT_PARAMETERS_SPEC
from webservice.algorithms.NexusCalcHandler import NexusCalcHandler
from webservice.webmodel import NexusResults, NoDataException
from netCDF4 import Dataset


# from mpl_toolkits.basemap import Basemap


# @nexus_handler
class TimeAvgMapCalcHandlerImpl(NexusCalcHandler):
    name = "Time Average Map"
    path = "/timeAvgMap"
    description = "Computes a Latitude/Longitude Time Average plot given an arbitrary geographical area and time range"
    params = DEFAULT_PARAMETERS_SPEC
    singleton = True

    def _find_native_resolution(self):
        # Get a quick set of tiles (1 degree at center of box) at 1 time stamp
        midLat = (self._minLat + self._maxLat) / 2
        midLon = (self._minLon + self._maxLon) / 2
        ntiles = 0
        t = self._endTime
        t_incr = 86400
        while ntiles == 0:
            nexus_tiles = self._get_tile_service().get_tiles_bounded_by_box(midLat - 0.5, midLat + 0.5, midLon - 0.5,
                                                                      midLon + 0.5, ds=self._ds, start_time=t - t_incr,
                                                                      end_time=t)
            ntiles = len(nexus_tiles)
            print('find_native_res: got %d tiles' % len(nexus_tiles))
            sys.stdout.flush()
            lat_res = 0.
            lon_res = 0.
            if ntiles > 0:
                for tile in nexus_tiles:
                    if lat_res < 1e-10:
                        lats = tile.latitudes.compressed()
                        if (len(lats) > 1):
                            lat_res = lats[1] - lats[0]
                    if lon_res < 1e-10:
                        lons = tile.longitudes.compressed()
                        if (len(lons) > 1):
                            lon_res = lons[1] - lons[0]
                    if (lat_res >= 1e-10) and (lon_res >= 1e-10):
                        break
            if (lat_res < 1e-10) or (lon_res < 1e-10):
                t -= t_incr

        self._latRes = lat_res
        self._lonRes = lon_res

    def _find_global_tile_set(self):
        ntiles = 0
        t = self._endTime
        t_incr = 86400
        while ntiles == 0:
            nexus_tiles = self._get_tile_service().get_tiles_bounded_by_box(self._minLat, self._maxLat, self._minLon,
                                                                      self._maxLon, ds=self._ds, start_time=t - t_incr,
                                                                      end_time=t)
            ntiles = len(nexus_tiles)
            print('find_global_tile_set got %d tiles' % ntiles)
            sys.stdout.flush()
            t -= t_incr
        return nexus_tiles

    def _prune_tiles(self, nexus_tiles):
        del_ind = np.where([np.all(tile.data.mask) for tile in nexus_tiles])[0]
        for i in np.flipud(del_ind):
            del nexus_tiles[i]

    # @staticmethod
    # def _map(tile_in):
    def _map(self, tile_in):
        print('Started tile %s' % tile_in.section_spec)
        print('tile lats = ', tile_in.latitudes)
        print('tile lons = ', tile_in.longitudes)
        print('tile = ', tile_in.data)
        sys.stdout.flush()
        lats = tile_in.latitudes
        lons = tile_in.longitudes
        if len(lats > 0) and len(lons > 0):
            min_lat = np.ma.min(lats)
            max_lat = np.ma.max(lats)
            min_lon = np.ma.min(lons)
            max_lon = np.ma.max(lons)
            good_inds_lat = np.where(lats.mask == False)
            good_inds_lon = np.where(lons.mask == False)
            min_y = np.min(good_inds_lat)
            max_y = np.max(good_inds_lat)
            min_x = np.min(good_inds_lon)
            max_x = np.max(good_inds_lon)
            tile_inbounds_shape = (max_y - min_y + 1, max_x - min_x + 1)
            days_at_a_time = 90
            t_incr = 86400 * days_at_a_time
            avg_tile = np.ma.array(np.zeros(tile_inbounds_shape,
                                            dtype=np.float64))
            cnt_tile = np.ma.array(np.zeros(tile_inbounds_shape,
                                            dtype=np.uint32))
            t_start = self._startTime
            while t_start <= self._endTime:
                t_end = min(t_start + t_incr, self._endTime)
                t1 = time()
                print('nexus call start at time %f' % t1)
                sys.stdout.flush()
                nexus_tiles = self._get_tile_service().get_tiles_bounded_by_box(min_lat - self._latRes / 2,
                                                                          max_lat + self._latRes / 2,
                                                                          min_lon - self._lonRes / 2,
                                                                          max_lon + self._lonRes / 2, ds=self._ds,
                                                                          start_time=t_start, end_time=t_end)
                t2 = time()
                print('nexus call end at time %f' % t2)
                print('secs in nexus call: ', t2 - t1)
                sys.stdout.flush()
                self._prune_tiles(nexus_tiles)
                print('t %d to %d - Got %d tiles' % (t_start, t_end,
                                                     len(nexus_tiles)))
                sys.stdout.flush()
                for tile in nexus_tiles:
                    tile.data.data[:, :] = np.nan_to_num(tile.data.data)
                    avg_tile.data[:, :] += tile.data[0,
                                           min_y:max_y + 1,
                                           min_x:max_x + 1]
                    cnt_tile.data[:, :] += (~tile.data.mask[0,
                                             min_y:max_y + 1,
                                             min_x:max_x + 1]).astype(np.uint8)
                t_start = t_end + 1

            print('cnt_tile = ', cnt_tile)
            cnt_tile.mask = ~(cnt_tile.data.astype(bool))
            avg_tile.mask = cnt_tile.mask
            avg_tile /= cnt_tile
            print('Finished tile %s' % tile_in.section_spec)
            print('Tile avg = ', avg_tile)
            sys.stdout.flush()
        else:
            avg_tile = None
            min_lat = None
            max_lat = None
            min_lon = None
            max_lon = None
            print('Tile %s outside of bounding box' % tile_in.section_spec)
            sys.stdout.flush()
        return (avg_tile, min_lat, max_lat, min_lon, max_lon)

    def _lat2ind(self, lat):
        return int((lat - self._minLatCent) / self._latRes)

    def _lon2ind(self, lon):
        return int((lon - self._minLonCent) / self._lonRes)

    def _create_nc_file(self, a):
        print('a=', a)
        print('shape a = ', a.shape)
        sys.stdout.flush()
        lat_dim, lon_dim = a.shape
        rootgrp = Dataset("tam.nc", "w", format="NETCDF4")
        rootgrp.createDimension("lat", lat_dim)
        rootgrp.createDimension("lon", lon_dim)
        rootgrp.createVariable("TRMM_3B42_daily_precipitation_V7", "f4",
                               dimensions=("lat", "lon",))
        rootgrp.createVariable("lat", "f4", dimensions=("lat",))
        rootgrp.createVariable("lon", "f4", dimensions=("lon",))
        rootgrp.variables["TRMM_3B42_daily_precipitation_V7"][:, :] = a
        rootgrp.variables["lat"][:] = np.linspace(self._minLatCent,
                                                  self._maxLatCent, lat_dim)
        rootgrp.variables["lon"][:] = np.linspace(self._minLonCent,
                                                  self._maxLonCent, lon_dim)
        rootgrp.close()

    def calc(self, computeOptions, **args):
        """

        :param computeOptions: StatsComputeOptions
        :param args: dict
        :return:
        """

        self._minLat = float(computeOptions.get_min_lat())
        self._maxLat = float(computeOptions.get_max_lat())
        self._minLon = float(computeOptions.get_min_lon())
        self._maxLon = float(computeOptions.get_max_lon())
        self._ds = computeOptions.get_dataset()[0]
        self._startTime = computeOptions.get_start_time()
        self._endTime = computeOptions.get_end_time()

        self._find_native_resolution()
        print('Using Native resolution: lat_res=%f, lon_res=%f' % (self._latRes, self._lonRes))
        self._minLatCent = self._minLat + self._latRes / 2
        self._minLonCent = self._minLon + self._lonRes / 2
        nlats = int((self._maxLat - self._minLatCent) / self._latRes) + 1
        nlons = int((self._maxLon - self._minLonCent) / self._lonRes) + 1
        self._maxLatCent = self._minLatCent + (nlats - 1) * self._latRes
        self._maxLonCent = self._minLonCent + (nlons - 1) * self._lonRes
        print('nlats=', nlats, 'nlons=', nlons)
        print('center lat range = %f to %f' % (self._minLatCent,
                                               self._maxLatCent))
        print('center lon range = %f to %f' % (self._minLonCent,
                                               self._maxLonCent))
        sys.stdout.flush()
        a = np.zeros((nlats, nlons), dtype=np.float64, order='C')

        nexus_tiles = self._find_global_tile_set()
        # print 'tiles:'
        # for tile in nexus_tiles:
        #     print tile.granule
        #     print tile.section_spec
        #     print 'lat:', tile.latitudes
        #     print 'lon:', tile.longitudes

        #                                                          nexus_tiles)
        if len(nexus_tiles) == 0:
            raise NoDataException(reason="No data found for selected timeframe")

        print('Initially found %d tiles' % len(nexus_tiles))
        sys.stdout.flush()
        self._prune_tiles(nexus_tiles)
        print('Pruned to %d tiles' % len(nexus_tiles))
        sys.stdout.flush()
        # for tile in nexus_tiles:
        #    print 'lats: ', tile.latitudes.compressed()
        #    print 'lons: ', tile.longitudes.compressed()

        avg_tiles = list(map(self._map, nexus_tiles))
        print('shape a = ', a.shape)
        sys.stdout.flush()
        # The tiles below are NOT Nexus objects.  They are tuples
        # with the time avg map data and lat-lon bounding box.
        for tile in avg_tiles:
            if tile is not None:
                (tile_data, tile_min_lat, tile_max_lat,
                 tile_min_lon, tile_max_lon) = tile
                print('shape tile_data = ', tile_data.shape)
                print('tile data mask = ', tile_data.mask)
                sys.stdout.flush()
                if np.any(np.logical_not(tile_data.mask)):
                    y0 = self._lat2ind(tile_min_lat)
                    y1 = self._lat2ind(tile_max_lat)
                    x0 = self._lon2ind(tile_min_lon)
                    x1 = self._lon2ind(tile_max_lon)
                    print('writing tile lat %f-%f, lon %f-%f, map y %d-%d, map x %d-%d' % \
                          (tile_min_lat, tile_max_lat,
                           tile_min_lon, tile_max_lon, y0, y1, x0, x1))
                    sys.stdout.flush()
                    a[y0:y1 + 1, x0:x1 + 1] = tile_data
                else:
                    print('All pixels masked in tile lat %f-%f, lon %f-%f, map y %d-%d, map x %d-%d' % \
                          (tile_min_lat, tile_max_lat,
                           tile_min_lon, tile_max_lon, y0, y1, x0, x1))
                    sys.stdout.flush()

        self._create_nc_file(a)

        return TimeAvgMapResults(results={}, meta={}, computeOptions=computeOptions)


class TimeAvgMapResults(NexusResults):

    def __init__(self, results=None, meta=None, computeOptions=None):
        NexusResults.__init__(self, results=results, meta=meta, stats=None, computeOptions=computeOptions)