Merge pull request #67 from pixalytics-ltd/plotting-update

Plotting update - upversioned python packages and added log scaling option with associated test case

tools/psy-maps/psy-maps.xml

@@ -1,4 +1,4 @@
-<tool id="psy_maps" name="map plot" version="1.2.1">
+<tool id="psy_maps" name="map plot" version="1.3.0" profile="23.0">
     <description>gridded (lat/lon) netCDF data</description>
     <edam_topics>
       <edam_topic>topic_3855</edam_topic>
@@ -9,15 +9,15 @@
     </edam_operations>
     <requirements>
         <requirement type="package" version="3">python</requirement>
-        <requirement type="package" version="1.2.1">psyplot</requirement>
-        <requirement type="package" version="1.2.0">psy-maps</requirement>
-        <requirement type="package" version="1.2.0">psy-reg</requirement>
-        <requirement type="package" version="1.4.2">netcdf4</requirement>
+        <requirement type="package" version="1.4.3">psyplot</requirement>
+        <requirement type="package" version="1.4.2">psy-maps</requirement>
+        <requirement type="package" version="1.4.0">psy-reg</requirement>
+        <requirement type="package" version="1.6.5">netcdf4</requirement>
     </requirements>
     <command detect_errors="exit_code"><![CDATA[
-       HOME=`pwd`  &&
        python3 '$__tool_directory__/psymap_simple.py'
             '$ifilename' '$variable'
+            --logscale '$adv.logscale'
             --proj '$adv.projection'
             --cmap '$adv.colormap'
             --output image.png
@@ -51,6 +51,10 @@
                     <param name="title" type="text" value="" label="plot title" />
                 </when>
             </conditional>
+            <param name="logscale" type="select" label="log scale the data">
+                <option value="yes">yes</option>
+                <option value="no" selected="true">no</option>
+            </param>
             <param name="projection" type="select">
                 <option value="" selected="true">PlateCarree</option>
                 <option value="robin">Robinson</option>
@@ -162,6 +166,14 @@
             <param name="colormap" value="RdBu_r" />
             <output name="ofilename" ftype="png" file="TS.f2000.T31T31.control.cam.h0.0014-12_ortho.png" compare="sim_size" delta="500"/>
         </test>
+        <test>
+            <param name="ifilename" value="ESACCI-OC-L3S-CHLOR_A-20220601-RESIZED-fv6.0.nc" />
+            <param name="variable" value="chlor_a" />
+            <param name="projection" value="" />
+            <param name="colormap" value="jet" />
+            <param name="logscale" value="yes" />
+            <output name="ofilename" ftype="png" file="ESACCI-OC-L3S-CHLOR_A-20220601-RESIZED-fv6.0.png" compare="sim_size" delta="1500"/>
+        </test>
     </tests>
     <help><![CDATA[
 
@@ -174,37 +186,40 @@ This tool wraps the functionality of ``psyplot.plot.mapplot``.
 .. class:: infomark
 
         The wrappers aims at providing the same functionality as ``psyplot.plot.mapplot`` but it has some limitations.
-        The input file must be in netCDF format and must contain 2D arrays with geographical coordinatates
-        (latitudes, longitudes) given in degrees.
-        TODO: Fill in help.
+        The input file must be in netCDF format and must contain 2D arrays with geographical coordinatates 
+	(latitudes, longitudes) given in degrees.
+        The code may fail when the data is significantly different from what has been tested.
 
 **What it does**
 ----------------
 
 This tools creates an image (png format) corresponding to the visualization on a geographical map of a variable extracted from a
-netCDF file (input file). By default, the projection is ``PlateCarree`` and colormap is ``jet``. These settings can be changed in *Advanced settings*.
+netCDF file (input file). 
+By default, the projection is ``PlateCarree`` and colormap is ``jet``. 
+These settings can be changed in *Advanced settings* alongside an option to log scale the plot.
 
 
 **Usage**
 
 ::
 
   usage: psymap_simple.py [-h] [--proj PROJ] [--cmap CMAP] [--output OUTPUT] [--time TIME]
-                          [--format FORMAT] [--title TITLE] [--ncol NCOL] [--nrow NROW] [-v] input varname
+                          [--format FORMAT] [--title TITLE] [--ncol NCOL] [--nrow NROW] [-l] [-v] input varname
 
 
 Positional arguments:
 ~~~~~~~~~~~~~~~~~~~~~
 
 - **input**:            input filename with geographical coordinates (netCDF format)
-- **varname**:          Specify which variable to plot (case sensitive)
+- **varname**:          specify which variable to plot (case sensitive)
 
 Optional arguments:
 ~~~~~~~~~~~~~~~~~~~~~
 
   -h, --help       show this help message and exit
-  --proj PROJ      Specify the projection on which we draw
-  --cmap CMAP      Specify which colormap to use for plotting
+  -l, --logscale   log scale the data
+  --proj PROJ      specify the projection on which we draw
+  --cmap CMAP      specify which colormap to use for plotting
   --output OUTPUT  output filename to store resulting image (png format)
   -v, --verbose    switch on verbose mode
   --time TIME      list of times to plot for multiple plots
@@ -218,5 +233,6 @@ psyplot and psy-maps can be found at https://psyplot.readthedocs.io/projects/psy
 
     ]]></help>
     <citations>
+        <citation type="doi">10.21105/joss.00363</citation>
     </citations>
 </tool>

tools/psy-maps/psymap_simple.py

@@ -4,6 +4,7 @@
 # usage: psymap_simple.py [-h] [--proj PROJ]
 #                              [--cmap CMAP]
 #                              [--output OUTPUT]
+#                              [-l]
 #                              [-v]
 #                              input varname
 #
@@ -14,6 +15,7 @@
 #
 # optional arguments:
 #  -h, --help       show this help message and exit
+#  -l, --logscale   log scale the data
 #  --proj PROJ      Specify the projection on which we draw
 #  --cmap CMAP      Specify which colormap to use for plotting
 #  --output OUTPUT  output filename to store resulting image (png format)
@@ -26,47 +28,63 @@
 #
 
 import argparse
+import math
 import warnings
 from pathlib import Path
 
 import matplotlib as mpl
+import psyplot.project as psy  # noqa: I202,E402
+import xarray
+
 mpl.use('Agg')
 from matplotlib import pyplot  # noqa: I202,E402
-
-import psyplot.project as psy  # noqa: I202,E402
 from psyplot import rcParams   # noqa: I202,E402
 
 
 class PsyPlot ():
-    def __init__(self, input, proj, varname, cmap, output, verbose=False,
-                 time=[], nrow=1, ncol=1, format="%B %e, %Y",
-                 title=""):
+    def __init__(self, input, varname, output=None, logscale=False, cmap=None,
+                 proj=None, verbose=False, time=None, nrow=None, ncol=None,
+                 format=None, title=None):
         self.input = input
-        self.proj = proj
         self.varname = varname
-        self.cmap = cmap
-        self.time = time
-        if format is None:
-            self.format = ""
-        else:
-            self.format = format.replace('X', '%')
-        if title is None:
-            self.title = ""
+        if proj is None or proj == "":
+            self.proj = "cyl"
         else:
+            self.proj = proj
+        self.cmap = cmap if cmap is not None else "jet"
+        self.time = time if time is not None else []
+        self.ncol = int(ncol) if ncol is not None else int(1)
+        self.nrow = int(nrow) if nrow is not None else int(1)
+
+        # Open dataset
+        ds = xarray.open_dataset(input)[varname]
+        minv = math.log2(ds.data.min())
+        maxv = math.log2(ds.data.max())
+        if title is not None:
             self.title = title
-        if ncol is None:
-            self.ncol = 1
         else:
-            self.ncol = int(ncol)
-        if nrow is None:
-            self.nrow = 1
+            self.title = ds.long_name
+            if len(self.title) > 60:
+                self.title = ds.standard_name
+        del ds
+
+        if logscale:
+            # Check that data range is sufficient for log scaling
+            if maxv < (minv * (10.0 ** 2.0)):
+                print("Not possible to log scale, switching to linear scale")
+                self.bounds = None
+            else:
+                self.bounds = ['log', 2]
+        else:
+            self.bounds = None
+        if format is None:
+            self.format = ""
         else:
-            self.nrow = int(nrow)
+            self.format = "%B %e, %Y"
         if output is None:
             self.output = Path(input).stem + '.png'
         else:
             self.output = output
-        self.verbose = verbose
         if verbose:
             print("input: ", self.input)
             print("proj: ", self.proj)
@@ -77,37 +95,32 @@ def __init__(self, input, proj, varname, cmap, output, verbose=False,
             print("nrow: ", self.nrow)
             print("title: ", self.title)
             print("date format: ", self.format)
+            print("logscale: ", self.bounds)
             print("output: ", self.output)
 
     def plot(self):
+        clabel = '{desc}'
         if self.title and self.format:
             title = self.title + "\n" + self.format
         elif not self.title and self.format:
             title = self.format
         elif self.title and not self.format:
             title = self.title
-        else:
-            title = '%(long_name)s'
+            clabel = self.title
 
-        if self.cmap is None and self.proj is None:
-            psy.plot.mapplot(self.input, name=self.varname,
-                             title=title,
-                             clabel='{desc}')
-        elif self.proj is None or not self.proj:
-            psy.plot.mapplot(self.input, name=self.varname,
-                             title=title,
-                             cmap=self.cmap, clabel='{desc}')
-        elif self.cmap is None or not self.cmap:
+        # Plot with chosen options
+        if self.bounds is None:
             psy.plot.mapplot(self.input, name=self.varname,
+                             cmap=self.cmap,
                              projection=self.proj,
                              title=title,
-                             clabel='{desc}')
+                             clabel=clabel)
         else:
             psy.plot.mapplot(self.input, name=self.varname,
-                             cmap=self.cmap,
+                             cmap=self.cmap, bounds=self.bounds,
                              projection=self.proj,
                              title=title,
-                             clabel='{desc}')
+                             clabel=clabel)
 
         pyplot.savefig(self.output)
 
@@ -119,50 +132,41 @@ def multiple_plot(self):
             title = self.format
         else:
             title = self.title + "\n" + self.format
+
         mpl.rcParams['figure.figsize'] = [20, 8]
         mpl.rcParams.update({'font.size': 8})
         rcParams.update({'plotter.maps.grid_labelsize': 8.0})
-        if self.cmap is None and self.proj is None:
-            m = psy.plot.mapplot(self.input, name=self.varname,
-                                 title=title,
-                                 ax=(self.nrow, self.ncol),
-                                 time=self.time, sort=['time'],
-                                 clabel='{desc}')
-            m.share(keys='bounds')
-        elif self.proj is None or not self.proj:
-            m = psy.plot.mapplot(self.input, name=self.varname,
-                                 title=title,
-                                 ax=(self.nrow, self.ncol),
-                                 time=self.time, sort=['time'],
-                                 cmap=self.cmap, clabel='{desc}')
-            m.share(keys='bounds')
-        elif self.cmap is None or not self.cmap:
+
+        # Plot using options
+        if self.bounds is None:
             m = psy.plot.mapplot(self.input, name=self.varname,
+                                 cmap=self.cmap,
                                  projection=self.proj,
                                  ax=(self.nrow, self.ncol),
                                  time=self.time, sort=['time'],
                                  title=title,
                                  clabel='{desc}')
-            m.share(keys='bounds')
         else:
             m = psy.plot.mapplot(self.input, name=self.varname,
-                                 cmap=self.cmap,
+                                 cmap=self.cmap, bounds=self.bounds,
                                  projection=self.proj,
                                  ax=(self.nrow, self.ncol),
                                  time=self.time, sort=['time'],
                                  title=title,
                                  clabel='{desc}')
-            m.share(keys='bounds')
+
+        m.share(keys='bounds')
 
         pyplot.savefig(self.output)
 
 
-def psymap_plot(input, proj, varname, cmap, output, verbose, time,
+def psymap_plot(input, proj, varname, logscale, cmap, output, verbose, time,
                 nrow, ncol, format, title):
     """Generate plot from input filename"""
 
-    p = PsyPlot(input, proj, varname, cmap, output, verbose, time,
+    p = PsyPlot(input, varname, output, logscale, cmap, proj, verbose, time,
                 nrow, ncol, format, title)
+
     if len(time) == 0:
         p.plot()
     else:
@@ -185,6 +189,10 @@ def psymap_plot(input, proj, varname, cmap, output, verbose, time,
         'varname',
         help='Specify which variable to plot (case sensitive)'
     )
+    parser.add_argument(
+        "--logscale",
+        help='Plot the log scaled data'
+    )
     parser.add_argument(
         '--cmap',
         help='Specify which colormap to use for plotting'
@@ -223,6 +231,12 @@ def psymap_plot(input, proj, varname, cmap, output, verbose, time,
         time = []
     else:
         time = list(map(int, args.time.split(",")))
-    psymap_plot(args.input, args.proj, args.varname, args.cmap,
+
+    if args.logscale == 'no':
+        logscale = False
+    else:
+        logscale = True
+
+    psymap_plot(args.input, args.proj, args.varname, logscale, args.cmap,
                 args.output, args.verbose, time,
                 args.nrow, args.ncol, args.format, args.title)