fix plotting from #95

HISTORY.md

@@ -60,7 +60,7 @@ Add the required tools to analyze time series and their properties
 *Next release will use CDO tools instead of xMip tools to do pre-processing for reliability reasons*
 
 * Lon and lat coords for clustering by @JoranAngevaare in https://github.com/JoranAngevaare/optim_esm_tools/pull/79
-* Queriable area field by @JoranAngevaare in https://github.com/JoranAngevaare/optim_esm_tools/pull/80
+* Queryable area field by @JoranAngevaare in https://github.com/JoranAngevaare/optim_esm_tools/pull/80
 
 
 **Full Changelog**: https://github.com/JoranAngevaare/optim_esm_tools/compare/v0.4.0...v0.5.0

optim_esm_tools/analyze/region_finding.py

@@ -259,12 +259,18 @@ def _plot_masks(self, masks, ax=None, legend=True):
             oet.plotting.plot.setup_map()
             ax = plt.gca()
         for i, (label, xy) in enumerate(zip(self.labels, points.values())):
-            ax.scatter(*xy, marker='oxv^'[i], label=f'Maximum {label}')
+            ax.scatter(
+                *xy,
+                marker='oxv^'[i],
+                label=f'Maximum {label}',
+                transform=oet.plotting.plot.get_cartopy_transform(),
+            )
         if legend:
             ax.legend(**oet.utils.legend_kw())
         plt.suptitle(self.title, y=0.95)
-        plt.ylim(-90, 90)
-        plt.xlim(-180, 180)
+        _xlim, _ylim = oet.plotting.plot.get_xy_lim_for_projection()
+        plt.xlim(*_xlim)
+        plt.ylim(*_ylim)
 
     def _mask_to_coord(self, mask_2d):
         arg_mask = np.argwhere(mask_2d)[0]
@@ -453,14 +459,25 @@ def _plot_masks(
 
         ds_dummy['area_square'] = (ds_dummy['cell_area'].dims, all_masks / 1e6)
 
-        ds_dummy['area_square'].plot(cbar_kwargs=mask_cbar_kw, vmin=0, extend='neither')
+        ds_dummy['area_square'].plot(
+            cbar_kwargs=mask_cbar_kw,
+            vmin=0,
+            extend='neither',
+            transform=oet.plotting.plot.get_cartopy_transform(),
+        )
         plt.title('')
         if scatter_medians:
             if cluster_kw is None:
                 cluster_kw = dict()
             for m_i, cluster in enumerate(clusters):
                 lat, lon = np.median(cluster, axis=0)
-                ax.scatter(lon, lat, label=f'cluster {m_i}', **cluster_kw)
+                ax.scatter(
+                    lon,
+                    lat,
+                    label=f'cluster {m_i}',
+                    **cluster_kw,
+                    transform=oet.plotting.plot.get_cartopy_transform(),
+                )
             if legend:
                 plt.legend(**oet.utils.legend_kw())
         plt.suptitle(f'Clusters {self.title}', y=0.97 if len(masks) < 4 else 0.99)

optim_esm_tools/optim_esm_conf.ini

@@ -20,7 +20,6 @@ lon_lat_dim = lon,lat
 
 # If any of these names are in the dataset, remove them as they break pre-processing and are calculated for the regridded file anyway
 remove_vars = area cell_area GEOLON GEOLAT
-cartopy_projection = Robinson
 
 # In the clustering method, assume points closer than clustering_fudge_factor * max_distance belong
 # to the same cluster. See analyze.clustering._infer_max_step_size for more information
@@ -97,3 +96,7 @@ sos = Salinity
 tas=None,None
 tas_detrend=-6,6
 siconc=-5,105
+
+[cartopy]
+projection = Robinson
+transform = PlateCarree

optim_esm_tools/plotting/map_maker.py

@@ -43,7 +43,7 @@ def set_kw(self):
             title=dict(fontsize=12),
             gridspec=dict(hspace=0.3),
             cbar=dict(orientation='horizontal', extend='both'),
-            plot=dict(transform=get_cartopy_projection()),
+            plot=dict(transform=get_cartopy_transform()),
             subplot=dict(projection=get_cartopy_projection()),
         )
 
@@ -172,8 +172,10 @@ def plot_i(self, label, ax=None, coastlines=True, **kw):
 
         plt_ax = prop.plot(**kw)
 
-        plt.xlim(-180, 180)
-        plt.ylim(-90, 90)
+        _xlim, _ylim = oet.plotting.plot.get_xy_lim_for_projection()
+        plt.xlim(*_xlim)
+        plt.ylim(*_ylim)
+
         description = self.conditions[label].long_description
         ax.set_title(label.upper() + '\n' + description, **self.kw['title'])
         gl = ax.gridlines(draw_labels=True)
@@ -463,14 +465,15 @@ def summarize_mask(
                 **mm_sel.kw.get('cbar', {}),
                 **dict(extend='neither', label='Sum of area [km$^2$]'),
             },
+            transform=oet.plotting.plot.get_cartopy_transform(),
         )
         ax.coastlines()
         exponent = int(np.log10(tot_area))
 
         plt.title(f'Area ${tot_area/(10**exponent):.1f}\\times10^{exponent}$ km$^2$')
         gl = ax.gridlines(draw_labels=True)
         gl.top_labels = False
-        gl.right_labels = False
+        # gl.right_labels = False
     else:
         ax = fig.add_subplot(1, 2, 2, projection=get_cartopy_projection())
         mm_sel.plot_i(label=plot, ax=ax, coastlines=True)

optim_esm_tools/plotting/plot.py

@@ -1,6 +1,7 @@
 import matplotlib.pyplot as plt
 import optim_esm_tools as oet
-from optim_esm_tools.config import config
+from optim_esm_tools.config import config, get_logger
+import typing as ty
 
 
 def setup_map(
@@ -64,10 +65,45 @@ def get_unit(ds, var):
     return ds[var].attrs.get('units', '?').replace('%', '\%')
 
 
-def get_cartopy_projection(projection=None, **projection_kwargs):
+def get_cartopy_projection(projection=None, _field='projection', **projection_kwargs):
     import cartopy.crs as ccrs
 
-    projection = projection or config['analyze']['cartopy_projection']
+    projection = projection or config['cartopy'][_field]
     if not hasattr(ccrs, projection):
         raise ValueError(f'Invalid projection {projection}')
     return getattr(ccrs, projection)(**projection_kwargs)
+
+
+def get_cartopy_transform(projection=None, **projection_kwargs):
+    return get_cartopy_projection(
+        projection=projection, _field='transform', **projection_kwargs
+    )
+
+
+def get_xy_lim_for_projection(
+    projection=None,
+) -> ty.Tuple[ty.Tuple[float, float], ty.Tuple[float, float]]:
+    projection = projection or config['cartopy']['projection']
+    lims = dict(
+        Robinson=(
+            (-17005833.33052523, 17005833.33052523),
+            (-8625154.6651, 8625154.6651),
+        ),
+        EqualEarth=(
+            (-17243959.06221695, 17243959.06221695),
+            (-8392927.59846645, 8392927.598466456),
+        ),
+        Mollweide=(
+            (-18040095.696147293, 18040095.696147293),
+            (-9020047.848073646, 9020047.848073646),
+        ),
+        PlateCarree=((0, 360), (-90, 90)),
+    )
+    if projection not in lims:
+        get_logger().warning(
+            f'No hardcoded x/y lims for {projection}, might yield odd figures.'
+        )
+    return lims.get(
+        projection,
+        ((0, 360), (-90, 90)),
+    )