| # 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) |
