diff --git a/bin/map_main.py b/bin/map_main.py index b174da71dc2f741c583a96b84c2eaee37585491a..a4c204ae6ea42ed0bffcf91b0c2433128606f5c9 100755 --- a/bin/map_main.py +++ b/bin/map_main.py @@ -29,7 +29,7 @@ mpl_logger = logging.getLogger('matplotlib') mpl_logger.setLevel('CRITICAL') from mapgenerator.plotting.config import ArgumentParser -from mapgenerator.plotting import plotmap +from mapgenerator.plotting import plotmap, timeseries import sys from datetime import datetime, timedelta @@ -76,7 +76,7 @@ class MapGenerator: if res['pltype'] == 'cross': plotmap.PlotCross(loglevel).plot(**res) elif res['pltype'] == 'timeseries': - plotmap.PlotSeries(loglevel).plot(**res) + timeseries.PlotSeries(loglevel).plot(**res) elif res['pltype'] == 'map': plotmap.PlotMap(loglevel).plot(**res) else: diff --git a/mapgenerator/plotting/config.py b/mapgenerator/plotting/config.py index 35cf6c3b7287c23b06338c558c6a48f8b537c005..c18887ae9cbdb79e8f278596fd4ffa780c6600eb 100644 --- a/mapgenerator/plotting/config.py +++ b/mapgenerator/plotting/config.py @@ -13,11 +13,11 @@ from mapgenerator import mg_exceptions import logging -logging.basicConfig(level=logging.DEBUG) +logging.basicConfig(level=logging.WARNING) log = logging.getLogger(__name__) -#def list_callback(option, opt_str, value, parser): +# def list_callback(option, opt_str, value, parser): # setattr(parser.values, option.dest, value.split(',')) # print("**************", value, type(value)) # print('\toption:', repr(option)) @@ -36,11 +36,11 @@ class ArgumentParser(object): self.parser = argparse.ArgumentParser(description='') self.parser.add_argument('-V', '--version', action='version', version=__version__, - help="returns MapGenerator version number and exit") - self.parser.add_argument('--config', #is_config_file=True, + help="returns version and exit") + self.parser.add_argument('--config', dest="config", help='specifies the config file to read' - ) #required=False) + ) self.parser.add_argument("--section", dest="section", @@ -59,23 +59,27 @@ class ArgumentParser(object): # self.parser.add_argument("--loglevel", dest="loglevel", - choices=['DEBUG', 'INFO', 'WARNING', 'ERROR', + choices=['DEBUG', + 'INFO', + 'WARNING', + 'ERROR', 'CRITICAL'], help="logging levelbug (default=WARNING)") self.parser.add_argument("--indir", dest="indir", - help="directory where is/are stored the input file(s)") + help="directory for input file(s)") self.parser.add_argument("--outdir", dest="outdir", - help="directory where is/are to be stored the output file(s)") + help="directory for output file(s)") self.parser.add_argument("--file_format", dest="filefmt", - help="format of the output file(s) (default=png)") + help="outputs format (default=png)") self.parser.add_argument("--timesteps", dest="timesteps", nargs="+", - help="timesteps to plot (list or 'all' - default=[0])") + help="""timesteps (list or 'all' - + default=[0])""") self.parser.add_argument("--bounds", dest="bounds", nargs="+", @@ -85,7 +89,9 @@ class ArgumentParser(object): help="colorbar ticks list (default=None)") self.parser.add_argument("--colors", dest="colors", - help="colors list or colormap name (default='jet'), for a list check here: https://matplotlib.org/users/colormaps.html") + help="""colors list or colormap name + (default='jet'), for a list check here: + https://matplotlib.org/users/colormaps.html""") self.parser.add_argument("--lon", dest="lon", nargs="+", @@ -115,7 +121,7 @@ class ArgumentParser(object): help="extend colormap extremes (default='neither')") self.parser.add_argument("--resolution", dest="resolution", - choices=['c', 'l', 'i', 'h', 'f'], + choices=['110m', '50m', '10m'], help="details level of coastlines/coutries") self.parser.add_argument("--dpi", dest="dpi", @@ -144,9 +150,6 @@ class ArgumentParser(object): self.parser.add_argument("--contours", dest="contours", help="additional contour variable") - self.parser.add_argument("--keep_aspect", - dest="keep_aspect", - help="preserve map aspect ratio when drawing (default=False)") self.parser.add_argument("--shapefiles", dest="shapefiles", help="additional shape files") @@ -162,12 +165,12 @@ class ArgumentParser(object): self.parser.add_argument("--ysize", dest="ysize", help="proportion of Y dimension (between 0 and 1)") + self.parser.add_argument("--orig_projection", + dest="orig_projection", + help="origin map projection (default=PlateCarree)") self.parser.add_argument("--projection", dest="projection", - help="map projection to draw (default=cyl)") - self.parser.add_argument("--area_thresh", - dest="area_thresh", - help="area threshold in surface drawing (including lakes)") + help="map projection to draw (default=PlateCarree)") self.parser.add_argument("--title", dest="title", @@ -189,7 +192,8 @@ class ArgumentParser(object): help="dimension name and operation in the form DIMENSION,OPERATION (mean, max, min, ...) as documented in the numpy reference guide. If OPERATION is a number it will extract the correspondent element") self.parser.add_argument("--alpha", dest="alpha", - help="alpha value for colormap netween 0 and 1 (default=1)") + help="""alpha value for colormap netween 0 + and 1 (default=1)""") ############ self.parser.add_argument("--img_template", dest="img_template", @@ -200,9 +204,9 @@ class ArgumentParser(object): self.parser.add_argument("--nomap", dest="nomap", help="hide map") -# self.parser.add_argument("--kml", -# dest="kml", -# help="generate KML. String with base url (default='')") + self.parser.add_argument("--kml", + dest="kml", + help="generate KML (default=False)") self.parser.add_argument("--kmz", dest="kmz", help="generate KMZ (default=False)") @@ -215,7 +219,8 @@ class ArgumentParser(object): help="background orography (default=None)") self.parser.add_argument("--logo", dest="logo", - help="put a logo (image path, x, y - default=None)") + help="""put a logo (image path, x, y - + default=None)""") # self.parser.add_argument("--noruntime", # dest="noruntime", # help="no runtime") @@ -264,18 +269,18 @@ class ArgumentParser(object): # help="max") except Exception as e: - log.error('Unhandled exception on MapGenerator: %s' % e, exc_info=True) + log.error('Unhandled exception on MapGenerator: %s' % e, + exc_info=True) - #----------------------------------------------------------------------- + # ----------------------------------------------------------------------- # Parse arguments and preprocess - #----------------------------------------------------------------------- + # ----------------------------------------------------------------------- def parse_args(self, args=None): """ Parse arguments given to an executable :param args: """ try: - #return self.parser.parse_args(args) return self.parser.parse_args(args) except Exception as e: print(e) diff --git a/mapgenerator/plotting/definitions.py b/mapgenerator/plotting/definitions.py index efcc3987251c57073a0d5902878e203613699c06..8607b5b74b4f9c43520ece80fd6b8f4293f38a9f 100644 --- a/mapgenerator/plotting/definitions.py +++ b/mapgenerator/plotting/definitions.py @@ -5,8 +5,8 @@ # @author: see AUTHORS file from netCDF4 import Dataset as nc -#from netCDF4 import MFDataset as mnc -#from grads import GaLab +# from netCDF4 import MFDataset as mnc +# from grads import GaLab import pandas as pd import numpy as np import numpy.ma as ma @@ -19,34 +19,34 @@ import copy import logging -#logging.basicConfig(level=logging.WARNING) -log = logging.getLogger(__name__) +# logging.basicConfig(level=logging.INFO) +LOG = logging.getLogger(__name__) def is_grid_regular(lon, lat): """ Check spacing if uniform. Returns True if grid is regular """ - #print("Checking if lon is strictly_increasing ...") - if not np.array([strictly_increasing(lon[x]) \ - for x in np.arange(lon.shape[0])]).all(): + # print("Checking if lon is strictly_increasing ...") + if not np.array([strictly_increasing(lon[x]) for x in + np.arange(lon.shape[0])]).all(): return False - #print("Checking if lat is strictly_increasing ...") - if not np.array([strictly_increasing(lat[:, y]) \ - for y in np.arange(lat.shape[1])]).all(): + # print("Checking if lat is strictly_increasing ...") + if not np.array([strictly_increasing(lat[:, y]) for y in + np.arange(lat.shape[1])]).all(): return False - #print("Checking lon spacing ...") - is_x_reg = np.array([(np.spacing(lon[x]) == np.spacing(lon[x+1])).all() \ - for x in np.arange(lon.shape[0]-1)]).all() + # print("Checking lon spacing ...") + is_x_reg = np.array([(np.spacing(lon[x]) == np.spacing(lon[x+1])).all() + for x in np.arange(lon.shape[0]-1)]).all() - print(is_x_reg) + LOG.debug(is_x_reg) if not is_x_reg: return is_x_reg - #print("Checking lat spacing ...") - return np.array([(np.spacing(lat[:, y]) == np.spacing(lat[:, y+1])).all() \ - for y in np.arange(lat.shape[1]-1)]).all() + # print("Checking lat spacing ...") + return np.array([(np.spacing(lat[:, y]) == np.spacing(lat[:, y+1])).all() + for y in np.arange(lat.shape[1]-1)]).all() def strictly_increasing(vector): @@ -60,12 +60,13 @@ def do_interpolation(data, lon, lat): avoid """ # lon as second dimension # lat as first dimension - log.info("Irregular grid, performing interpolation ...") + LOG.info("Irregular grid, performing interpolation ...") lat_s = lat.shape[0] lon_s = lon.shape[1] reglon = np.linspace(lon.min(), lon.max(), lon_s) reglat = np.linspace(lat.min(), lat.max(), lat_s) - log.info("LAT: %s, LON: %s, DATA: %s" % (reglon.shape, reglat.shape, data.shape)) + LOG.info("LAT: %s, LON: %s, DATA: %s", reglon.shape, reglat.shape, + data.shape) regx, regy = np.meshgrid(reglon, reglat) result = np.empty(data.shape)*np.nan if ma.is_masked(data): @@ -73,17 +74,20 @@ def do_interpolation(data, lon, lat): # assuming that lat and lon are the latest dimensions if len(data.shape) == 2: - result = griddata((lon.ravel(), lat.ravel()), data.ravel(), (regx, regy)) + result = griddata((lon.ravel(), lat.ravel()), data.ravel(), (regx, + regy)) elif len(data.shape) == 3: for i in range(data.shape[0]): - tmp = griddata((lon.ravel(), lat.ravel()), data[i].ravel(), (regx, regy)) + tmp = griddata((lon.ravel(), lat.ravel()), data[i].ravel(), (regx, + regy)) result[i] = tmp.reshape(result[i].shape) elif len(data.shape) == 4: for i in range(data.shape[0]): for j in range(data.shape[1]): - tmp = griddata((lon.ravel(), lat.ravel()), data[i, j].ravel(), (regx, regy)) + tmp = griddata((lon.ravel(), lat.ravel()), data[i, j].ravel(), + (regx, regy)) result[i, j, :, :] = tmp.reshape(result[i, j, :, :].shape) - #elif len(data.shape) == 5: + # elif len(data.shape) == 5: return np.ma.masked_where(np.isnan(result), result), regx, regy @@ -91,32 +95,35 @@ def extract_coords(self, var, coord='lons'): """ Given a variable and a coordinate name (lons, lats) extract values related to given coordinate interval """ var = var[:].squeeze() - log.info("VAR SHAPE %s" % str(var.shape)) + LOG.info("VAR SHAPE %s" % str(var.shape)) if not self.subsetting: - log.info("Not subsetting") + LOG.info("Not subsetting") return var[:], None cur_coord = self.__getattribute__(coord) - log.info("Coord: %s, cur_coord: %s" % (coord, str(cur_coord))) - if cur_coord:# and len(cur_coord) in (2,3): + LOG.info("Coord: %s, cur_coord: %s" % (coord, str(cur_coord))) + if cur_coord: # and len(cur_coord) in (2,3): if len(var.shape) == 2: if coord == 'lons': var1 = var[0] else: var1 = var[:, 0] - else: var1 = var + else: + var1 = var if strictly_increasing(var1[:]): - idxs = np.where((var1[:] >= cur_coord[0]) & \ - (var1[:] <= cur_coord[-1])) - log.info("idxs %s %s %s %s" % (str(idxs), var1.shape, len(idxs), idxs[0].shape)) + idxs = np.where((var1[:] >= cur_coord[0]) & + (var1[:] <= cur_coord[-1])) + LOG.info("idxs %s %s %s %s", str(idxs), var1.shape, len(idxs), + idxs[0].shape) return var1[idxs], idxs[0] - log.info("Not subsetting") + LOG.info("Not subsetting") return var[:], None def parse_parameters_list(plist): """ Parser for list parameter to be able to enter intervals and number of steps (examples: 0,1,2 0-6,2 ...""" - if isinstance(plist, list) or isinstance(plist, tuple) or isinstance(plist, np.ndarray): + if isinstance(plist, list) or isinstance(plist, tuple) or \ + isinstance(plist, np.ndarray): if len(plist) == 1: plist = plist[0] else: @@ -134,21 +141,21 @@ def parse_parameters_list(plist): steps = (end-start)/10 else: steps = float(steps) - return list(np.arange(start, end+int(steps), steps)) #, endpoint=True)) + return list(np.arange(start, end+int(steps), steps)) else: res = eval(str(plist).replace('m', '-')) if isinstance(res, int) or isinstance(res, float): - return [res,] + return [res] return list(res) def parse_parameter(var, val): """ Parameter parser """ # parse numeric list arguments - if var in ('lon', 'lat', 'timesteps', 'bounds'): + if var in ('lon', 'lat', 'timesteps', 'bounds', 'ticks'): if var == 'timesteps' and val in ('all', ['all']): return 'all' - if val != None: + if val is not None: return parse_parameters_list(val) # parse float arguments @@ -159,7 +166,7 @@ def parse_parameter(var, val): # parse simple list arguments elif var in ('coordsopts', 'coastsopts', 'countropts'): if val and isinstance(val, str): - #log.info(val, type(val)) + LOG.debug("%s: %s :: %s", var, val, type(val)) return val.split(',') return val @@ -169,29 +176,28 @@ def parse_parameters(params): """ Parser for parameters """ ret = copy.deepcopy(params) for var in params: - #log.info("{}: {}".format(var, params[var])) + LOG.debug("%s: %s", var, params[var]) val = parse_parameter(var, params[var]) ret[var] = val - #log.info("{}: {}".format(var, val)) + LOG.debug("%s: %s", var, val) return ret def parse_path(directory, dfile): """ Path parser """ - #print "Opening data file", f + # print "Opening data file", f if os.path.isabs(dfile): return dfile else: - #log.info("Input: %s", f) + LOG.debug("Input: %s", dfile) return os.path.join(directory, dfile) class MapData(object): """Wrapper class to handle all the data needed to draw a map""" - def __init__(self, loglevel='WARNING', **kwargs): - log.setLevel(loglevel) - log.info("map_data initializer") + def __init__(self, **kwargs): + LOG.info("map_data initializer") self.reset() self.max_reset() # Process optional arguments @@ -209,11 +215,9 @@ class MapData(object): self.max_data = None def set_grid_data(self, data): - #print "Setting grid data" self.grid_data = data def set_map_data(self, data): - #print "Setting map data" self.map_data = data def set_max_data(self, data): @@ -247,38 +251,26 @@ class MapData(object): return self.scatter_data def has_map_data(self): - return (self.map_data != None) + return (self.map_data is not None) def has_wind_data(self): - return (self.wind_data != None) + return (self.wind_data is not None) def has_contour_data(self): return self.contour_data is not None def has_max_data(self): - return (self.max_data != None) + return (self.max_data is not None) def has_scatter_data(self): - return (self.scatter_data != None) + return (self.scatter_data is not None) class MapDataHandler(object): -# ncdf = None -# data = None -# srcs = [] -# srcgaps = [] -# windsrc = [] -# windopts = [] -# indir = None -# lats = [] -# lons = [] -# max = None -# transf = False def __init__(self, sources, srcvars, indir, srcgaps, windsource, windopts, lats, lons, transf, subsetting, dimension, **kwargs): - #self.galab = GaLab(Bin='grads', Window=False, Echo=False) - self.ncdf = [[], []] # list of lists, srcfiles and windfiles + self.ncdf = [[], []] # list of lists, srcfiles and windfiles self.srcs = sources self.srcvars = srcvars self.srcgaps = srcgaps @@ -302,14 +294,14 @@ class MapDataHandler(object): if ('winds' in kwargs) \ and (kwargs['winds']['src']) \ and (kwargs['winds']['opts']): - log.info("Will do winds") + LOG.info("Will do winds") self.do_winds = True self.windsrc = kwargs['winds']['src'] self.windopts = kwargs['winds']['opts'] if ('contours' in kwargs) \ and (kwargs['contours']['var']): - log.info("Will do contours") + LOG.info("Will do contours") self.do_contours = True self.contoursrc = self.srcs[0] self.contourvar = kwargs['contours']['var'] @@ -319,7 +311,6 @@ class MapDataHandler(object): if vcond: self.varconds = eval(vcond) - #self.galab.cmd("reinit") # Open sources for sfile in self.srcs: fpath = parse_path(self.indir, sfile) @@ -330,7 +321,8 @@ class MapDataHandler(object): if vname in ('longitude', 'lon', 'nav_lon'): x, x_idxs = extract_coords(self, fin.variables[vname]) elif vname in ('latitude', 'lat', 'nav_lat'): - y, y_idxs = extract_coords(self, fin.variables[vname], coord='lats') + y, y_idxs = extract_coords(self, fin.variables[vname], + coord='lats') self.grid = (x, y) if x_idxs is not None and y_idxs is not None: @@ -338,9 +330,9 @@ class MapDataHandler(object): # Open wind source if available for w in self.windsrc: - log.info("WINDS: %s" % w) + LOG.info("WINDS: %s", w) wpath = parse_path(self.indir, w) - log.info("WPATH: %s" % wpath) + LOG.info("WPATH: %s", wpath) self.ncdf[1].append(nc(wpath)) # if self.transf: @@ -353,12 +345,11 @@ class MapDataHandler(object): f.close() for w in self.ncdf[1]: w.close() - #print "Resetting GraDs engine" def dim_operation(self, srcvar, srcfile): """ Operations over dimensions """ dim, oper = self.dimension.split(',') - log.info("Dimension %s, operation or index/value: %s" % (dim, oper)) + LOG.info("Dimension %s, operation or index/value: %s" % (dim, oper)) dimidx = srcvar.dimensions.index(dim) # string try: @@ -377,21 +368,22 @@ class MapDataHandler(object): def get_data_for_tstep(self, step): """ Get data per timestep """ self.data.reset() - #print "Getting map data for step ", step + # print "Getting map data for step ", step self.data.set_map_data(self.get_map_data_for_tstep(step)) if self.do_winds: - log.info("Getting wind data for step '%s'" % step) + LOG.info("Getting wind data for step '%s'" % step) self.data.set_wind_data(self.get_wind_data_for_tstep(step)) if self.do_contours: - log.info("Getting contour data for step '%s'" % step) + LOG.info("Getting contour data for step '%s'" % step) self.data.set_contour_data(self.get_contour_data_for_tstep(step)) # Calculate maximum over data # FIXME should this step be optional? # if (self.data.has_max_data()): -# self.data.set_max_data(maximum(self.data.get_max_data(), self.data.get_map_data())) +# self.data.set_max_data(maximum(self.data.get_max_data(), +# self.data.get_map_data())) # else: # self.data.set_max_data(self.data.get_map_data()) # self.last_tstep = step @@ -401,7 +393,8 @@ class MapDataHandler(object): '''Get current max data. NOTE: A call to this method will reset current max data, so it can be called only once! ''' - log.info("Retrieving max (last calculated tstep: %d); this will reset current MAX data!!" % self.last_tstep) + LOG.info("""Retrieving max (last calculated tstep: %d); this will reset + current MAX data!!""", self.last_tstep) max_data = MapData(map_data=self.data.get_max_data()) self.data.max_reset() return max_data @@ -410,10 +403,11 @@ class MapDataHandler(object): """ Get map data for tstep """ map_data = [] for src_n in range(len(self.srcs)): - tstep = step + int(self.srcgaps[src_n]) + tstep = int(step + self.srcgaps[src_n]) + LOG.info("TSTEP: %s", tstep) srcfile = self.ncdf[0][src_n] - #TODO zoom on lat lon according to the domain - #FIXME variable operations + # TODO zoom on lat lon according to the domain + # FIXME variable operations # if self.varconds and len(self.varconds)>=src_n: # varcond = self.varconds[src_n] # cond = varcond[0] @@ -443,29 +437,36 @@ class MapDataHandler(object): srcvar = self.dim_operation(srcvar, srcfile) if self.grid_idxs: # if no expression - #log.info("SRCVAR: %s - GRID_IDXS: %s" % (str(srcvar.shape), str(self.grid_idxs))) - values = srcvar[tstep, :, :][self.grid_idxs[1], :][:, self.grid_idxs[0]] + # LOG.info("SRCVAR: %s - GRID_IDXS: %s" % + # (str(srcvar.shape), str(self.grid_idxs))) + values = srcvar[tstep, :, :][self.grid_idxs[1], + :][:, self.grid_idxs[0]] else: values = srcvar[tstep, :, :] - log.info("1. SRCVAR: %s - GRID_IDXS: %s" % (str(values.shape), str(self.grid_idxs))) + LOG.info("1. SRCVAR: %s - GRID_IDXS: %s", str(values.shape), + str(self.grid_idxs)) else: # otherwise search for a variable of the dataset into the # expression given for vname in srcfile.variables: if vname in expvar: expvar = expvar.replace(vname, 'newvar') + LOG.info("EXPVAR: %s", expvar) srcvar = srcfile.variables[vname] + LOG.info("SRCVAR: %s", srcvar) if self.dimension: srcvar = self.dim_operation(srcvar, srcfile) if self.grid_idxs: - newvar = srcvar[tstep, :, :][self.grid_idxs[1], :][:, self.grid_idxs[0]] + newvar = srcvar[tstep, :, :][ + self.grid_idxs[1], :][:, self.grid_idxs[0]] else: newvar = srcvar[tstep, :, :] values = eval(expvar) break - log.info("2. SRCVAR: %s - GRID_IDXS: %s" % (str(values.shape), str(self.grid_idxs))) + LOG.info("2. SRCVAR: %s - GRID_IDXS: %s", str(values.shape), + str(self.grid_idxs)) if map_data: - #exp_res = self.ncdf.exp(EXPVAR) + # exp_res = self.ncdf.exp(EXPVAR) map_data.append(values) else: map_data = [values] @@ -475,7 +476,7 @@ class MapDataHandler(object): def get_wind_data_for_tstep(self, tstep): """ Get wind data for tstep """ wind_data = {'u': None, 'v': None} - log.info("WINDOPTS: %s" % str(self.windopts)) + LOG.info("WINDOPTS: %s" % str(self.windopts)) srcfile = self.ncdf[1][0] srcvar_u = srcfile.variables[self.windopts[0]] srcvar_v = srcfile.variables[self.windopts[1]] @@ -483,8 +484,10 @@ class MapDataHandler(object): srcvar_u = self.dim_operation(srcvar_u, srcfile) srcvar_v = self.dim_operation(srcvar_v, srcfile) if self.grid_idxs: - wind_data['u'] = srcvar_u[tstep, :, :][self.grid_idxs[1], :][:, self.grid_idxs[0]] - wind_data['v'] = srcvar_v[tstep, :, :][self.grid_idxs[1], :][:, self.grid_idxs[0]] + wind_data['u'] = srcvar_u[tstep, :, :][self.grid_idxs[1], + :][:, self.grid_idxs[0]] + wind_data['v'] = srcvar_v[tstep, :, :][self.grid_idxs[1], + :][:, self.grid_idxs[0]] else: wind_data['u'] = srcvar_u[tstep, :, :] wind_data['v'] = srcvar_v[tstep, :, :] @@ -495,33 +498,35 @@ class MapDataHandler(object): def get_contour_data_for_tstep(self, tstep): """ Get contour data for tstep """ cdata = None - log.info("Contours variable: %s" % self.contourvar) + LOG.info("Contours variable: %s" % self.contourvar) srcfile = self.ncdf[0][0] srcvar = srcfile.variables[self.contourvar] if self.dimension: srcvar = self.dim_operation(srcvar, srcfile) if self.grid_idxs: - cdata = srcvar[tstep, :, :][self.grid_idxs[1], :][:, self.grid_idxs[0]] + cdata = srcvar[tstep, :, :][self.grid_idxs[1], + :][:, self.grid_idxs[0]] else: cdata = srcvar[tstep, :, :] return cdata def get_scatter_data_for_tstep(self, tstep): """ """ - log.info("Scatter variable: %s" % self.scatter) + LOG.info("Scatter variable: %s" % self.scatter) srcfile = self.ncdf[0][0] srcvar = srcfile.variables[self.contourvar] if self.dimension: srcvar = self.dim_operation(srcvar, srcfile) if self.grid_idxs: - cdata = srcvar[tstep, :, :][self.grid_idxs[1], :][:, self.grid_idxs[0]] + cdata = srcvar[tstep, :, :][self.grid_idxs[1], + :][:, self.grid_idxs[0]] else: cdata = srcvar[tstep, :, :] return cdata def get_dims(self): - log.info("Getting dimensions from file '%s'" % self.srcs[0]) - #self.ncdf.cmd("set dfile 1") + LOG.info("Getting dimensions from file '%s'" % self.srcs[0]) + # self.ncdf.cmd("set dfile 1") for vname in self.ncdf[0][0].variables: if vname in ('time', 'time_counter'): tvar = vname @@ -546,7 +551,7 @@ class MapDataHandler(object): cdate = {(np.where(tdata == i)[0][0], (dtime + timedelta(seconds=int(i))).strftime(tfmt)) for i in tdata} - else: #typ == 'hours': + else: # typ == 'hours': cdate = {(np.where(tdata == i)[0][0], (dtime + timedelta(hours=int(i))).strftime(tfmt)) for i in tdata} @@ -558,7 +563,7 @@ class DataFrameHandler(object): def __init__(self, dframe, separator=',', ): if isinstance(dframe, str): - log.info("opening '%s'" % dframe) + LOG.debug("opening '%s'" % dframe) self.reader = pd.read_csv(dframe, sep=separator, parse_dates=True) else: self.reader = dframe @@ -566,10 +571,10 @@ class DataFrameHandler(object): def filter(self, **kwargs): """ Filtering dataframes """ + latcol = 'latcol' in kwargs and kwargs['latcol'] or 'lat' + loncol = 'loncol' in kwargs and kwargs['loncol'] or 'lon' + colcol = 'colcol' in kwargs and kwargs['colcol'] or 'color' datecol = 'datecol' in kwargs and kwargs['datecol'] or 'date' - latcol = 'latcol' in kwargs and kwargs['latcol'] or 'lat' - loncol = 'loncol' in kwargs and kwargs['loncol'] or 'lon' - colcol = 'colcol' in kwargs and kwargs['colcol'] or 'color' sizecol = 'sizecol' in kwargs and kwargs['sizecol'] or 'size' ret = {} @@ -585,9 +590,12 @@ class DataFrameHandler(object): gdf = self.reader.groupby(datecol) for gname in gdf.groups: group = gdf.get_group(gname) - dtime = datetime.strptime(gname, "%Y-%m-%d %H:%M:00") + if isinstance(gname, str): + dtime = datetime.strptime(gname, "%Y-%m-%d %H:%M:00") + else: + dtime = gname.to_pydatetime() ret[dtime] = group[[loncol, latcol, colcol, sizecol]] - #log.info('filtering on date: %s and hour %s' % (date, hour)) + # LOG.info('filtering on date: %s and hour %s', date, hour) return ret @@ -602,9 +610,9 @@ class DataTransform(object): Colors and bounds must have the same size''' - log.info('bounds: %s' % bounds) - log.info('colors: %s' % colors) - log.info('values: %s' % values) + LOG.info('bounds: %s' % bounds) + LOG.info('colors: %s' % colors) + LOG.info('values: %s' % values) if len(colors) != len(bounds) + 1: raise Exception("Wrong size for colors/bounds") @@ -612,31 +620,31 @@ class DataTransform(object): out = [] for idx in range(0, len(values)): - log.info("processing %s" % str(values[idx])) + LOG.info("processing %s" % str(values[idx])) if values[idx] < bounds[0]: out.append(colors[0]) elif values[idx] >= bounds[len(bounds) - 1]: out.append(colors[len(colors) - 1]) else: for bidx in range(0, len(bounds) - 1): - if values[idx] >= bounds[bidx] and values[idx] < bounds[bidx + 1]: + if values[idx] >= bounds[bidx] and \ + values[idx] < bounds[bidx + 1]: out.append(colors[bidx + 1]) - log.info("data transform output: %s" % str(out)) + LOG.info("data transform output: %s" % str(out)) return out class MapDrawOptions(object): """ Map draw options """ - def __init__(self, loglevel='WARNING', **kwargs): - log.setLevel(loglevel) + def __init__(self, loglevel='INFO', **kwargs): + LOG.setLevel(loglevel) ## General Options self.area_thresh = None # Default area threshold for maps self.resolution = None # Default resolution is intermediate self.xsize = '1.' # (when keep_aspect = False) Default horizontal size [value between 0 and 1] - self.ysize = '1.' # (when keep_aspect = False) Default vertical size [value between 0 and 1] - self.dpi = '200' # Default DPI + self.ysize = '1.' # (when keep_aspect = False) Default vertical size [value between 0 and 1] # Default DPI self.coordsopts = ".3,8,grey" self.coastsopts = ".5,grey" self.countropts = ".3,grey" @@ -658,29 +666,29 @@ class MapDrawOptions(object): self.wind_head_width = 8 # Wind arrowhead width default self.wind_width = 0.008 # Wind arrow width default self.wind_minshaft = 0.25 # Wind arrow minimum shaft default - self.wind_label_xpos = 0.8 # Wind label x position (relative to quiver plot) - self.wind_label_ypos = -0.075 # Wind label y position (relative to quiver plot) - self.wind_label_scale = 20 # Size of the label (represents m/s) + self.wind_label_xpos = 0.8 # Wind label x position (relative to quiver plot) + self.wind_label_ypos = -0.075 # Wind label y position (relative to quiver plot) + self.wind_label_scale = 20 # Size of the label (represents m/s) - ## CONTOUR OPTIONS + # CONTOUR OPTIONS self.contours_int = 1 # Default contour interval - self.contours_color = 'k' # Default contour color + self.contours_color = 'k' # Default contour color self.contours_linewidth = 0.8 # Default line width for contours - self.contours_label = False # if True labels of contours are showed - self.contours_label_fontsize = 13 # Default font size for contour label - self.contours_label_format = '%d' # Default format for contour label + self.contours_label = False # if True labels of contours are showed + self.contours_label_fontsize = 13 # Default font size for contour label + self.contours_label_format = '%d' # Default format for contour label self.contours_exclude_vals = [0] # Default list of values to exclude from drawing # self.nocontourf = False # if True contourf is hidden # self.nomap = False # if True map is hidden # self.nointerp = False # if False use contourf, else pcolor - ## SHAPEFILES + # SHAPEFILES self.shapefiles = [] # Default shapefiles list is empty self.shapef_width_init = 0.75 # Initial line width for shapefiles self.shapef_width_step = 0.25 # Line width difference between two consecutive shapefiles - log.setLevel(loglevel) + LOG.setLevel(loglevel) def has_shape_files(self): return (len(self.shapefiles) > 0) @@ -701,9 +709,12 @@ class MapGenerator(object): self.outdir = kwargs.get('outdir', '.') self.srcgaps = kwargs.get('srcgaps', [0]) self.title = kwargs.get('title', '') - self.fontsize = kwargs.get('fontsize', 12) - self.keep_aspect = kwargs.get('keep_aspect', False) - self.projection = kwargs.get('projection', 'cyl') + self.fontsize = kwargs.get('fontsize', 10) + self.loglevel = kwargs.get('loglevel', 'WARNING') + self.orig_projection = kwargs.get('orig_projection', 'PlateCarree') + self.orig_projection_kwargs = kwargs.get('orig_projection_kwargs', dict()) + self.projection = kwargs.get('projection', 'PlateCarree') + self.projection_kwargs = kwargs.get('projection_kwargs', dict()) self.subsetting = kwargs.get('subsetting', True) self.dimension = kwargs.get('dimension', None) self.continents = kwargs.get('continents', None) @@ -720,24 +731,35 @@ class MapGenerator(object): self.config_dir = kwargs.get('config_dir', \ os.path.join(os.environ['HOME'], '.mapgenerator')) self.config_file = kwargs.get('config_file', 'config.cfg') - self.debug = kwargs.get('debug', False) self.alpha = kwargs.get('alpha', None) self.shapefiles = kwargs.get('shapefiles', None) + self.dpi = kwargs.get('dpi', 200) + + @staticmethod + def _get_default_title(cube): + if cube.attributes.get('plot_name', None): + show_name = cube.attributes['plot_name'] + elif cube.long_name and len(cube.long_name) < 45: + show_name = cube.long_name + elif cube.standard_name and len(cube.standard_name) < 45: + show_name = cube.standard_name + else: + show_name = cube.var_name + return f"{show_name} ({cube.units})" class MapCross(MapGenerator): """ Define common attributes to maps and cross sections """ - def __init__(self, loglevel='WARNING', **kwargs): + def __init__(self, loglevel='INFO', **kwargs): """ Initialize class with MapGenerator attributes plus some others """ - log.setLevel(loglevel) + LOG.setLevel(loglevel) MapGenerator.__init__(self, **kwargs) self.bounds = kwargs.get('bounds', None) self.ticks = kwargs.get('ticks', None) self.colors = kwargs.get('colors', 'jet') -# self.over = kwargs.get('over', None) -# self.under = kwargs.get('under', None) self.bad = kwargs.get('bad', None) + self.N = kwargs.get('N', None) self.lat = kwargs.get('lat', []) self.lon = kwargs.get('lon', []) self.wind = kwargs.get('wind', []) diff --git a/mapgenerator/plotting/plotmap.py b/mapgenerator/plotting/plotmap.py index c443c6c72fbc292150d39171fe29dfa0596565f3..5dc08363a82506662d32a455c3341eb4e48c3753 100644 --- a/mapgenerator/plotting/plotmap.py +++ b/mapgenerator/plotting/plotmap.py @@ -4,26 +4,20 @@ # @license: https://www.gnu.org/licenses/gpl-3.0.html # @author: see AUTHORS file - import matplotlib as mpl -mpl.use('Agg') - import matplotlib.pyplot as plt -from mpl_toolkits.basemap import Basemap +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import cartopy.io.shapereader as shpreader from matplotlib.cbook import delete_masked_points -from matplotlib.image import imread +from matplotlib.path import Path +# from matplotlib.image import imread import numpy as np from datetime import datetime -#from datetime import timedelta -#import time import os.path import subprocess -#import sys from collections import Iterable from PIL import Image -#from mapgenerator.plotting.mapgenerator import MapCross -#from mapgenerator.plotting.mapgenerator import MapDrawOptions -#from mapgenerator.plotting.mapgenerator import MapDataHandler from .definitions import MapData from .definitions import MapCross from .definitions import MapDrawOptions @@ -35,21 +29,20 @@ from .definitions import parse_parameter from .definitions import parse_parameters from .definitions import parse_parameters_list from .definitions import parse_path -from urllib.request import urlopen +from .tools import set_title +from .tools import set_resolution +from .tools import gen_kml +from .tools import gen_anim +# try: +# from urllib2 import urlopen +# except ImportError: +# from urllib.request import urlopen import copy from glob import glob - import logging -#logging.basicConfig(level=logging.WARNING) -log = logging.getLogger(__name__) - - -class PlotSeries(object): - """ Main class for plotting time series """ - - def __init__(self, loglevel='WARNING', **kwargs): - pass +mpl.use('Agg') +LOG = logging.getLogger(__name__) class PlotCross(object): @@ -76,45 +69,70 @@ class PlotMap(MapCross, MapDrawOptions): self.kmz = False self.norm = None self.cmap = None - #self.formats = None - self.map = None #Map - #self.first_image = True # Is is the first image of the batch? - self.map_name = None # name of the figure to map animation colors + self._orig_crs = ccrs.PlateCarree() + self._crs = ccrs.PlateCarree() + self.mgaxis = None + self.mgplot = None + self.stamen_terrain = None + self.zoom_level = None + # self.map = None #Map + self.first_image = False # Is is the first image of the batch? + self.map_name = None # name of the figure to map animation colors # self.drawopts = {'coastlines': # {'linewidth': 0.5, 'color': 'grey', }, # 'countries': # {'linewidth': 0.3, 'color': 'grey', }, # } - log.setLevel(loglevel) def __setattr__(self, key, value): - #log.info("SETATTR: {} - {}".format(key, value)) + LOG.debug("SETATTR: %s - %s", key, value) super(PlotMap, self).__setattr__(key, parse_parameter(key, value)) - def set_resolution(self): - """ 'c','l','i','h','f' """ - # lats = 36, 72, 108, 144, 180 - # lons = 72, 144, 216, 288, 360 - minlat, maxlat = self.lat[0], self.lat[-1] - minlon, maxlon = self.lon[0], self.lon[-1] - lats = maxlat + abs(minlat) - lons = maxlon + abs(minlon) - if lats > 130 and lons > 260: - self.resolution = 'c' - elif lats > 100 and lons > 200: - self.resolution = 'l' - elif lats > 60 and lons > 120: - self.resolution = 'i' - elif lats > 30 and lons > 60: - self.resolution = 'h' - elif lats <= 30 and lons <= 60: - self.resolution = 'f' + def _build_orig_projection(self): + self._orig_crs = getattr(ccrs, self.orig_projection)( + **self.orig_projection_kwargs) + + def _build_projection(self): + if self.background: + import cartopy.io.img_tiles as cimgt + self.stamen_terrain = cimgt.Stamen('terrain-background') + self._crs = self.stamen_terrain.crs + self.projection = 'Stamen' else: - self.resolution = 'c' # default + self._crs = getattr(ccrs, self.projection)( + **self.projection_kwargs) def set_color_map(self): """ Create color map """ - if not isinstance(self.colors, str) and isinstance(self.colors, Iterable): + + extend_opts = { + 'max': { + 'cmap_arg': lambda colors: colors[:-1], + 'cmap_ext': ('set_over'), + 'ext_arg': (lambda colors: colors[-1]), + 'norm_arg': (lambda bounds: bounds+[np.inf], + lambda cmap: cmap.N+1) + }, + 'min': { + 'cmap_arg': lambda colors: colors[1:], + 'cmap_ext': ('set_under'), + 'ext_arg': (lambda colors: colors[0]), + 'norm_arg': (lambda bounds: [-np.inf]+bounds, + lambda cmap: cmap.N+1) + }, + 'both': { + 'cmap_arg': self.colors[1:-1], + 'cmap_ext': ('set_over', 'set_under'), + 'ext_arg': (lambda colors: colors[-1], + lambda colors: colors[0]), + 'norm_arg': (lambda bounds: [-np.inf]+bounds+[np.inf], + lambda cmap: cmap.N + 2) + + } + } + + if not isinstance(self.colors, str) and \ + isinstance(self.colors, Iterable): # adjust number of colors to number of bounds if self.bounds: if self.extend in ('min', 'max'): @@ -123,23 +141,23 @@ class PlotMap(MapCross, MapDrawOptions): elif self.extend == 'both': if len(self.bounds)+1 < len(self.colors): self.colors = self.colors[:len(self.bounds)+1] - if self.extend == 'max': - self.cmap = mpl.colors.ListedColormap(self.colors[:-1]) - self.cmap.set_over(self.colors[-1]) - elif self.extend == 'min': - self.cmap = mpl.colors.ListedColormap(self.colors[1:]) - self.cmap.set_under(self.colors[0]) - elif self.extend == 'both': - self.cmap = mpl.colors.ListedColormap(self.colors[1:-1]) - self.cmap.set_over(self.colors[-1]) - self.cmap.set_under(self.colors[0]) + + if self.extend in extend_opts: + self.cmap = mpl.colors.ListedColormap( + extend_opts[self.extend]['cmap_arg'](self.colors)) + for ext, arg in zip(extend_opts[self.extend]['cmap_ext'], + extend_opts[self.extend]['ext_arg']): + getattr(self.cmap, ext)(arg(self.colors)) else: self.cmap = mpl.colors.ListedColormap(self.colors) custom_cmap = True else: try: - self.cmap = mpl.cm.get_cmap(self.colors) - except: + if self.N: + self.cmap = mpl.cm.get_cmap(self.colors, self.N) + else: + self.cmap = mpl.cm.get_cmap(self.colors) + except ValueError: self.cmap = mpl.cm.get_cmap('jet') custom_cmap = False @@ -147,28 +165,25 @@ class PlotMap(MapCross, MapDrawOptions): self.cmap.set_bad(self.bad) if self.background or self.kml or self.kmz: - try: - self.cmap.colors = [(0, 0, 0, 0),]+[c for c in self.cmap.colors[1:]] - except: # if self.extend in ('min','both'): + if self.extend in ('min', 'both'): self.cmap.set_under(color=(0, 0, 0, 0)) -# else: -# self.cmap.colors = [(0,0,0,0),]+[c for c in self.cmap.colors[1:]] + else: + self.cmap.colors = [(0, 0, 0, 0), ] + [c for c in + self.cmap.colors[1:]] # normalize colormap if self.bounds and self.smooth and not custom_cmap: - if self.extend == 'min': - self.norm = mpl.colors.BoundaryNorm([-np.inf]+self.bounds, self.cmap.N+1) - elif self.extend == 'max': - self.norm = mpl.colors.BoundaryNorm(self.bounds+[np.inf], self.cmap.N+1) - elif self.extend == 'both': - self.norm = mpl.colors.BoundaryNorm([-np.inf]+self.bounds+[np.inf], self.cmap.N+2) + if self.extend in ('min', 'max', 'both'): + self.norm = mpl.colors.BoundaryNorm( + extend_opts[self.extend]['norm_arg'][0](self.bounds), + extend_opts[self.extend]['norm_arg'][1](self.cmap) + ) else: self.norm = mpl.colors.BoundaryNorm(self.bounds, self.cmap.N) elif self.bounds: self.norm = mpl.colors.BoundaryNorm(self.bounds, self.cmap.N) - - def set_color_bar(self, mco, location='left', drawedges=False, cax=None): + def set_color_bar(self, mco, location='right', drawedges=False, cax=None): """ Create color bar """ # xs = self.xsize # nap = .8-((1-xs)*0.4) @@ -177,36 +192,30 @@ class PlotMap(MapCross, MapDrawOptions): mpl.rcParams['axes.linewidth'] = 0.1 mpl.rcParams['axes.formatter.useoffset'] = False if self.ticks: - #log.info("ticks %s" % str(self.ticks)) + LOG.debug("ticks %s" % str(self.ticks)) self.ticks = parse_parameters_list(self.ticks) elif self.bounds: self.ticks = self.bounds - else: - try: - self.ticks = mco.levels - except: - pass - - log.info("***** {} *****".format(self.formats)) - if self.map: - cb = self.map.colorbar(mco, ticks=self.ticks, format=self.formats, - pad="2%", size="3%", extend=self.extend, - drawedges=drawedges) - else: - cb = plt.colorbar(mco, ticks=self.ticks, format=self.formats, - pad=.06, extend=self.extend, - drawedges=drawedges) - - cb.ax.tick_params(labelsize=float(self.coordsopts[1])) - for lin in cb.ax.yaxis.get_ticklines(): + elif hasattr(mco, 'levels'): + self.ticks = mco.levels + + LOG.debug("***** Formats: %s *****", self.formats) + cax, kwargs = mpl.colorbar.make_axes(self.mgaxis, location=location, + pad=0.02, shrink=0.8) + cbar = self.mgplot.colorbar(mco, cax=cax, ax=self.mgaxis, + ticks=self.ticks, format=self.formats, + pad=.06, extend=self.extend, + drawedges=drawedges, **kwargs) + + cbar.ax.tick_params(labelsize=float(self.coordsopts[1])) + for lin in cbar.ax.yaxis.get_ticklines(): lin.set_visible(False) - #self.print_time("colorbar") - ## If required, draw arrow pointing to the specified limits + # If required, draw arrow pointing to the specified limits if self.has_limits(): # The following two values are absolute shrink = 0.75 - upper = 0.863 - 0.45 * (0.8 - shrink) # formula seems to work ok + upper = 0.863 - 0.45 * (0.8 - shrink) # formula seems to work ok lower = 0.137 + 0.45 * (0.8 - shrink) span = upper - lower arrows = self.limits @@ -216,114 +225,39 @@ class PlotMap(MapCross, MapDrawOptions): else: ticks = mco.levels - f = span/float(len(ticks) - 1) + fac = span/float(len(ticks) - 1) for i in arrows: for j in range(0, len(ticks) - 1): - if (i >= ticks[j] and i <= ticks[j+1]): - if (i == ticks[j+1]): - p = j+1 + if i >= ticks[j] and i <= ticks[j+1]: + if i == ticks[j+1]: + idx = j+1 add = 0 else: - p = j - add = (float(i - ticks[j])/float(ticks[j+1] - ticks[j]))*f - pos = f*p + lower + idx = j + add = (float(i - ticks[j])/float(ticks[j+1] - + ticks[j]))*fac + pos = fac*idx + lower plt.arrow(0.88, pos + add, 0.10, 0, length_includes_head=True, head_length=0.10, head_width=0.025, fill=True, color='k') - #plt.colorbar(drawedges=drawedges, cax=cax, ticks=self.bounds) - - def set_title(self, title, sdate, cdate, step, stime): - """ Set the title of the current image to the one provided, substituting the - patterns """ - - title_dic = { - 'syear': sdate.strftime("%Y"), - 'smonth': sdate.strftime("%m"), - 'sMONTH': sdate.strftime("%b"), - 'sday': sdate.strftime("%d"), - 'shour': sdate.strftime("%H"), - 'year': cdate.strftime("%Y"), - 'month': cdate.strftime("%m"), - 'MONTH': cdate.strftime("%b"), - 'day': cdate.strftime("%d"), - 'hour': cdate.strftime("%H"), - 'step': step, - 'simday': "%d" % (int(stime)/24), - 'simhh': stime, - #'simhh': "%02d" % (int(stime)%24) - } - - try: - p_title = str(title) % title_dic #("%sh" % s_time, valid) - except Exception as e: - print("Title error: %s", str(e)) - p_title = title % title_dic #("%sh" % s_time, valid) - - return p_title - - def print_time(self, msg=None): - """ Print time """ - # initial time - if 'start_t' not in globals(): - global start_t, last_t - print("Creating start_t") - start_t = datetime.now() - last_t = start_t - - now_t = datetime.now() - diff_t = now_t - last_t - if msg != None: - print('TIME:', msg, ' done in ', diff_t.seconds, ' s') - last_t = now_t - - def run_command(self, comm_str, fatal=True): - """ Run command """ - print(comm_str) - p = subprocess.Popen(comm_str.split(), stderr=subprocess.PIPE) - output, err = p.communicate() - print(output, '-----', err) -# st, out = commands.getstatusoutput(comm_str) -# if st != 0: -# print "Error: %s" % str(out) -# if (fatal == True): -# sys.exit(1) - - def gen_anim(self, inpattern, outfile): - ''' Generate an animation starting from a set of - images, specified by the inpattern parameter ''' - #print "Animation: ", indir, inpattern, outdir, outfile - # Create animation - self.run_command( - "/usr/bin/convert -delay %s -loop 0 -layers Optimize %s/%s %s/%s" % \ - (self.anim_delay, self.outdir, inpattern, self.outdir, outfile)) - # Remove intermediate files - #self.run_command("rm %s/%s %s" % (indir, inpattern, self.map_name), False) - -# ccc = [] -# ccc.extend(self.colors) -# print "Filtering scatter data for date/hour: ", curr_date.strftime("%Y%m%d"), s_time24 -# csv_lat, csv_lon, csv_val = csv_handler.filter(curr_date.strftime("%Y%m%d"), s_time24) -# #csv_col = self.dataTransform.setColorsFromBuckets(ccc, self.bounds, csv_con) -# if(len(csv_lat) > 0): -# scatter_data = [csv_lon, csv_lat, 40, csv_col] -# print "**************** SCATTER_DATA", scatter_data -# print "Filtered %d records for scatterplot" % len(scatter_data) - + # plt.colorbar(drawedges=drawedges, cax=cax, ticks=self.bounds) def init_map(self, grid): - """ Initialize a Basemap map. Initialization should be performed only once + """ Initialize a map. Initialization should be performed only once the beginning of a serie of images """ glon, glat = grid - #print("***",self.lon, self.lat,"***") - if not self.lat: # or len(self.lat) not in (2,3): + # print("***",self.lon, self.lat,"***") + if not self.lat: # or len(self.lat) not in (2,3): self.lat = "%s-%s" % (str(round(glat.min(), 1)).replace('-', 'm'), str(round(glat.max(), 1)).replace('-', 'm')) - if not self.lon: # or len(self.lon) not in (2,3): + if not self.lon: # or len(self.lon) not in (2,3): self.lon = "%s-%s" % (str(round(glon.min(), 1)).replace('-', 'm'), str(round(glon.max(), 1)).replace('-', 'm')) - #print(self.lon, self.lat) + self._build_orig_projection() + self._build_projection() + # print(self.lon, self.lat) # # Fix the printout of tick values to avoid .0 decimals in integers # strs = [] # if self.bounds: @@ -336,26 +270,30 @@ class PlotMap(MapCross, MapDrawOptions): # Set the colormap self.set_color_map() if not self.resolution: - self.set_resolution() + self.resolution, self.zoom_level = \ + set_resolution(self.lon, self.lat) + else: + _, self.zoom_level = set_resolution(self.lon, self.lat) self.first_image = True - - def gen_image_map(self, fig_name, grid, data, img_title, cur_scatter_data=None, **kwargs): + def gen_image_map(self, fig_name, grid, data, img_title, + cur_scatter_data=None, **kwargs): """ Generate image map """ # overwrite option if os.path.exists(fig_name + ".png") and not self.overwrite: - print(fig_name, " already exists.") + LOG.info(fig_name, " already exists.") if self.subplot is None: plt.clf() return fig_name - map_data = data.getMapData() + map_data = data.map_data # FIXME - scatter_data = data.getScatterData() or cur_scatter_data - if self.subplot is None: - plt.clf() + scatter_data = data.scatter_data or cur_scatter_data + # if self.subplot is None: + # self.mgplot.clear() + # params = { # 'font.size': 14, # 'text.fontsize': 28, @@ -365,62 +303,76 @@ class PlotMap(MapCross, MapDrawOptions): # # mpl.rcParams.update(params) - if not self.nomap: - if self.projection != 'cyl': - lon_0 = (self.lon[-1]-abs(self.lon[0]))/2 - lat_0 = (self.lat[-1]-abs(self.lat[0]))/2 - #print("-----------", lon_0, lat_0, "-----------") - else: - lon_0, lat_0 = None, None - if self.subplot: - self.map = Basemap( - ax=self.mgaxis, - projection=self.projection, resolution=self.resolution, - llcrnrlon=self.lon[0], llcrnrlat=self.lat[0], - urcrnrlon=self.lon[-1], urcrnrlat=self.lat[-1], - lon_0=lon_0, lat_0=lat_0, - fix_aspect=self.keep_aspect, - area_thresh=self.area_thresh, - boundinglat=self.boundinglat, - ) - else: - self.map = Basemap( - projection=self.projection, resolution=self.resolution, - llcrnrlon=self.lon[0], llcrnrlat=self.lat[0], - urcrnrlon=self.lon[-1], urcrnrlat=self.lat[-1], - lon_0=lon_0, lat_0=lat_0, - fix_aspect=self.keep_aspect, - area_thresh=self.area_thresh, - boundinglat=self.boundinglat, - ) + # Draw filled contour + self.mgaxis.set_aspect(self.xsize/self.ysize) + + if self.projection.lower().endswith('polarstereo'): + self.mgaxis.set_extent([-180, 179.9999999999999, + self.lat[0], self.lat[-1]], + self._orig_crs) + # Compute a circle in axes coordinates, which we can use as a + # boundary for the map. We can pan/zoom as much as we like - the + # boundary will be permanently circular. + theta = np.linspace(0, 2*np.pi, 100) + center, radius = [0.5, 0.5], 0.5 + verts = np.vstack([np.sin(theta), np.cos(theta)]).T + circle = Path(verts * radius + center) + self.mgaxis.set_boundary(circle, transform=self.mgaxis.transAxes) + else: + self.mgaxis.set_extent([self.lon[0], self.lon[-1], + self.lat[0], self.lat[-1]], + self._orig_crs) + +# if self.subplot: +# self.map = Basemap( +# ax=self.mgaxis, +# projection=self.projection, resolution=self.resolution, +# llcrnrlon=self.lon[0], llcrnrlat=self.lat[0], +# urcrnrlon=self.lon[-1], urcrnrlat=self.lat[-1], +# lon_0=lon_0, lat_0=lat_0, +# fix_aspect=self.keep_aspect, +# area_thresh=self.area_thresh, +# boundinglat=self.boundinglat, +# ) +# else: +# self.map = Basemap( +# projection=self.projection, resolution=self.resolution, +# llcrnrlon=self.lon[0], llcrnrlat=self.lat[0], +# urcrnrlon=self.lon[-1], urcrnrlat=self.lat[-1], +# lon_0=lon_0, lat_0=lat_0, +# fix_aspect=self.keep_aspect, +# area_thresh=self.area_thresh, +# boundinglat=self.boundinglat, +# ) glon, glat = grid - #log.info("GLON: %s, GLAT: %s" % (str(glon.shape), str(glat.shape))) + LOG.debug("0. GLON: %s, GLAT: %s", str(glon.shape), str(glat.shape)) + if len(glon.shape) == 1 and len(glat.shape) == 1: glon, glat = np.meshgrid(glon, glat) - #print(glon, glat) - log.info("1. GLON: %s, GLAT: %s" % (str(glon.shape), str(glat.shape))) + LOG.debug("1. GLON: %s, GLAT: %s", str(glon.shape), str(glat.shape)) # if curvilinear grid do interpolation, return already gridded coords if not is_grid_regular(glon, glat): map_data[0], glon, glat = do_interpolation(map_data[0], glon, glat) - log.info("2. GLON: %s, GLAT: %s" % (str(glon.shape), str(glat.shape))) + LOG.info("2. GLON: %s, GLAT: %s", str(glon.shape), str(glat.shape)) - if not self.nomap: - x, y = self.map(*(glon, glat)) - else: - x, y = glon, glat +# if not self.nomap: +# x, y = self.map(*(glon, glat)) +# else: + xloc, yloc = glon, glat - log.info("3. GLON: %s, GLAT: %s" % (str(x.shape), str(y.shape))) + LOG.info("3. GLON: %s, GLAT: %s", str(xloc.shape), str(yloc.shape)) # print "map_data", map_data[0].shape - #self.print_time("meshgrid") - #log.info("X: %s, Y: %s" % (str(x.shape), str(y.shape))) + # self.print_time("meshgrid") + # LOG.info("X: %s, Y: %s" % (str(x.shape), str(y.shape))) # # nomap option -# if self.nomap and os.path.exists("%s-%s/%s.png" % (run_date, varName, fig_name)) and not self.overwrite: +# if self.nomap and os.path.exists("%s-%s/%s.png" % (run_date, var_name, +# fig_name)) and not self.overwrite: # print fig_name, " already exists." # plt.clf() # return fig_name @@ -428,299 +380,308 @@ class PlotMap(MapCross, MapDrawOptions): if self.nomap: self.mgplot.frameon = False - # Draw filled contour - if not self.keep_aspect: -# log.info("Not keeping original aspect") -# log.info("XSIZE: %s" % self.xsize) -# log.info("YSIZE: %s" % self.ysize) - ysize_norm = self.ysize*0.8 - xsize_norm = self.xsize*0.8 - rows = img_title.count("\n") - blines = max(rows - 2, 0) - param = self.fontsize/400 - ystart_base_norm = 0.1 - param*(blines)#img_title - xstart_base_norm = 0.1 - ystart_norm = (0.8 - ysize_norm) / 2 + ystart_base_norm - xstart_norm = (0.8 - xsize_norm) / 2 + xstart_base_norm - #ax_map = self.mgaxis.axes([xstart_norm, ystart_norm, xsize_norm, ysize_norm]) -# log.info("Ynorm: %s, Xnorm: %s, Ystart: %s, Xstart: %s" % (ysize_norm, xsize_norm, ystart_norm, xstart_norm)) - self.mgplot.add_axes([xstart_norm, ystart_norm, xsize_norm, ysize_norm]) - #print "Switching AXES for MAP:", ax_map - - if self.alpha != None and self.bounds != None: - if self.extend in ('min', 'both'): - map_data[0] = np.ma.masked_where(map_data[0] < self.bounds[0], map_data[0]) - else: - map_data[0] = np.ma.masked_where((map_data[0] >= - self.bounds[0]) & - (map_data[0] < - self.bounds[1]), - map_data[0]) +# if self.alpha is not None and self.bounds is not None: +# if self.extend in ('min', 'both'): +# map_data[0] = np.ma.masked_where(map_data[0] < +# self.bounds[0], map_data[0]) +# else: +# map_data[0] = np.ma.masked_where((map_data[0] >= +# self.bounds[0]) & +# (map_data[0] < +# self.bounds[1]), +# map_data[0]) mco = None - log.info(type(map_data[0])) + LOG.info(type(map_data[0])) map_data[0] = np.ma.filled(map_data[0], np.nan) - log.info(type(map_data[0])) - - if not self.nocontourf and not self.nomap and self.smooth: - log.info("X: %s, Y: %s, DATA: %s" % (x.shape, y.shape, map_data[0])) - mco = self.map.contourf(x, y, - map_data[0], - cmap=self.cmap, - norm=self.norm, - levels=self.bounds, - extend=self.extend, - horizontalalignment='center', - alpha=self.alpha, - antialiased=True) - - elif not self.nocontourf and not self.nomap and not self.smooth: - log.info("X: %s, Y: %s, DATA: %s" % (x.shape, y.shape, map_data[0])) - mco = self.map.pcolormesh(x, y, - map_data[0], - cmap=self.cmap, - norm=self.norm, - alpha=self.alpha, - antialiased=True) - - elif not self.nocontourf and self.nomap and self.smooth: - mco = plt.contourf(x, y, - map_data[0], - cmap=self.cmap, - norm=self.norm, - levels=self.bounds, - extend=self.extend, - horizontalalignment='center', - alpha=self.alpha, - antialiased=True) - - plt.axis([self.lon[0], self.lon[-1], self.lat[0], self.lat[-1]]) - - elif not self.nocontourf and self.nomap and not self.smooth: - mco = plt.pcolormesh(x, y, - map_data[0], - cmap=self.cmap, - norm=self.norm, - alpha=self.alpha, - antialiased=True) - - plt.axis([self.lon[0], self.lon[-1], self.lat[0], self.lat[-1]]) + LOG.info(type(map_data[0])) + +# if not self.nocontourf and not self.nomap and self.smooth: +# LOG.info("X: %s, Y: %s, DATA: %s" % (x.shape, y.shape, +# map_data[0])) +# mco = self.map.contourf(x, y, +# map_data[0], +# cmap=self.cmap, +# norm=self.norm, +# levels=self.bounds, +# extend=self.extend, +# horizontalalignment='center', +# alpha=self.alpha, +# antialiased=True) +# +# elif not self.nocontourf and not self.nomap and not self.smooth: +# LOG.info("X: %s, Y: %s, DATA: %s" % (x.shape, y.shape, +# map_data[0])) +# mco = self.map.pcolormesh(x, y, +# map_data[0], +# cmap=self.cmap, +# norm=self.norm, +# alpha=self.alpha, +# antialiased=True) + + if not self.nocontourf and self.smooth: + mco = self.mgaxis.contourf(xloc, yloc, + map_data[0], + cmap=self.cmap, + norm=self.norm, + levels=self.bounds, + extend=self.extend, + alpha=self.alpha, + transform=self._orig_crs, + antialiased=True) + + elif not self.nocontourf and not self.smooth: + mco = self.mgaxis.pcolormesh(xloc, yloc, + map_data[0], + cmap=self.cmap, + norm=self.norm, + alpha=self.alpha, + transform=self._orig_crs, + antialiased=True) # Draw shape files if self.has_shape_files(): line_w = float(self.countropts[0]) for shapef in self.shapefiles: - log.info("Processing shape file: {} with line width: {}".format(shapef, line_w)) - self.map.readshapefile(shapef, "%s" % os.path.basename(shapef), linewidth=line_w, color=self.countropts[1]) - line_w = max(self.shapef_width_step, line_w - self.shapef_width_step) - #self.print_time("contourf") - - ## FIXME Modify to use scatter inside DATA - ## if DATA.hasScatterData()... etc... + LOG.info("Processing shape file: %s with line width: %s", + shapef, line_w) + adm1_shapes = list(shpreader.Reader(shapef).geometries()) + self.mgaxis.add_geometries(adm1_shapes, self._crs, + linewidth=line_w, + edgecolor=self.countropts[1]) + # self.map.readshapefile(shapef, "%s" % + # os.path.basename(shapef), linewidth=line_w, + # color=self.countropts[1]) + line_w = max(self.shapef_width_step, line_w - + self.shapef_width_step) + + # FIXME Modify to use scatter inside DATA + # if DATA.hasScatterData()... etc... if scatter_data is not None: - log.info("Plotting scatter data: {} of keys {}".format(str(scatter_data), str(scatter_data.keys()))) + LOG.info("Plotting scatter data: %s of keys %s", + str(scatter_data), str(scatter_data.keys())) if mco and self.bounds: - self.map.scatter(scatter_data['lon'].tolist(), - scatter_data['lat'].tolist(), - s=scatter_data['size'].tolist(), - c=scatter_data['color'].tolist(), - marker='o', - linewidth=0.3, - vmin=self.bounds[0], - vmax=self.bounds[-1], - cmap=self.cmap, - norm=self.norm, - zorder=10) + self.mgaxis.scatter(scatter_data['lon'].tolist(), + scatter_data['lat'].tolist(), + s=scatter_data['size'].tolist(), + c=scatter_data['color'].tolist(), + marker='o', + linewidth=0.3, + vmin=self.bounds[0], + vmax=self.bounds[-1], + cmap=self.cmap, + norm=self.norm, + transform=self._orig_crs, + zorder=10) elif mco: - self.map.scatter(scatter_data['lon'].tolist(), - scatter_data['lat'].tolist(), - s=scatter_data['size'].tolist(), - c=scatter_data['color'].tolist(), - marker='o', - linewidth=0.3, - cmap=self.cmap, - norm=self.norm, - zorder=10) + self.mgaxis.scatter(scatter_data['lon'].tolist(), + scatter_data['lat'].tolist(), + s=scatter_data['size'].tolist(), + c=scatter_data['color'].tolist(), + marker='o', + linewidth=0.3, + cmap=self.cmap, + norm=self.norm, + transform=self._orig_crs, + zorder=10) elif self.bounds: - mco = self.map.scatter(scatter_data['lon'].tolist(), - scatter_data['lat'].tolist(), - s=scatter_data['size'].tolist(), - c=scatter_data['color'].tolist(), - marker='o', - linewidth=0.3, - vmin=self.bounds[0], - vmax=self.bounds[-1], - cmap=self.cmap, - norm=self.norm, - zorder=10) + mco = self.mgaxis.scatter(scatter_data['lon'].tolist(), + scatter_data['lat'].tolist(), + s=scatter_data['size'].tolist(), + c=scatter_data['color'].tolist(), + marker='o', + linewidth=0.3, + vmin=self.bounds[0], + vmax=self.bounds[-1], + cmap=self.cmap, + norm=self.norm, + transform=self._orig_crs, + zorder=10) else: - mco = self.map.scatter(scatter_data['lon'].tolist(), - scatter_data['lat'].tolist(), - s=scatter_data['size'].tolist(), - c=scatter_data['color'].tolist(), - marker='o', - linewidth=0.3, - cmap=self.cmap, - norm=self.norm, - zorder=10) - - if data.hasWindData(): - winds = data.getWindData() - X, Y, U, V = delete_masked_points(x.ravel(), y.ravel(), winds['u'].ravel(), winds['v'].ravel()) - #print("Wind scale is", self.wind_units, self.wind_scale) - if winds.has_key('barbs'): - self.map.barbs(X, Y, U, V, #units=self.wind_units, - #headlength=self.wind_head_length, - #headwidth=self.wind_head_width, - #width=self.wind_width, - #minshaft=self.wind_minshaft, - #scale=self.wind_scale, - color='k') - #self.print_time("barbs") + mco = self.mgaxis.scatter(scatter_data['lon'].tolist(), + scatter_data['lat'].tolist(), + s=scatter_data['size'].tolist(), + c=scatter_data['color'].tolist(), + marker='o', + linewidth=0.3, + cmap=self.cmap, + norm=self.norm, + transform=self._orig_crs, + zorder=10) + + if data.wind_data: + winds = data.wind_data + x_vec, y_vec, u_vec, v_vec = delete_masked_points( + xloc.ravel(), + yloc.ravel(), + winds['u'].ravel(), + winds['v'].ravel() + ) + + if 'barbs' in winds: + self.mgaxis.barbs(x_vec, y_vec, u_vec, v_vec, + units=self.wind_units, + headlength=self.wind_head_length, + headwidth=self.wind_head_width, + width=self.wind_width, + minshaft=self.wind_minshaft, + scale=self.wind_scale, + transform=self._orig_crs, + color='k') else: - Q = self.map.quiver(X, Y, U, V, units=self.wind_units, - headlength=self.wind_head_length, - headwidth=self.wind_head_width, - width=self.wind_width, - minshaft=self.wind_minshaft, - scale=self.wind_scale, - color='gray') + quiv = self.mgaxis.quiver(x_vec, y_vec, u_vec, v_vec, + units=self.wind_units, + headlength=self.wind_head_length, + headwidth=self.wind_head_width, + width=self.wind_width, + minshaft=self.wind_minshaft, + scale=self.wind_scale, + transform=self._orig_crs, + color='gray') # Draw the key - plt.quiverkey(Q, + plt.quiverkey(quiv, self.wind_label_xpos, self.wind_label_ypos, self.wind_label_scale, label='%s m/s' % self.wind_label_scale, coordinates='axes', labelpos='S', labelsep=0.05) - #self.print_time("quivers") - if data.hasContourData(): + if data.contour_data: interval = self.contours_int exclude = self.contours_exclude_vals - cLowBound = -99999 - cUppBound = 99999 - cdata = data.getContourData() - #print(cdata, type(cdata)) + clow_bound = -99999 + cupp_bound = 99999 + cdata = data.contour_data try: - cMin = min(filter (lambda a: a > cLowBound, cdata.ravel())) - adjcMin = int(cMin - (cMin % interval) - interval * 2) + # cmin = min(filter (lambda a: a > clow_bound, cdata.ravel())) + cmin = cdata.where(cdata > clow_bound).min() + adjcmin = int(cmin - (cmin % interval) - interval * 2) except ValueError: - cMin =0 + cmin = 0 try: - cMax = max(filter (lambda a: a < cUppBound, cdata.ravel())) - adjcMax = int(cMax - (cMax % interval) + interval * 2) + # cmax = max(filter (lambda a: a < cupp_bound, cdata.ravel())) + cmax = cdata.where(cdata < cupp_bound).max() + adjcmax = int(cmax - (cmax % interval) + interval * 2) except ValueError: - cMax = 0 - lvls = np.arange(adjcMin, adjcMax, interval) - for ex in exclude: - lvls = filter (lambda ex: ex != 0, lvls) - if (len(lvls) > 0) and map_data != None: - mco = self.map.contourf(x, y, - map_data[0], - cmap=self.cmap, - norm=self.norm, - levels=self.bounds, - extend=self.extend, - horizontalalignment='center', - alpha=self.alpha) - - #print "-----------------", map_data, cdata, "--------------------" - if map_data != None and (map_data[0] == cdata).all(): - print(":::::::::::::::::::::::: SAME !!! :::::::::::::::::::") - CS = plt.contour(x, y, - cdata, - levels=self.bounds, - colors=self.contours_color, - linewidths=self.contours_linewidth, - alpha=self.alpha) + cmax = 0 + lvls = np.arange(adjcmin, adjcmax, interval) + for exc in exclude: + lvls = [exc for exc in lvls if exc != 0] + if (len(lvls) > 0) and map_data is not None: + mco = plt.contourf(xloc, yloc, + map_data[0], + cmap=self.cmap, + norm=self.norm, + levels=self.bounds, + extend=self.extend, + horizontalalignment='center', + alpha=self.alpha) + + if map_data is not None and (map_data[0] == cdata).all(): + LOG.debug(":::::::::::::::::::::: SAME !!! ::::::::::::::::") + mcs = plt.contour(xloc, yloc, + cdata, + levels=self.bounds, + colors=self.contours_color, + linewidths=self.contours_linewidth, + alpha=self.alpha) else: - print(":::::::::::::::::::::::: DIFFERENT !!! :::::::::::::::::::") - print("MIN:", cMin, "MAX:", cMax) - CS = plt.contour(x, y, - cdata, - levels=lvls, - colors=self.contours_color, - linewidths=self.contours_linewidth, - alpha=self.alpha) + LOG.debug(":::::::::::::::::::: DIFFERENT !!! :::::::::::::") + LOG.debug("MIN: %s - MAX: %s", cmin, cmax) + mcs = plt.contour(xloc, yloc, + cdata, + levels=lvls, + colors=self.contours_color, + linewidths=self.contours_linewidth, + alpha=self.alpha) if self.contours_label: - self.mgplot.clabel(CS, + self.mgplot.clabel(mcs, inline=1, fontsize=self.contours_label_fontsize, - #backgroundcolor='r', + # backgroundcolor='r', fmt=self.contours_label_format) -# else: -# print "Not drawing contours since levels has ", len(lvls), " elements" - #coords normalization + # coords normalization # lat_offset = abs(self.lat[0]) % self.lat[2] # lon_offset = abs(self.lon[0]) % self.lon[2] if not self.nomap and not self.kml and not self.kmz: -# self.map.drawparallels(np.arange(self.lat[0]+lat_offset, self.lat[1], self.lat[2]),labels=[1,0,0,0],linewidth=self.coords_linewidth, fontsize=self.coords_fontsize) -# self.map.drawmeridians(np.arange(self.lon[0]+lon_offset, self.lon[1], self.lon[2]),labels=[0,0,0,1],linewidth=self.coords_linewidth, fontsize=self.coords_fontsize) - #self.print_time("coords") - self.map.drawparallels(self.lat, labels=[1,0,0,0],linewidth=float(self.coordsopts[0]), fontsize=float(self.coordsopts[1]), color=str(self.coordsopts[2])) - self.map.drawmeridians(self.lon, labels=[0,0,0,1],linewidth=float(self.coordsopts[0]), fontsize=float(self.coordsopts[1]), color=str(self.coordsopts[2])) - self.map.drawcountries(linewidth=float(self.countropts[0]), color=str(self.countropts[1])) - self.map.drawcoastlines(linewidth=float(self.coastsopts[0]), color=str(self.coastsopts[1])) + if self.continents: - self.map.fillcontinents(color=self.continents, zorder=0) - -# drawopts = self.drawopts -# if drawopts: -# # Drawing countries -# if 'countries' in drawopts: -# self.map.drawcountries(**drawopts['countries']) -# #self.print_time("countries") -# # Drawing coastlines (land/sea, lakes, etc...) -# if 'coastlines' in drawopts: -# self.map.drawcoastlines(**drawopts['coastlines']) -# #self.print_time("coastlines") -## # Drawing states (FIXME does it work outside the U.S.?) -## if 'states' in drawopts: -## self.map.drawstates(**drawopts['states']) -## #self.print_time("states") + self.mgaxis.add_feature(cfeature.LAND, color=self.continents, + zorder=10) + + self.mgaxis.coastlines(resolution=self.resolution, + linewidth=float(self.coastsopts[0]), + color=str(self.coastsopts[1]), + zorder=15) + + self.mgaxis.add_feature(cfeature.BORDERS + .with_scale(self.resolution), + linewidth=float(self.countropts[0]), + edgecolor=str(self.countropts[1]), + zorder=15) + + draw_labels = bool(self.projection == 'PlateCarree') + + grl = self.mgaxis.gridlines(xlocs=self.lon, ylocs=self.lat, + crs=self._crs, + draw_labels=draw_labels, + linestyle='--', + linewidth=float(self.coordsopts[0]), + color=str(self.coordsopts[2]), + zorder=20) + + if draw_labels: + grl.xlabels_top = False + grl.ylabels_right = False + grl.xlabel_style = {'size': float(self.coordsopts[1])} + grl.ylabel_style = grl.xlabel_style # Change axes for colorbar if self.colorbar and not self.kml and not self.kmz and \ (not self.nocontourf or scatter_data is not None): self.set_color_bar(mco) - #m.bluemarble() - if self.background == 'shadedrelief': - self.map.shadedrelief() - if self.background == 'bluemarble': - self.map.bluemarble() - if self.background == 'etopo': - self.map.etopo() - if self.background == 'GIS': - h0 = float(self.lat[-1] - self.lat[0]) - w0 = float(self.lon[-1] - self.lon[0]) - ff = h0/w0 # form factor - height = 1024*ff # height - size = "%d,%d" % (1024, height) - basemap_url = "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/export?\ -bbox=%s,%s,%s,%s&\ -bboxSR=4326&\ -imageSR=4326&\ -size=%s&\ -dpi=%d&\ -format=png24&\ -f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi) - self.map.imshow(imread(urlopen(basemap_url)), origin='upper') - - #logo + if self.background: + self.mgaxis.add_image(self.stamen_terrain, self.zoom_level) +# if self.background == 'shadedrelief': +# #m.bluemarble() +# if self.background == 'shadedrelief': +# self.map.shadedrelief() +# if self.background == 'bluemarble': +# self.map.bluemarble() +# if self.background == 'etopo': +# self.map.etopo() +# if self.background == 'GIS': +# h0 = float(self.lat[-1] - self.lat[0]) +# w0 = float(self.lon[-1] - self.lon[0]) +# ff = h0/w0 # form factor +# height = 1024*ff # height +# size = "%d,%d" % (1024, height) +# basemap_url = +# "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/export?\ +# bbox=%s,%s,%s,%s&\ +# bboxSR=4326&\ +# imageSR=4326&\ +# size=%s&\ +# dpi=%d&\ +# format=png24&\ +# f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, +# self.dpi) self.map.imshow(imread(urlopen(basemap_url)), origin='upper') + + # logo if self.logo: - im = Image.open(self.logo[0]) + img = Image.open(self.logo[0]) # height = im.size[1] # width = im.size[0] # We need a float array between 0-1, rather than # a uint8 array between 0-255 - nim = np.array(im).astype(np.float) / 255 + nim = np.array(img).astype(np.float) / 255 # With newer (1.0) versions of matplotlib, you can # use the "zorder" kwarg to make the image overlay @@ -729,74 +690,54 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi # SAVEFIG fullname = "{}.{}".format(fig_name, self.filefmt) - log.warning("printing {}".format(fullname)) + LOG.info("printing %s", fullname) if self.kml or self.kmz: - #self.mgplot.axes(frameon=0) + # self.mgplot.axes(frameon=0) if self.save: - self.mgplot.savefig(fullname, - bbox_inches='tight', - frameon=0, - pad_inches=0, dpi=self.dpi, transparent=True) + self.mgplot.savefig(fullname, bbox_inches='tight', + frameon=0, pad_inches=0, dpi=self.dpi, + transparent=True) else: - #a = plt.axes([0.0, 0.0, 1.0, 1.0]) - #plt.axis('off') - plt.setp(self.mgaxis.spines.values(), linewidth=.1) - xpos = self.xsize/2 - if not self.keep_aspect: - ypos = ystart_norm + ysize_norm + 0.05 - else: - ypos = self.ysize*0.95 + self.mgaxis.set_title(img_title, fontsize=self.fontsize, zorder=0) - if self.subplot is None: - self.mgplot.text( - xpos, - ypos, - img_title, - horizontalalignment='center', - verticalalignment='top', - fontsize=self.fontsize, - zorder=0, - ) - else: - self.mgaxis.set_title(img_title, fontsize=self.fontsize, zorder=0,) if self.save: - self.mgplot.savefig(fullname, bbox_inches='tight', pad_inches=.2, dpi=self.dpi) - return fig_name + self.mgplot.savefig(fullname, bbox_inches='tight', + pad_inches=.2, dpi=self.dpi) + return fig_name def compare_maps(self, map_names, outdir, date, steps, tpl=''): """ Generate maps comparison """ os.chdir(outdir) - #print "MAPNAMES", map_names - #for map_name in map_names: + # print "MAPNAMES", map_names + # for map_name in map_names: # newMapNames = map_names[0] # items = {} #[i[0] for i in map_names] # for i in range(0, len(map_names)-1): -# items, newMapNames, newMapNameTpl = self.twoMapsCompare([newMapNames, map_names[i+1]], items) +# items, newMapNames, newMapNameTpl = +# self.twoMapsCompare([newMapNames, map_names[i+1]], items) - mn = np.array(map_names) - ks = mn.T + mnames = np.array(map_names) + k_s = mnames.T idx = 0 - for item in ks: -# if idx >= len(newMapNames): -# continue + for item in k_s: + # if idx >= len(newMapNames): + # continue imgs = ' '.join([str(i) + '.gif' for i in item]) if tpl: - newimg = tpl % {'date':date, 'step':steps[idx]} + newimg = tpl % {'date': date, 'step': steps[idx]} else: newimg = steps[idx] - #comm = "montage %s -tile 2x -geometry 800x600+1+1 %s.gif" % (imgs, newMapNames[idx]) comm = "montage %s -tile 2x -geometry +1+1 %s.gif" % (imgs, newimg) - #print "COMPARE COMMAND:", comm - st, out = subprocess.getstatusoutput(comm) + _, _ = subprocess.getstatusoutput(comm) idx += 1 - #montage -tile 2x -geometry 800x600+1+1 + # montage -tile 2x -geometry 800x600+1+1 - #rename from 00-06-... to 00-01 ... + # rename from 00-06-... to 00-01 ... # st, out = commands.getstatusoutput("ls %s*" % newMapNameTpl) # if st != 0: # print "Error: %s" % str(out) @@ -805,8 +746,10 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi # #print names # for num in range(0, len(names)): # snum = "%02d" % num -# #print "mv %s %s.gif" % (names[num], newMapNameTpl.replace("??", snum)) -# st2, out2 = commands.getstatusoutput("mv -f %s %s.gif" % (names[num], newMapNameTpl.replace("??", snum))) +# #print "mv %s %s.gif" % (names[num], +# newMapNameTpl.replace("??", snum)) st2, out2 = +# commands.getstatusoutput("mv -f %s %s.gif" % (names[num], +# newMapNameTpl.replace("??", snum))) # if st2 != 0: # print "Error: %s" % str(out2) @@ -814,153 +757,61 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi nimg0 = "??" nimg1 = "loop" if tpl: - nimg0 = tpl % {'date':date, 'step':nimg0} - nimg1 = tpl % {'date':date, 'step':nimg1} - st, out = subprocess.getstatusoutput("convert -delay 75 -loop 0 %s.gif %s.gif" % (nimg0, nimg1)) - if st != 0: - print("Error: %s" % str(out)) - - - def genKML(self, fig_names, srcvars, dims, online=True): - """ """ - from lxml import etree -# from pykml.parser import Schema - from pykml.factory import KML_ElementMaker as KML -# from datetime import datetime, timedelta - import zipfile - - lon = self.lon - lat = self.lat - - print(fig_names) - #kmlName = '-'.join(os.path.basename(fig_names[0]).split('-')[:-1]) - kmlName = os.path.basename(fig_names[0])[:-3] - varName = srcvars[0] #fig_names[0].split('-')[-2] - run_date0 = datetime.strptime(dims[0], "%HZ%d%b%Y") #.strftime("%Y%m%d%H") - run_date1 = datetime.strptime(dims[1], "%HZ%d%b%Y") #.strftime("%Y%m%d%H") - begindate = run_date0.strftime("%Y%m%d%H") - dirName = os.path.join(self.outdir, begindate + '-' + varName).replace('./','') - interval = run_date1 - run_date0 - datefmt = "%Y-%m-%dT%H:00:00Z" - - doc = KML.kml( - KML.Folder( - KML.name(kmlName), - KML.LookAt( - KML.longitude((lon[1]+lon[-1])/2), - KML.latitude((lat[1]+lat[-1])/2), - KML.range(8500000.0), - KML.tilt(0), - KML.heading(0), - ), - KML.ScreenOverlay( - KML.name("Legend"), - KML.open(1), - KML.Icon( - KML.href("%s/%s-colorbar.png" % (dirName, begindate)), - ), - KML.overlayXY(x="0",y="1",xunits="fraction",yunits="fraction"), - KML.screenXY(x=".01",y="0.55",xunits="fraction",yunits="fraction"), - KML.rotationXY(x="0",y="0",xunits="fraction",yunits="fraction"), - KML.size(x="0",y="0.5",xunits="fraction",yunits="fraction"), - id="%s-ScreenOverlay" %(begindate), - ), - ) - ) + nimg0 = tpl % {'date': date, 'step': nimg0} + nimg1 = tpl % {'date': date, 'step': nimg1} + status, out = subprocess.getstatusoutput( + "convert -delay 75 -loop 0 %s.gif %s.gif" % (nimg0, nimg1)) + if status != 0: + LOG.error("Error: %s", str(out)) - try: - os.mkdir(dirName) - except: - pass - - if not online: - zf = zipfile.ZipFile("%s.kmz" % kmlName, 'w') - - for fig_name, dat in zip(fig_names, sorted(dims.keys())): - dt = dims[dat] - fig_name = fig_name + ".png" - figPath = os.path.join(dirName,os.path.basename(fig_name)) - startdate = datetime.strptime(dt, "%HZ%d%b%Y") - begdate = startdate.strftime(datefmt) - enddate = (startdate + interval).strftime(datefmt) - starth = startdate.hour - doc.Folder.append(KML.GroundOverlay( - KML.name("%02d:00:00Z" % starth), - KML.TimeSpan( - KML.begin(begdate), - KML.end(enddate), - ), - KML.Icon( - KML.href(figPath), - KML.viewBoundScale(1.0), - ), - KML.altitude(0.0), - KML.altitudeMode("relativeToGround"), - KML.LatLonBox( - KML.south(lat[0]), - KML.north(lat[-1]), - KML.west(lon[0]), - KML.east(lon[-1]), - KML.rotation(0.0), - ), - )) - - if os.path.exists(fig_name): - os.rename(fig_name, figPath) - zf.write(figPath) - - outf = file("%s.kml" % kmlName, 'w') - outf.write('\n') - outf.write(etree.tostring(doc, pretty_print=True)) - outf.close() - - if not online: - zf.write("%s.kml" % kmlName) - zf.write("%s/%s-colorbar.png" % (dirName, begindate)) - zf.close() - - for img in os.listdir(dirName): - os.remove("%s/%s" % (dirName, img)) - os.rmdir(dirName) - os.remove("%s.kml" % kmlName) - - - def aplot(self, x, y, map_data=None, wind_data=None, contour_data=None, scatter_data=None, **kwargs): - """ Plot one or more maps directly from numerical arrays """ + def init_environ(self, **kwargs): + """ initialize environment with all arguments """ - # store options to be restored at the end - localvars = copy.deepcopy(vars(self)) # update from config file if ('config' in kwargs) and ('section' in kwargs): self.load_conf(kwargs['section'], kwargs['config']) + # update from command line vars(self).update(parse_parameters(kwargs)) # arguments are only local + LOG.setLevel(self.loglevel) + # parse wind arguments if self.windopts and (len(self.windopts) == 4): self.wind_scale = int(self.windopts[3]) - log.info("SHP {}".format(str(self.shapefiles))) + LOG.info("SHP %s", str(self.shapefiles)) if self.shapefiles: - self.shapefiles = [f.replace(".shp", "") for f in glob(self.shapefiles)] - log.info("SHP {}".format(str(self.shapefiles))) + self.shapefiles = [f.replace(".shp", "") for f in + glob(self.shapefiles)] + LOG.info("SHP %s", str(self.shapefiles)) + + def aplot(self, xloc, yloc, map_data=None, wind_data=None, + contour_data=None, scatter_data=None, **kwargs): + """ Plot one or more maps directly from numerical arrays """ + + # store options to be restored at the end + localvars = copy.deepcopy(vars(self)) + + self.init_environ(**kwargs) if scatter_data is not None: if isinstance(scatter_data, str): - scatter_data = DataFrameHandler(parse_path(self.indir, scatter_data)).filter() + scatter_data = DataFrameHandler( + parse_path(self.indir, scatter_data)).filter() else: scatter_data = DataFrameHandler(scatter_data).filter() else: scatter_data = None - log.info("scatter data: {}".format(str(scatter_data))) + LOG.info("scatter data: %s", str(scatter_data)) data = MapData( - map_data=[map_data], #temporarily + map_data=[map_data], # temporarily wind_data=wind_data, contour_data=contour_data, scatter_data=scatter_data ) - grid = (x, y) + grid = (xloc, yloc) self.init_map(grid) if self.img_template: @@ -970,24 +821,29 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi fig_name = "%s/%s" % (self.outdir, f_name) if self.subplot is None: plt.clf() - self.mgplot, self.mgaxis = plt.subplots() + self.mgaxis = plt.axes(projection=self._crs) + self.mgplot = plt.gcf() else: - self.mgaxis = plt.subplot(self.subplot[0], self.subplot[1], self.subplot[2]) + self.mgaxis = plt.subplot(self.subplot[0], + self.subplot[1], + self.subplot[2], + projection=self._crs) self.mgplot = plt.gcf() - fn = self.gen_image_map( + self.gen_image_map( fig_name, grid, data, self.title, - #scatter_data=scatter_data + # scatter_data=scatter_data ) # # Create Max at required intervals # if self.maxdata and n_time in self.maxdata: # if (self.maxtitle != None): # fig_name = "%s_%s" % (f_name, 'MAXD%d' % (n_time/24)) -# p_title = self.set_title(self.maxtitle, s_date, curr_date, s_time) +# p_title = self.set_title(self.maxtitle, s_date, curr_date, +# s_time) # self.gen_image_map( # fig_name, # grid, @@ -1004,7 +860,7 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi """ Plot cube """ if 'title' not in kwargs: - kwargs['title'] = '{0.long_name} ({0.units})'.format(cube) + kwargs['title'] = self._get_default_title(cube) self.aplot( cube.coord('longitude').points, @@ -1018,41 +874,30 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi # store options to be restored at the end localvars = copy.deepcopy(vars(self)) - # update from config file - if ('config' in kwargs) and ('section' in kwargs): - self.load_conf(kwargs['section'], kwargs['config']) - # update from command line - vars(self).update(parse_parameters(kwargs)) - # arguments are only local + + self.init_environ(**kwargs) if isinstance(self.srcvars, str): self.srcvars = [self.srcvars] elif not isinstance(self.srcvars, list): - log.error('Error') + LOG.error('Error') if isinstance(self.srcfiles, str): self.srcfiles = [self.srcfiles] elif not isinstance(self.srcfiles, list): - log.error('Error') - - # parse wind arguments - if (self.windopts and (len(self.windopts) == 4)): - self.wind_scale = int(self.windopts[3]); - log.info("SHP {}".format(str(self.shapefiles))) - if self.shapefiles: - self.shapefiles = [f.replace(".shp","") for f in glob(self.shapefiles)] - log.info("SHP {}".format(str(self.shapefiles))) + LOG.error('Error') if self.scatter is not None: if isinstance(self.scatter, str): - scatter_data = DataFrameHandler(parse_path(self.indir, self.scatter)).filter() + scatter_data = DataFrameHandler( + parse_path(self.indir, self.scatter)).filter() else: scatter_data = DataFrameHandler(self.scatter).filter() else: scatter_data = None - log.info("scatter data: {}".format(str(scatter_data))) + LOG.info("scatter data: %s", str(scatter_data)) - #run = '' + # run = '' run_tmp = '' map_names = [] @@ -1082,71 +927,84 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi if self.timesteps == 'all': self.timesteps = list(sorted(dims.keys())) - run_tmp = dims[0] #12:30Z30NOV2010 + run_tmp = dims[0] # 12:30Z30NOV2010 if run_tmp[:-10] == "01:30": run_tmp = run_tmp.replace("01:30", "00") - #run = "%sh %s %s %s" % (run_tmp[:-10], run_tmp[-9:-7], run_tmp[-7:-4], run_tmp[-4:]) + # run = "%sh %s %s %s" % (run_tmp[:-10], run_tmp[-9:-7], + # run_tmp[-7:-4], run_tmp[-4:]) # return grid grid = nc_handler.grid - #if not NOMAP: + # if not NOMAP: self.init_map(grid) s_date = datetime.strptime("%s %s %s %s" % (run_tmp[-4:], run_tmp[-7:-4], run_tmp[-9:-7], - run_tmp[0:2]), "%Y %b %d %H") + run_tmp[0:2]), + "%Y %b %d %H") ss_date = s_date.strftime("%Y%m%d") steps = [] - START = self.timesteps[0] + start = self.timesteps[0] - #log.info("Start @ %s - s_date: %s - run_tmp: %s" % (START, s_date, run_tmp)) + # LOG.info("Start @ %s - s_date: %s - run_tmp: %s" % (START, s_date, + # run_tmp)) - if START is None: - START = 0 + if start is None: + start = 0 - for n_time in self.timesteps: #range(START,int(TOTAL),INTERVAL): + for n_time in self.timesteps: # range(START,int(TOTAL),INTERVAL): if self.subplot is None: plt.clf() - self.mgplot, self.mgaxis = plt.subplots() + self.mgaxis = plt.axes(projection=self._crs) + self.mgplot = plt.gcf() else: - self.mgplot, self.mgaxis = plt.subplot(self.subplot[0], - self.subplot[1], - self.subplot[3]) - + self.mgaxis = plt.subplot(self.subplot[0], + self.subplot[1], + self.subplot[3], + projection=self._crs) + self.mgplot = plt.gcf() valid_tmp = dims[n_time] - curr_date = datetime.strptime("%s %s %s %s" % (valid_tmp[-4:], valid_tmp[-7:-4], valid_tmp[-9:-7], valid_tmp[0:2]), "%Y %b %d %H" ) + curr_date = datetime.strptime("%s %s %s %s" % (valid_tmp[-4:], + valid_tmp[-7:-4], + valid_tmp[-9:-7], + valid_tmp[0:2]), + "%Y %b %d %H") s_time = "%02d" % n_time steps.append(s_time) p_time = (curr_date - s_date).total_seconds()/3600 - s_time24 = "%02d" % (p_time%24) - #s_time3d = "%03d" % p_time # Currently unused, only for Valentina who has more than 99 img in a set + # s_time24 = "%02d" % (p_time%24) + # s_time3d = "%03d" % p_time # Currently unused, only for + # Valentina who has more than 99 img in a set if self.img_template: f_name = self.img_template % {'date': ss_date, 'step': s_time} fig_name = "%s/%s" % (self.outdir, f_name) - loop_name = self.img_template % {'date': ss_date, 'step': 'loop'} + loop_name = self.img_template % {'date': ss_date, + 'step': 'loop'} else: - f_name = '.'.join(os.path.basename(self.srcfiles[0]).split('.')[:-1]) + f_name = '.'.join(os.path.basename(self.srcfiles[0]) + .split('.')[:-1]) fig_name = "%s/%s_%s" % (self.outdir, f_name, s_time) loop_name = "%s_loop" % (f_name) - #valid = "%02dUTC %s" % (p_time % 24, curr_date.strftime("%d %b %Y")) + # valid = "%02dUTC %s" % (p_time % 24, curr_date.strftime("%d %b + # %Y")) - #stime = "%02d" % (int(run_tmp.split('Z')[0])+p_time) + # stime = "%02d" % (int(run_tmp.split('Z')[0])+p_time) stime = "%02d" % (p_time) - p_title = self.set_title(self.title, s_date, curr_date, s_time, stime) + p_title = set_title(self.title, s_date, curr_date, s_time, stime) cur_scatter_data = None - if scatter_data != None and (curr_date in scatter_data): + if scatter_data is not None and (curr_date in scatter_data): cur_scatter_data = scatter_data[curr_date] - fn = self.gen_image_map( + fname = self.gen_image_map( fig_name, grid, nc_handler.get_data_for_tstep(n_time), @@ -1154,14 +1012,15 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi cur_scatter_data=cur_scatter_data ) - fig_names.append(fn) + fig_names.append(fname) plt_names.append((self.mgplot, self.mgaxis)) # Create Max at required intervals if self.maxdata and n_time in self.maxdata: - if self.maxtitle != None: + if self.maxtitle is not None: fig_name = "%s_%s" % (f_name, 'MAXD%d' % (n_time/24)) - p_title = self.set_title(self.maxtitle, s_date, curr_date, s_time, stime) + p_title = set_title(self.maxtitle, s_date, curr_date, + s_time, stime) self.gen_image_map( fig_name, grid, @@ -1169,19 +1028,20 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi p_title ) else: - print("Missing MAXTITLE, thus not generating MAX images!") + LOG.info("Missing MAXTITLE, not generating MAX images!") if self.kml or self.kmz: plt.clf() run_date = s_date.strftime("%Y%m%d%H") - #separate colorbar + # separate colorbar mpl.rcParams['axes.linewidth'] = 0.1 fig = plt.figure(figsize=(.1, 12)) - ax = fig.add_axes([0.05, 0.80, 0.9, 0.15], axisbg=(1, 1, 1, 0)) + cax = fig.add_axes([0.05, 0.80, 0.9, 0.15]) + # , axisbg=(1, 1, 1, 0)) print("...", self.bounds, "...") - cb = mpl.colorbar.ColorbarBase( - ax, + cbar = mpl.colorbar.ColorbarBase( + cax, # =self.mgaxis, cmap=self.cmap, norm=self.norm, values=self.bounds, @@ -1189,16 +1049,16 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi extend=self.extend, drawedges=False) if self.bounds: - cb.set_ticklabels(self.bounds) - else: - try: - cb.set_ticklabels(mco.levels) - except: - pass - cb.ax.tick_params(labelsize=6) - plt.setp(plt.getp(ax, 'yticklabels'), color='w') - plt.setp(plt.getp(ax, 'yticklabels'), fontsize=6) - for lin in cb.ax.yaxis.get_ticklines(): + cbar.set_ticklabels(self.bounds) +# else: +# try: +# cbar.set_ticklabels(mco.levels) +# except: +# pass + cbar.ax.tick_params(labelsize=6) + plt.setp(plt.getp(cax, 'yticklabels'), color='w') + plt.setp(plt.getp(cax, 'yticklabels'), fontsize=6) + for lin in cbar.ax.yaxis.get_ticklines(): lin.set_visible(False) tit = str(p_title) @@ -1209,46 +1069,46 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi var_unit = '' plt.xlabel( p_title, - #"%s\n%s" % (self.srcvars[0], var_unit), + # "%s\n%s" % (self.srcvars[0], var_unit), horizontalalignment='left', color='w', fontsize=6, ) - ax.xaxis.set_label_coords(-0.025, -0.025) + cax.xaxis.set_label_coords(-0.025, -0.025) + + os.makedirs("%s-%s" % (run_date, self.srcvars[0])) - try: - os.mkdir("%s-%s" % (run_date, self.srcvars[0])) - except: - pass if self.save: fig.savefig( - "%s/%s-%s/%s-colorbar.png" % (self.outdir, run_date, self.srcvars[0], run_date), + "%s/%s-%s/%s-colorbar.png" % (self.outdir, run_date, + self.srcvars[0], run_date), bbox_inches='tight', pad_inches=0, dpi=self.dpi, transparent=True ) - #generate KMZ - Offline - if self.kmz: - print("Generating KMZ ...") - self.genKML(fig_names, self.srcvars, dims, online=False) - - #generate KML - Online - if self.kml: - print("Generating KML ...") - self.genKML(fig_names, self.srcvars, dims) + # generate KMZ/KML - Offline/Online - online priority + if self.kmz or self.kml: + online = self.kml + LOG.info("Generating KMZ ...") + gen_kml(fig_names, self.srcvars[0], self.lon, self.lat, dims, + self.outdir, online=online) - #generate animation. + # generate animation. if self.anim: - self.gen_anim( + gen_anim( "%s.%s" % (loop_name.replace('loop', '*'), self.filefmt), - "%s.gif" % loop_name + "%s.gif" % loop_name, + self.outdir, + self.anim_delay ) -# fulloutdir = os.path.join(self.outdir, os.path.dirname(self.img_template)) +# fulloutdir = os.path.join(self.outdir, +# os.path.dirname(self.img_template)) # loop_name = os.path.basename(loop_name) # self.gen_anim( # fulloutdir, -# "%s.%s" % (loop_name.replace('loop','*'), self.filefmt), +# "%s.%s" % (loop_name.replace('loop','*'), +# self.filefmt), # fulloutdir, # "%s.gif" % loop_name # ) @@ -1258,20 +1118,18 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi # restore options without local function parameters vars(self).update(localvars) - #print "Returning ", plt_names - #return plt_names + # print "Returning ", plt_names + # return plt_names # num = len(map_names) # if num > 1: # mg.compare_maps(map_names, OUTDIR, ss_date, steps, JOINT_TEMPLATE ) - def reset_conf(self): """ Back to the initial conditions. """ self.__init__() - def load_conf(self, section=None, fpath=None, reset=False): """ Load existing configurations from file. """ @@ -1283,7 +1141,7 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi fpath = parse_path(self.indir, fpath) if not os.path.isfile(fpath): - print("Error %s" % fpath) + LOG.error("Error %s", fpath) return opts = read_conf(section, fpath) @@ -1294,7 +1152,6 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi self.reset_conf() vars(self).update(parse_parameters(opts)) - def write_conf(self, section, fpath=None): """ Write configurations on file. """ @@ -1316,7 +1173,8 @@ f=image" % (self.lon[0], self.lat[0], self.lon[-1], self.lat[-1], size, self.dpi self.reset_conf() defaults = copy.deepcopy(vars(self)) # calculate diff - diff = {k : curropts[k] for k in curropts if curropts[k] != defaults[k]} + diff = {k: curropts[k] for k in curropts if curropts[k] != + defaults[k]} # reload diff vars(self).update(diff) # save diff diff --git a/mapgenerator/plotting/timeseries.py b/mapgenerator/plotting/timeseries.py new file mode 100644 index 0000000000000000000000000000000000000000..5fbf7b9b22e2176b76a09edbfe273b12f1013ec0 --- /dev/null +++ b/mapgenerator/plotting/timeseries.py @@ -0,0 +1,206 @@ +from ast import parse +import math +import os +import logging +from matplotlib.dates import date2num, num2date + +import matplotlib.pyplot as plt +from matplotlib.gridspec import GridSpec +import matplotlib.dates as mdates +import netCDF4 + +from mapgenerator.plotting.definitions import MapGenerator + +logger = logging.getLogger(__name__) + +class PlotSeries(MapGenerator): + """ Main class for plotting time series """ + + def __init__(self, loglevel='WARNING', **kwargs): + """ Initialize class with attributes """ + super().__init__(loglevel=loglevel, kwargs=kwargs) + self._current_fig = None + self.scalex = kwargs.get('scalex', None) + self.scaley = kwargs.get('scaley', None) + self.suptitle_fontsize = kwargs.get('suptitle_fontsize', 24) + self.title_fontsize = kwargs.get('title_fontsize', 16) + self.axis_fontsize = kwargs.get('axis_fontsize', 12) + self.plot_size = kwargs.get("plot_size", (8.0, 6.0)) + + def _close(self): + plt.close(self._current_fig) + self._current_fig = None + + def plot_cube(self, cube, coord=None, **kwargs): + if coord: + coord = cube.coord(coord) + if not coord and cube.dim_coords: + coord = cube.dim_coords[0] + if not coord and cube.dim_coords: + coord = cube.aux_coords[0] + + if coord.units.calendar: + points = date2num(coord.units.num2date(coord.points)) + else: + points = coord.points + + if 'xlabel' not in kwargs: + kwargs['xlabel'] = self._get_default_title(cube) + self._current_fig = plt.figure(figsize=self.plot_size) + self._plot_array( + points, + cube.data, + title=self._get_default_title(cube), + **kwargs, + ) + self._set_time_axis(coord) + self._current_fig.tight_layout() + suptitle = kwargs.pop('suptitle', None) + if suptitle: + plt.suptitle(suptitle, y=1.08, fontsize=self.suptitle_fontsize) + self._save_fig(self.img_template) + self._close() + + def _plot_array(self, x, y, **kwargs): + invertx = kwargs.pop('invertx', False) + inverty = kwargs.pop('inverty', False) + xlabel = kwargs.pop('xlabel', None) + ylabel = kwargs.pop('ylabel', None) + xlimits = kwargs.pop('xlimits', None) + ylimits = kwargs.pop('ylimits', None) + title = kwargs.pop('title', None) + plt.plot(x, y) + ax = plt.gca() + if xlabel: + ax.set_xlabel(xlabel, fontsize=self.axis_fontsize) + if ylabel: + ax.set_ylabel(ylabel, fontsize=self.axis_fontsize) + if title: + ax.set_title(title, fontsize=self.title_fontsize, y=1.04) + if invertx: + logger.debug('Invert x axis') + ax.invert_xaxis() + if inverty: + logger.debug('Invert y axis') + ax.invert_yaxis() + if self.scalex: + ax.set_xscale(self.scalex) + if self.scaley: + ax.set_yscale(self.scaley) + if xlimits: + if xlimits == 'tight': + ax.set_xlim(left=x.min(), right=x.max()) + elif xlimits == 'auto': + ax.set_autoscalex_on(True) + else: + ax.set_xlim(left=xlimits[0], right=xlimits[1]) + if ylimits: + if ylimits == 'tight': + ax.set_ylim(left=y.min(), right=y.max()) + elif ylimits == 'auto': + ax.set_autoscaley_on(True) + else: + ax.set_ylim(left=ylimits[0], right=ylimits[1]) + ax.grid(b=True, which='major', axis='both', alpha=0.6) + ax.grid(b=True, which='minor', axis='both', alpha=0.3) + ax.tick_params(axis='both', which='major', labelsize=self.axis_fontsize) + + + def multiplot_cube(self, cube, coord, multi_coord, ncols=2, invert=False, **kwargs): + coord = cube.coord(coord) + multi_coord = cube.coord(multi_coord) + if multi_coord.shape[0] == 1: + self.plot_cube(cube, coord, **kwargs) + return + sharex = kwargs.get('sharex', False) + sharey = kwargs.get('sharey', False) + + self._current_fig = plt.figure( + figsize=( + ncols * self.plot_size[0], + math.ceil(multi_coord.shape[0] / ncols) * self.plot_size[1]) + ) + suptitle = kwargs.pop('suptitle', '') + if suptitle: + suptitle = f"{suptitle}\n\n{self._get_default_title(cube)}" + else: + suptitle = self._get_default_title(cube) + self._current_fig.suptitle( + suptitle, y=1.0, fontsize=self.suptitle_fontsize + ) + gs = GridSpec(math.ceil(multi_coord.shape[0] / ncols), ncols) + for i, plot_cube in enumerate(cube.slices_over(multi_coord)): + if i == 0 or not (sharex or sharey): + self._current_fig.add_subplot(gs[i]) + elif sharex: + if sharey: + self._current_fig.add_subplot(gs[i], sharex=plt.gca(), sharey=plt.gca()) + kwargs.pop('invertx', None) + kwargs.pop('inverty', None) + else: + self._current_fig.add_subplot(gs[i], sharex=plt.gca()) + kwargs.pop('invertx', None) + elif sharey: + self._current_fig.add_subplot(gs[i], sharey=plt.gca()) + kwargs.pop('inverty', None) + title = plot_cube.coord(multi_coord).cell(0).point + if isinstance(title, bytes): + title = title.decode() + plt.title( + title, + fontsize=self.title_fontsize + ) + if coord.units.calendar: + points = date2num(coord.units.num2date(coord.points)) + else: + points = coord.points + if invert: + x = plot_cube.data + y = points + if 'ylabel' not in kwargs: + kwargs['ylabel'] = self._get_default_title(coord) + else: + x = points + y = plot_cube.data + if 'xlabel' not in kwargs: + kwargs['xlabel'] = self._get_default_title(coord) + + self._plot_array( + x, y, **kwargs, + ) + self._set_time_axis(coord) + + self._current_fig.tight_layout(pad=2.0) + self._save_fig(self.img_template) + self._close() + + def _save_fig(self, name): + fullname = os.path.join(self.outdir ,f"{name}.{self.filefmt}") + self._current_fig.savefig(fullname, bbox_inches='tight', pad_inches=.2, dpi=self.dpi) + + + + @staticmethod + def _set_time_axis(coord): + if not coord.units.calendar: + return + axis = plt.gca().xaxis + years = coord.cell(-1).point.year - coord.cell(0).point.year + if years < 10: + major_locator = 1 + minor_locator = None + elif years < 50: + major_locator = 5 + minor_locator = 1 + elif years < 100: + major_locator = 10 + minor_locator = 2 + else: + major_locator = 20 + minor_locator = 5 + axis.set_major_locator(mdates.YearLocator(major_locator)) + if minor_locator: + axis.set_minor_locator(mdates.YearLocator(minor_locator)) + axis.set_major_formatter(mdates.DateFormatter('%Y')) + axis.label._text = f"{coord.name()} (years)" + diff --git a/mapgenerator/plotting/tools.py b/mapgenerator/plotting/tools.py new file mode 100644 index 0000000000000000000000000000000000000000..ce3708c9b4a11cf41ae3596eaa7929e5ad18687f --- /dev/null +++ b/mapgenerator/plotting/tools.py @@ -0,0 +1,213 @@ +from datetime import datetime +import subprocess +import os +import logging + +LOG = logging.getLogger(__name__) + + +def print_time(msg=None): + """ Print time """ + # initial time + if 'start_t' not in globals(): + global start_t, last_t + LOG.info("Creating start_t") + start_t = datetime.now() + last_t = start_t + + now_t = datetime.now() + diff_t = now_t - last_t + if msg is not None: + LOG.info('TIME: %s done in %s s', msg, diff_t.seconds) + last_t = now_t + + +def run_command(comm_str): # , fatal=True): + """ Run command """ + LOG.info(comm_str) + proc = subprocess.Popen(comm_str.split(), stderr=subprocess.PIPE) + output, err = proc.communicate() + LOG.info(output, '-----', err) +# st, out = commands.getstatusoutput(comm_str) +# if st != 0: +# print "Error: %s" % str(out) +# if (fatal == True): +# sys.exit(1) + + +def gen_anim(inpattern, outfile, outdir, anim_delay): + ''' Generate an animation starting from a set of + images, specified by the inpattern parameter ''' + # print "Animation: ", indir, inpattern, outdir, outfile + # Create animation + run_command( + "/usr/bin/convert -delay %s -loop 0 -layers Optimize %s/%s %s/%s" + % (anim_delay, outdir, inpattern, outdir, outfile)) + # Remove intermediate files + # run_command("rm %s/%s %s" % (indir, inpattern, map_name), + # False) + +# ccc = [] +# ccc.extend(colors) +# print "Filtering scatter data for date/hour: ", +# curr_date.strftime("%Y%m%d"), s_time24 csv_lat, csv_lon, csv_val = +# csv_handler.filter(curr_date.strftime("%Y%m%d"), s_time24) +# #csv_col = dataTransform.setColorsFromBuckets(ccc, +# bounds, csv_con) +# if(len(csv_lat) > 0): +# scatter_data = [csv_lon, csv_lat, 40, csv_col] +# print "**************** SCATTER_DATA", scatter_data +# print "Filtered %d records for scatterplot" % len(scatter_data) + + +def set_resolution(lon, lat): + """ 110m, 50m, 10m """ + # lats = 36, 72, 108, 144, 180 + # lons = 72, 144, 216, 288, 360 + minlat, maxlat = lat[0], lat[-1] + minlon, maxlon = lon[0], lon[-1] + lats = maxlat + abs(minlat) + lons = maxlon + abs(minlon) + if lats > 100 and lons > 200: + resolution = '110m' + zoom = 1 + elif lats > 45 and lons > 75: + resolution = '50m' + elif lats <= 45 and lons <= 75: + resolution = '10m' + zoom = 4 + else: + resolution = '110m' # default + zoom = 8 + return resolution, zoom + + +def set_title(title, sdate, cdate, step, stime): + """ Set the title of the current image to the one provided, substituting the + patterns """ + + title_dic = { + 'syear': sdate.strftime("%Y"), + 'smonth': sdate.strftime("%m"), + 'sMONTH': sdate.strftime("%b"), + 'sday': sdate.strftime("%d"), + 'shour': sdate.strftime("%H"), + 'year': cdate.strftime("%Y"), + 'month': cdate.strftime("%m"), + 'MONTH': cdate.strftime("%b"), + 'day': cdate.strftime("%d"), + 'hour': cdate.strftime("%H"), + 'step': step, + 'simday': "%d" % (int(stime)/24), + 'simhh': stime, + # 'simhh': "%02d" % (int(stime)%24) + } + + try: + p_title = str(title) % title_dic # ("%sh" % s_time, valid) + except Exception as err: + LOG.error("Title error: %s", str(err)) + p_title = title % title_dic # ("%sh" % s_time, valid) + + return p_title + + +def gen_kml(fig_names, var_name, lon, lat, dims, outdir='.', online=True): + """ Generate KML/KMZ files """ + from lxml import etree + from pykml.factory import KML_ElementMaker as KML + import zipfile + + LOG.info(fig_names) + # kml_name = '-'.join(os.path.basename(fig_names[0]).split('-')[:-1]) + kml_name = os.path.basename(fig_names[0])[:-3] + run_date0 = datetime.strptime(dims[0], "%HZ%d%b%Y") + run_date1 = datetime.strptime(dims[1], "%HZ%d%b%Y") + begindate = run_date0.strftime("%Y%m%d%H") + dir_name = os.path.join(outdir, begindate + '-' + + var_name).replace('./', '') + interval = run_date1 - run_date0 + datefmt = "%Y-%m-%dT%H:00:00Z" + + doc = KML.kml( + KML.Folder( + KML.name(kml_name), + KML.LookAt( + KML.longitude((lon[1]+lon[-1])/2), + KML.latitude((lat[1]+lat[-1])/2), + KML.range(8500000.0), + KML.tilt(0), + KML.heading(0), + ), + KML.ScreenOverlay( + KML.name("Legend"), + KML.open(1), + KML.Icon( + KML.href("%s/%s-colorbar.png" % (dir_name, begindate)), + ), + KML.overlayXY(x="0", y="1", xunits="fraction", + yunits="fraction"), + KML.screenXY(x=".01", y="0.55", xunits="fraction", + yunits="fraction"), + KML.rotationXY(x="0", y="0", xunits="fraction", + yunits="fraction"), + KML.size(x="0", y="0.5", + xunits="fraction", + yunits="fraction"), + id="%s-ScreenOverlay" % begindate, + ), + ) + ) + + os.makedirs(dir_name) + + if not online: + zfile = zipfile.ZipFile("%s.kmz" % kml_name, 'w') + + for fig_name, dat in zip(fig_names, sorted(dims.keys())): + dtime = dims[dat] + fig_name = fig_name + ".png" + fig_path = os.path.join(dir_name, os.path.basename(fig_name)) + startdate = datetime.strptime(dtime, "%HZ%d%b%Y") + begdate = startdate.strftime(datefmt) + enddate = (startdate + interval).strftime(datefmt) + starth = startdate.hour + doc.Folder.append(KML.GroundOverlay( + KML.name("%02d:00:00Z" % starth), + KML.TimeSpan( + KML.begin(begdate), + KML.end(enddate), + ), + KML.Icon( + KML.href(fig_path), + KML.viewBoundScale(1.0), + ), + KML.altitude(0.0), + KML.altitudeMode("relativeToGround"), + KML.LatLonBox( + KML.south(lat[0]), + KML.north(lat[-1]), + KML.west(lon[0]), + KML.east(lon[-1]), + KML.rotation(0.0), + ), + )) + + if os.path.exists(fig_name): + os.rename(fig_name, fig_path) + zfile.write(fig_path) + + outf = open("%s.kml" % kml_name, 'w') + outf.write('\n') + outf.write(etree.tostring(doc, pretty_print=True)) + outf.close() + + if not online: + zfile.write("%s.kml" % kml_name) + zfile.write("%s/%s-colorbar.png" % (dir_name, begindate)) + zfile.close() + + for img in os.listdir(dir_name): + os.remove("%s/%s" % (dir_name, img)) + os.rmdir(dir_name) + os.remove("%s.kml" % kml_name) diff --git a/setup.py b/setup.py index 81b93a9f8f5dd3f3ee98ad7688680ab4734723fa..966d4f7e3b0b2e86cae19efb341b0366f5fd27b7 100644 --- a/setup.py +++ b/setup.py @@ -72,7 +72,7 @@ setup( install_requires=[ "matplotlib", "pandas", - "Basemap", + "Cartopy", "numpy", "netCDF4", "ConfigArgParse",