rough outline of how to plot, save, and pass bytes buffer

pyaerocom/aeroval/_processing_base.py

@@ -208,3 +208,15 @@ def read_ungridded_obsdata(self, obs_name, var_name):
 
         data = col._read_ungridded(var_name)
         return data
+
+    def read_gridded_obsdata(self, obs_name, var_name):
+        """
+        Import gridded observation data, usually satellite data
+
+        Args:
+            obs_name (str): Name of observation network in :attr:`cfg`
+            var_name (str): Name of variable to be read.
+        """
+        col = self.get_colocator(obs_name=obs_name)
+        data = col._read_gridded(var_name)
+        return data

pyaerocom/aeroval/modelmaps_engine.py

@@ -3,7 +3,12 @@
 
 from pyaerocom import GriddedData, TsType
 from pyaerocom.aeroval._processing_base import DataImporter, ProcessingEngine
-from pyaerocom.aeroval.modelmaps_helpers import calc_contour_json, CONTOUR, OVERLAY
+from pyaerocom.aeroval.modelmaps_helpers import (
+    calc_contour_json,
+    plot_overlay_pixel_maps,
+    CONTOUR,
+    OVERLAY,
+)
 from pyaerocom.aeroval.varinfo_web import VarinfoWeb
 from pyaerocom.exceptions import (
     DataCoverageError,
@@ -107,7 +112,9 @@ def _run_model(self, model_name: str, var_list):
                     or OVERLAY in self.cfg.modelmaps_opts.plot_types.get(model_name, False)
                 ):
                     # create overlay (pixel) plots
-                    _files = self._process_overlay_map_var(model_name, var)
+                    _files = self._process_overlay_map_var(
+                        model_name, var, self.reanalyse_existing
+                    )
 
             except ModelVarNotAvailable as ex:
                 logger.warning(f"{ex}")
@@ -178,7 +185,7 @@ def _process_contour_map_var(self, model_name, var, reanalyse_existing):
 
         if not reanalyse_existing:
             if os.path.exists(fp_geojson):
-                logger.info(f"Skipping processing of {outname}: data already exists.")
+                logger.info(f"Skipping contour processing of {outname}: data already exists.")
                 return []
 
         maps_freq = TsType(self.cfg.modelmaps_opts.maps_freq)
@@ -209,12 +216,80 @@ def _process_contour_map_var(self, model_name, var, reanalyse_existing):
 
         return fp_geojson
 
-    def _process_overlay_map_var(self, model_name, var):
+    def _process_overlay_map_var(self, model_name, var, reanalyse_existing):
         """Process overlay map (pixels) for either model or obserations
         argument model_name is a misnomer because this can also be applied to observation networks
 
         Args:
             model_name (str): name of model or obs to make overlay pixel maps of
             var (str): variable name
         """
+        breakpoint()
+
+        try:
+            data = self.read_gridded_obsdata(model_name, var)
+        except EntryNotAvailable:
+            try:
+                data = self.read_model_data(model_name, var)
+            except Exception as e:
+                raise ModelVarNotAvailable(
+                    f"Cannot read data for model {model_name} (variable {var}): {e}"
+                )
+
+        var_ranges_defaults = self.cfg.var_scale_colmap
+
+        if var in var_ranges_defaults.keys():
+            cmapinfo = var_ranges_defaults[var]
+            varinfo = VarinfoWeb(var, cmap=cmapinfo["colmap"], cmap_bins=cmapinfo["scale"])
+        else:
+            cmapinfo = var_ranges_defaults["default"]
+            varinfo = VarinfoWeb(var, cmap=cmapinfo["colmap"], cmap_bins=cmapinfo["scale"])
+
+        data = self._check_dimensions(data)
+
+        outdir = self.cfg.path_manager.get_json_output_dirs()["contour/overlay"]
+        outname = f"{var}_{model_name}"
+
+        fp_overlay = os.path.join(
+            outdir, f"{outname}.{self.cfg.modelmaps_opts.overlay_save_format}"
+        )
+
+        if not reanalyse_existing:
+            if os.path.exists(fp_overlay):
+                logger.info(f"Skipping overlay processing of {outname}: data already exists.")
+                return []
+
+        maps_freq = TsType(self.cfg.modelmaps_opts.maps_freq)
+
+        if maps_freq == "coarsest":  # TODO: Implement this in terms of a TsType object. #1267
+            freq = min(TsType(fq) for fq in self.cfg.time_cfg.freqs)
+            freq = min(freq, self.cfg.time_cfg.main_freq)
+        else:
+            freq = maps_freq
+        tst = TsType(data.ts_type)
+
+        if tst < freq:
+            raise TemporalResolutionError(f"need {freq} or higher, got{tst}")
+        elif tst > freq:
+            data = data.resample_time(str(freq))
+
+        data.check_unit()
+
+        # https://github.com/matplotlib/matplotlib/issues/23319
+        overlay_plot = plot_overlay_pixel_maps(
+            data, cmap=varinfo.cmap, cmap_bins=varinfo.cmap_bins, outpath=fp_overlay
+        )
+
         raise NotImplementedError
+
+        # LB: Below will require some coordination with Thorbjørn
+        with self.avdb.lock():
+            self.avdb.put_map_overlay(
+                overlay_plot,
+                self.exp_output.proj_id,
+                self.exp_output.exp_id,
+                var,
+                model_name,
+            )
+
+        return fp_overlay

pyaerocom/aeroval/modelmaps_helpers.py

@@ -5,6 +5,7 @@
 from matplotlib.axes import Axes
 from matplotlib.colors import ListedColormap, to_hex
 from seaborn import color_palette
+import io
 
 try:
     from geojsoncontour import contourf_to_geojson
@@ -103,8 +104,33 @@ def calc_contour_json(data, cmap, cmap_bins):
     return geojson
 
 
-def plot_pixel_maps(
-    data,
-    cmap,
-):
-    pass
+def plot_overlay_pixel_maps(data, cmap, cmap_bins, outpath):
+    fig, axis = plt.subplots(
+        1,
+        1,
+        subplot_kw=dict(projection=ccrs.Mercator()),
+        figsize=(8, 8),
+    )
+    # LB: see if we can use xarray plotting or need to modify somewhow
+    data.plot(
+        ax=axis,
+        transform=ccrs.PlateCarree(),
+        add_colorbar=False,
+        add_labels=False,
+        vmin=0,
+        vmax=120,
+        cmap=cmap,
+    )
+
+    with io.BytesIO() as buffer:  # use buffer memory
+        plt.savefig(
+            outpath,
+            bbox_inches="tight",
+            transparent=True,
+        )
+        buffer.seek(0)
+        image = buffer.getvalue()
+
+    plt.close()
+
+    return image

pyaerocom/aeroval/setup_classes.py

@@ -79,6 +79,7 @@ class OutputPaths(BaseModel):
         "hm/ts",
         "contour",
         "profiles",
+        "contour/overlay",
     ]
     avdb_resource: Path | str | None = None