first changes

common/__init__.py

@@ -0,0 +1,3 @@
+from .inout import read_input
+from .boundarybox import BoundaryBox
+from .vector import Vector, NumpyVector

common/boundarybox.py

@@ -0,0 +1,20 @@
+class BoundaryBox:
+    lat_min: float
+    lat_max: float
+    lon_min: float
+    lon_max: float
+
+    def __init__(self, lat_min=-90., lat_max=90., lon_min=-180., lon_max=180.):
+        self.lat_min = lat_min
+        self.lat_max = lat_max
+        self.lon_min = lon_min
+        self.lon_max = lon_max
+
+    def middle_lat(self):
+        x = self.lat_min + 0.5 * (self.lat_max - self.lat_min)
+        return x
+
+    def middle_lon(self):
+        return self.lon_min + 0.5 * (self.lon_max - self.lon_min)
+
+

common/drawcurrents.py

@@ -0,0 +1,81 @@
+
+
+def drawcurrents(coordinates_rank, nx, ny, scale, resolution,
+                 level, time, lat, lon, ust, vst, mod,
+                 file_out, title, style, boundary_box):
+    """
+
+    :return:
+    """
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from mpl_toolkits.basemap import Basemap
+
+    fig = plt.figure(figsize=(10, 10))
+    ax = fig.add_subplot(111)
+
+    dx = 0.01
+    middle_lon = boundary_box.middle_lon()
+    middle_lat = boundary_box.middle_lat()
+    m = Basemap(llcrnrlon=boundary_box.lon_min - dx,
+                llcrnrlat=boundary_box.lat_min - dx,
+                urcrnrlon=boundary_box.lon_max + dx,
+                urcrnrlat=boundary_box.lat_max + dx,
+                resolution=resolution, projection='tmerc', lon_0=middle_lon, lat_0=middle_lat)
+
+    m.drawcoastlines()
+    m.fillcontinents(color='grey', lake_color='aqua')
+
+
+    #m.drawmapboundary(fill_color='aqua')
+
+    if coordinates_rank == 1:
+        lon, lat = np.meshgrid(lon, lat)
+
+    x, y = m(lon, lat)
+
+    # draw filled contours.
+    clevs = [0, 0.1, 0.2, 0.3, 0.4, 0.5]
+    # clevs = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
+    if style == 'contour':
+        cs = m.contourf(x, y, mod, clevs, cmap=plt.cm.jet)
+        plt.quiver(x[::nx, ::ny], y[::nx, ::ny], ust[::nx, ::ny], vst[::nx, ::ny], scale=scale)
+
+    if style == 'cvector':
+        clim = [0., 0.5]
+        cs = plt.quiver(x[::nx, ::ny], y[::nx, ::ny], ust[::nx, ::ny], vst[::nx, ::ny], mod[::nx, ::ny], clim=clim, scale=scale)
+
+    if style == 'windbarbs':
+        clim = [0., 10]
+        cs = plt.barbs(x[::nx, ::ny], y[::nx, ::ny], ust[::nx, ::ny], vst[::nx, ::ny], mod[::nx, ::ny], clim=clim,cmap=plt.cm.jet)
+
+    if style == 'stream':
+        print(lons.shape)
+        print(x[::, 1])
+        print(ust.shape)
+        xx = x[1::, ::]
+        xxx = xx[0]
+        yyy = y[::, 1::][1]
+
+        u = ust[::, ::]
+        v = vst[::, ::]
+        print(len(xxx))
+        print(len(yyy))
+        print(len(u))
+        print(len(v))
+
+        cs = plt.streamplot(lon, lat, u, v)
+
+    cbar = m.colorbar(cs, location='bottom', pad="5%")
+    cbar.ax.tick_params(labelsize=9)
+    cbar.set_label('m/s')
+
+    # add title
+    plt.title(title)
+
+    fig.savefig(file_out, dpi=100, facecolor='w', edgecolor='w', format='png',
+                transparent=False,  bbox_inches='tight', pad_inches=0.1)
+    plt.clf()
+    plt.close(fig)
+
+    return

common/drawmap.json

@@ -0,0 +1,28 @@
+{
+    "path_in": "",
+    "path_out": "",
+    "file_in": "",
+    "file_out": "",
+    "nx": 1,
+    "ny": 1,
+    "vector": 0,
+    "scalar": 1,
+    "u": "eastward_velocity",
+    "v": "northward_velocity",
+    "scalar_magnitude": "",
+    "resolution": "i",
+    "scale": 5,
+    "n_time": "all",
+    "n_level": 9,
+    "title": "",
+    "style": "cvector",
+    "limits": [
+        42.08,
+        42.45,
+        -9.004,
+        -8.601
+    ],
+    "wms_url": "https://ows.emodnet-bathymetry.eu/wms?",
+    "wms_layers": ["emodnet:mean_atlas_land", "coastlines"]
+
+}

