changes for mapping

* if global model, just plot extent around input dataset locations, not global map, and put note in title
* can input extent now
* code to shift longitudes is now a utils function and has tests added

ocean_model_skill_assessor/__init__.py

@@ -23,6 +23,7 @@
     compute_root_mean_square_error,
     compute_stats,
 )
+from .utils import shift_longitudes
 
 
 try:

ocean_model_skill_assessor/main.py

@@ -26,15 +26,7 @@
 
 from ocean_model_skill_assessor.plot import map, time_series
 
-from .utils import kwargs_search_from_model
-
-
-try:
-    import cartopy
-
-    CARTOPY_AVAILABLE = True
-except ImportError:  # pragma: no cover
-    CARTOPY_AVAILABLE = False  # pragma: no cover
+from .utils import kwargs_search_from_model, shift_longitudes
 
 
 def make_local_catalog(
@@ -447,16 +439,17 @@ def run(
     with cfx.set_options(custom_criteria=vocab.vocab):
         dam = dsm.cf[key_variable]
 
+    # include matching static mask if present
+    masks = dsm.filter_by_attrs(flag_meanings="land water")
+    if len(masks.data_vars) > 0:
+        mask_name = [mask for mask in masks.data_vars if dsm[mask].encoding["coordinates"] in dam.encoding["coordinates"]][0]
+        dam = xr.merge([dam, dsm[mask_name]])
+    else:
+        # still need it to be a dataset
+        dam = dam.to_dataset()
+
     # shift if 0 to 360
-    if dam.cf["longitude"].max() > 180:
-        lkey = dam.cf["longitude"].name
-        dam = dam.assign_coords(lon=(((dam[lkey] + 180) % 360) - 180))
-        # rotate arrays so that the locations and values are -180 to 180
-        # instead of 0 to 180 to -180 to 0
-        dam = dam.roll(lon=int((dam[lkey] < 0).sum()), roll_coords=True)
-        print(
-            "Longitudes are being shifted because they look like they are not -180 to 180."
-        )
+    dam = shift_longitudes(dam)
 
     # loop over catalogs and sources to pull out lon/lat locations for plot
     maps = []
@@ -577,12 +570,15 @@ def run(
             count += 1
 
     # map of model domain with data locations
-    if CARTOPY_AVAILABLE and len(maps) > 0:
-        figname = omsa.PROJ_DIR(project_name) / "map.png"
-        omsa.plot.map.plot_map(np.asarray(maps), figname, dsm, **kwargs_map)
+    if len(maps) > 0:
+        try:
+            figname = omsa.PROJ_DIR(project_name) / "map.png"
+            omsa.plot.map.plot_map(np.asarray(maps), figname, dsm, **kwargs_map)
+        except ModuleNotFoundError:
+            pass
     else:
         print(
-            "Not plotting map since cartopy is not installed or no datasets to work with."
+            "Not plotting map since no datasets to plot."
         )
     print(
         f"Finished analysis. Find plots and stats summaries in {omsa.PROJ_DIR(project_name)}."

ocean_model_skill_assessor/plot/map.py

@@ -5,15 +5,22 @@
 import pathlib
 
 from pathlib import PurePath
-from typing import Union
+from typing import Optional, Sequence, Union
 
 import matplotlib.patches as mpatches
 import matplotlib.pyplot as plt
 import numpy as np
 
 from xarray import DataArray, Dataset
 
-from ..utils import find_bbox
+from ..utils import find_bbox, shift_longitudes
+
+try:
+    import cartopy
+
+    CARTOPY_AVAILABLE = True
+except ImportError:  # pragma: no cover
+    CARTOPY_AVAILABLE = False  # pragma: no cover
 
 
 def plot_map(
@@ -22,6 +29,7 @@ def plot_map(
     ds: Union[DataArray, Dataset],
     alpha: int = 5,
     dd: int = 2,
+    extent: Optional[Sequence] = None,
 ):
     """Plot and save to file map of model domain and data locations.
 
@@ -33,8 +41,17 @@ def plot_map(
         Map will be saved here.
     ds : Union[DataArray, Dataset]
         Model output.
+    extent: 
+        [min longitude, max longitude, min latitude, max latitude]
     """
-
+
+    ds = shift_longitudes(ds)
+
+    if not CARTOPY_AVAILABLE:
+        raise ModuleNotFoundError(  # pragma: no cover
+            "Cartopy is not available so map will not be plotted."
+        )
+
     import cartopy
 
     pc = cartopy.crs.PlateCarree()
@@ -83,7 +100,22 @@ def plot_map(
         ax.annotate(i, xy=xyproj, xytext=xyproj, color=col_label)
 
     # [min lon, max lon, min lat, max lat]
-    extent = [bbox[0] - 0.1, bbox[2] + 0.1, bbox[1] - 0.1, bbox[3] + 0.1]
-    ax.set_extent(extent, pc)
+    if extent is None:
+        extent_use = [bbox[0] - 0.1, bbox[2] + 0.1, bbox[1] - 0.1, bbox[3] + 0.1]
+
+    # if model is global - based on extent - write that it is global and use smaller extent
+    if np.allclose(bbox[0], -180, atol=2) and np.allclose(bbox[2], 180, atol=2) and np.allclose(bbox[1], -90, atol=2) and np.allclose(bbox[3], 90, atol=2):
+        # explain global model
+        ax.set_title("Only showing part of global model")
+
+        # change delta deg for extent to max(10% of total diff lons/lats, 1 deg)
+        if extent is None:
+            dlon, dlat = 0.1*(min(min_lons) - max(max_lons)), 0.1*(min(min_lats) - max(max_lats))
+            ddlon, ddlat = max(dlon, 5), max(dlat, 2)
+            extent_use = [min(min_lons) - ddlon, max(max_lons) + ddlon,
+                          min(min_lats) - ddlat, max(max_lats) + ddlat]
+
+
+    ax.set_extent(extent_use, pc)
 
     fig.savefig(figname, dpi=100, bbox_inches="tight")

ocean_model_skill_assessor/utils.py

@@ -11,12 +11,40 @@
 import numpy as np
 import shapely.geometry
 import xarray as xr
+from xarray import Dataset, DataArray
 
 from intake.catalog import Catalog
 
 import ocean_model_skill_assessor as omsa
 
 
+def shift_longitudes(dam: Union[DataArray,Dataset]) -> Union[DataArray,Dataset]:
+    """Shift longitudes from 0 to 360 to -180 to 180 if necessary.
+
+    Parameters
+    ----------
+    dam : Union[DataArray,Dataset]
+        Object with model output to check
+
+    Returns
+    -------
+    Union[DataArray,Dataset]
+        Return model output with shifted longitudes, if it was necessary.
+    """
+
+    if dam.cf["longitude"].max() > 180:
+        lkey = dam.cf["longitude"].name
+        nlon = int((dam[lkey] >= 180).sum())  # number of longitudes to roll by
+        dam = dam.assign_coords(lon=(((dam[lkey] + 180) % 360) - 180))
+        # rotate arrays so that the locations and values are -180 to 180
+        # instead of 0 to 180 to -180 to 0
+        dam = dam.roll(lon=nlon, roll_coords=True)
+        print(
+            "Longitudes are being shifted because they look like they are not -180 to 180."
+        )
+    return dam
+
+
 def find_bbox(ds: xr.DataArray, dd: int = 1, alpha: int = 5) -> tuple:
     """Determine bounds and boundary of model.
 

tests/test_utils.py

@@ -125,3 +125,22 @@ def test_find_bbox():
     assert latkey == "lat"
     assert bbox == [0.0, 0.0, 9.0, 9.0]
     assert isinstance(p1, shapely.geometry.polygon.Polygon)
+
+
+def test_shift_longitudes():
+
+    ds = xr.Dataset()
+    ds["lon"] = (
+        "lon",
+        np.linspace(0,360,5)[:-1],
+        {"units": "degrees_east", "standard_name": "longitude"},
+    )
+    assert all(omsa.shift_longitudes(ds).cf["longitude"] == [-180.,  -90.,    0.,   90.])
+
+    ds = xr.Dataset()
+    ds["lon"] = (
+        "lon",
+        np.linspace(-180,180,5)[:-1],
+        {"units": "degrees_east", "standard_name": "longitude"},
+    )
+    assert all(omsa.shift_longitudes(ds).cf["longitude"] == ds.cf["longitude"])