common/drawmap.py

@@ -0,0 +1,252 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# Para el tratamiento de datos
+import numpy as np
+import pandas as pd
+# Para las representaciones
+import matplotlib.pyplot as plt
+from matplotlib.cm import jet
+# Para los directorios
+import os
+# Importamos funciones de otros programas
+from common.readers.reader_factory import read_factory
+from common import read_input
+from common.boundarybox import BoundaryBox
+from common.drawcurrents import drawcurrents
+from common.map import Map
+
+def read_inputs(input_file):
+    """Read keywords for options"""
+    input_keys = ['path_in',
+                  'file_in',
+                  'path_out'
+                  'file_out',
+                  'nx',
+                  'ny',
+                  'resolution',
+                  'scale',
+                  'n_time',
+                  'n_level',
+                  'title',
+                  'style',
+                  'limits',
+                  'vector',
+                  'scalar',
+                  'wms_url',
+                  'wms_layers']
+    return read_input(input_file, input_keys)
+
+
+def draw_map_vector(inputs, n):
+    """draw 1 maps of a day"""
+    draw_map = Vector(inputs)
+    draw_map.read_head()
+    draw_map.create_title(n)
+    draw_map.reader_by_time()
+    draw_map.draw()
+
+
+def draw_map_scalar(inputs, n):
+    draw_map = Scalar(inputs)
+    draw_map.read_head()
+    draw_map.create_title(n)
+
+    if n is None:
+        draw_map.reader_no_time()
+    else:
+        draw_map.reader_by_time()
+    draw_map.draw()
+
+
+def draw_map_24(inputs,start,end,ref_lon, ref_lat, key, draw_scale):
+    """draw 24+1 maps of a day"""
+
+    if inputs['scalar']:
+        draw_map = Scalar(inputs)
+    elif inputs['vector']:
+        draw_map = Vector(inputs)
+
+    draw_map.read_head()
+    draw_map.options.time = 0 # creo que esto sobra
+
+    for n in range(draw_map.reader.ini_ntime+start, draw_map.reader.ini_ntime+end):
+        draw_map.create_title(n)
+        draw_map.options.time = n
+        draw_map.reader_by_time()
+        draw_map.draw(ref_lon, ref_lat, key, draw_scale)
+
+
+class OptionsMap:
+
+    def __init__(self, inputs):
+
+        self.file_path_in = inputs['path_in']
+        self.file_path_out = inputs['path_out']
+        self.file_in = inputs['file_in']
+        self.file_name = os.path.join(self.file_path_in, self.file_in)
+        self.file_hdf_out = inputs['file_out']
+        self.file_out = os.path.join(self.file_path_out, self.file_hdf_out)
+
+        self.nx = inputs['nx']
+        self.ny = inputs['ny']
+        self.scale = inputs['scale']
+        self.resolution = inputs['resolution']
+        self.style = inputs['style']
+        self.title = inputs['title']
+
+        self.level = inputs['n_level']
+        self.time = inputs['n_time']
+        limits = inputs['limits']
+        self.boundary_box = BoundaryBox(limits[0], limits[1], limits[2], limits[3])
+
+        self.vector_bool = inputs['vector']
+        self.scalar_bool = inputs['scalar']
+
+        self.wms_url = inputs['wms_url']
+        self.wms_layers = inputs['wms_layers']
+
+
+class DrawMap:
+
+    def __init__(self, inputs):
+        self.options = OptionsMap(inputs)
+        self.reader = self.get_reader(self.options.file_name)
+
+    def read_head(self):
+
+        with self.reader.open():
+
+            lat = self.reader.latitudes
+            lon = self.reader.longitudes
+
+            if self.reader.coordinates_rank == 1:
+                self.lats = lat[0:self.reader.n_latitudes - 2]  # ATENCIóN:Esto y lo de abajo era -1, revisar
+                self.lons = lon[0:self.reader.n_longitudes - 2]
+            elif self.reader.coordinates_rank == 2:
+                self.lats = lat[0:self.reader.n_longitudes - 2, 0:self.reader.n_latitudes - 2]
+                self.lons = lon[0:self.reader.n_longitudes - 2, 0:self.reader.n_latitudes - 2]
+
+    def create_title(self, n_time):
+        if n_time is None:
+            n_time = 0
+        with self.reader.open():
+            data = self.reader.get_date(n_time)
+            data_str = data.strftime("%Y-%m-%d %H:%M UTC")
+            data_comp = data.strftime("%Y%m%d%H%M")
+            self.title_full = self.options.title + " " + data_str
+            self.file_out_full = self.options.file_out + '_' + data_comp + '.png'
+
+    def get_reader(self, file_in):
+        print('Opening: {0}'.format(file_in))
+        factory = read_factory(file_in)
+        return factory.get_reader()
+
+
+class Vector(DrawMap):
+
+    def __init__(self, inputs):
+
+        super().__init__(inputs)
+
+        self.u_name = inputs['u']
+        self.v_name = inputs['v']
+        self.us = None
+        self.vs = None
+        self.modules = None
+
+    def reader_by_time(self):
+
+        with self.reader.open():
+            u = self.reader.get_variable(self.u_name, self.options.time)
+            v = self.reader.get_variable(self.v_name, self.options.time)
+
+            if len(u.shape) == 3:
+                self.us = u[self.options.level, :-1, :- 1]
+                self.vs = v[self.options.level, :-1, :-1]
+
+            elif len(u.shape) == 2:
+                self.us = u[:-1, :-1]
+                self.vs = v[:-1, :-1]
+
+            self.modules = pow((pow(self.us, 2) + pow(self.vs, 2)), .5)
+
+    def draw(self):
+        drawcurrents(self.reader.coordinates_rank, self.options.nx, self.options.ny,
+                     self.options.scale, self.options.resolution,
+                     self.options.level, self.options.time,
+                     self.lats, self.lons, self.us, self.vs, self.modules,
+                     self.file_out_full, self.title_full, self.options.style,
+                     self.options.boundary_box)
+
+
+class Scalar(DrawMap):
+
+    def __init__(self, inputs):
+
+        super().__init__(inputs)
+
+        self.scalar_name = inputs['scalar_magnitude']
+        self.scalars = None
+
+    def reader_by_time(self):
+
+        with self.reader.open():
+            scalar = self.reader.get_variable(self.scalar_name, self.options.time)
+
+            if len(scalar.shape) == 3:
+                self.scalars = scalar[self.options.level, :-1, :- 1]
+
+            elif len(scalar.shape) == 2:
+                self.scalars = scalar[:-1, :-1]
+
+    def reader_no_time(self):
+        with self.reader.open():
+            scalar = self.reader.get_var(self.scalar_name)
+
+            if len(scalar.shape) == 3:
+                self.scalars = scalar[self.options.level, :-1, :- 1]
+
+            elif len(scalar.shape) == 2:
+                self.scalars = scalar[:-1, :-1]
+
+    def draw(self,ref_lon, ref_lat, key, draw_scale):
+        """
+        :return:
+        """
+        # a Map object is created to deal with
+        mapa = Map()
+        mapa.set_extension(self.options.boundary_box.lat_min, self.options.boundary_box.lat_max,
+                           self.options.boundary_box.lon_min, self.options.boundary_box.lon_max)
+        fig = mapa.add_figure()
+        mapa.add_wms(self.options.wms_url, self.options.wms_layers)
+        #mapa.add_wms("https://wms.mapama.gob.es/sig/Agua/PVERT/2021/wms.aspx?", 
+                       #["Puntos de vertido de depuradoras urbanas (Q2021. Dir 91/271/CEE)"])
+
+
+        #  Draw the scalars
+        if self.reader.coordinates_rank == 1:
+            longitudes, latitudes = np.meshgrid(self.lons, self.lats)
+        else:
+            longitudes, latitudes = self.lons, self.lats
+
+        im = mapa.ax.pcolormesh(longitudes, latitudes, self.scalars[:-1, :-1], shading='auto')  # since lats and lons are n-1xn-1 matrix
+        if draw_scale=='jet':
+            im.set_cmap(jet)
+
+        # add title
+        mapa.ax.set_title(self.title_full)
+        # add point
+        path_discharge = os.getcwd() +'\\common\\Rias_UWWTP.xlsx'
+        discharge_lats = (pd.read_excel(path_discharge, sheet_name = key[1])['Lat'].values)
+        discharge_lons = (pd.read_excel(path_discharge, sheet_name = key[1])['Lon'].values)
+        mapa.ax.scatter(discharge_lons, discharge_lats, marker='.', s=150, color='darkorange',
+                        edgecolors = "black", linewidths = 1.)
+        mapa.ax.scatter(ref_lon, ref_lat, marker='*', s=100, color='red')
+        fig.colorbar(im, ax=mapa.ax)
+        #
+        fig.savefig(self.file_out_full, dpi=100, facecolor='w', edgecolor='w', format='png',
+                    transparent=False, bbox_inches='tight', pad_inches=0.1)
+        plt.clf()
+        plt.close(fig)
+        return