DuraMAT pref maps

Author: tobin-ford

pvdeg/geospatial.py

@@ -970,6 +970,7 @@ def plot_sparse_analysis(
     method="nearest", 
     resolution:complex=100j,
     figsize:tuple=(10,8),
+    show_plot:bool=False,
 ) -> None:
     """
     Plot the output of a sparse geospatial analysis using interpolation.
@@ -999,7 +1000,7 @@ def plot_sparse_analysis(
     )
 
     fig = plt.figure(figsize=figsize)
-    ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal(), frameon=False)
+    ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal(), frameon=False) # these should be the same ccrs
     ax.patch.set_visible(False)
 
     extent = [lon.min(), lon.max(), lat.min(), lat.max()]
@@ -1009,7 +1010,7 @@ def plot_sparse_analysis(
         extent=extent,
         origin="lower",
         cmap="viridis",
-        transform=ccrs.PlateCarree(),
+        transform=ccrs.PlateCarree(), # why are ccrs different
     )
 
     shapename = "admin_1_states_provinces_lakes"
@@ -1030,6 +1031,70 @@ def plot_sparse_analysis(
     plt.title(f"Interpolated Sparse Analysis, {data_var}")
     plt.xlabel("Longitude")
     plt.ylabel("Latitude")
-    plt.show()
+    
+    if show_plot:
+        plt.show()
+
+    return fig, ax
+
+def plot_sparse_analysis_land(
+    result: xr.Dataset, 
+    data_var: str, 
+    method="nearest", 
+    resolution:complex=100j,
+    figsize:tuple=(10,8),
+    show_plot:bool=False,
+    proj=ccrs.PlateCarree(),
+):
+
+    import matplotlib.path as mpath
+    from cartopy.mpl.patch import geos_to_path
+
+    grid_values, lat, lon = interpolate_analysis(
+        result=result, data_var=data_var, method=method, resolution=resolution
+    )
+
+    fig = plt.figure(figsize=figsize)
+    ax = fig.add_axes([0, 0, 1, 1], projection=proj, frameon=False)
+    ax.patch.set_visible(False)
+
+    extent = [lon.min(), lon.max(), lat.min(), lat.max()]
+    ax.set_extent(extent, crs=proj)
+
+    mesh = ax.pcolormesh(lon, lat, grid_values, transform=proj, cmap='viridis')
+
+    land_path = geos_to_path(list(cfeature.LAND.geometries()))
+    land_path = mpath.Path.make_compound_path(*land_path)
+    plate_carre_data_transform = proj._as_mpl_transform(ax)
+    mesh.set_clip_path(land_path, plate_carre_data_transform)
+
+    shapename = "admin_1_states_provinces_lakes"
+    states_shp = shpreader.natural_earth(
+        resolution="110m", category="cultural", name=shapename
+    )
+
+    ax.add_geometries(
+        shpreader.Reader(states_shp).geometries(),
+        proj,
+        facecolor="none",
+        edgecolor="black",
+        linestyle=':'
+    )
+
+    cbar = plt.colorbar(mesh, ax=ax, orientation="vertical", fraction=0.02, pad=0.04)
+    cbar.set_label("Value")
+
+    utilities._add_cartopy_features(
+        ax=ax, 
+        features = [
+            cfeature.BORDERS,
+            cfeature.COASTLINE,
+            cfeature.LAND,
+            cfeature.OCEAN,
+        ],
+    )
+
+    if show_plot:
+        plt.show()
 
     return fig, ax

pvdeg/scenario.py

@@ -1767,7 +1767,7 @@ def plot_coords(
         coord_2: Optional[tuple[float]] = None,
         coords: Optional[np.ndarray[float]] = None,
         size: Union[int, float] = 1,
-    ) -> None:
+    ):
         """
         Plot lat-long coordinate pairs on blank map. Quickly view
         geospatial datapoints before your analysis.
@@ -1789,6 +1789,11 @@ def plot_coords(
         size : float
             matplotlib scatter point size. Without any downsampling NSRDB
             points will siginficantly overlap.
+
+        Returns:
+        --------
+        fig, ax
+            matplotlib figure and axis
         """
         fig = plt.figure(figsize=(15, 10))
         ax = plt.axes(projection=ccrs.PlateCarree())
@@ -1811,6 +1816,9 @@ def plot_coords(
         plt.title(f"Coordinate Pairs from '{self.name}' Meta Data")
         plt.show()
 
+        return fig, ax
+
+
     def plot_meta_classification(
         self,
         col_name: str = None,
@@ -1844,6 +1852,11 @@ def plot_meta_classification(
         size : float
             matplotlib scatter point size. Without any downsampling NSRDB
             points will siginficantly overlap.
+
+        Returns:
+        --------
+        fig, ax
+            matplotlib figure and axis
         """
         if not col_name:
             raise ValueError("col_name cannot be none")
@@ -1892,6 +1905,8 @@ def plot_meta_classification(
         plt.legend()
         plt.show()
 
+        return fig, ax
+
     def plot_world(
         self,
         data_variable: str,

pvdeg/utilities.py

@@ -11,7 +11,7 @@
 from collections import OrderedDict
 import xarray as xr
 from subprocess import run
-
+import cartopy.feature as cfeature
 
 def gid_downsampling(meta, n):
     """
@@ -1012,7 +1012,7 @@ def _find_bbox_corners(coord_1=None, coord_2=None, coords=None):
         min and max latitude and longitudes. Minimum latitude at lats[0].
         Maximum latitude at lats[1]. Same pattern for longs.
     """
-    if coord_1 and coord_2:
+    if coord_1 is not None and coord_2 is not None:
         lats = [coord_1[0], coord_2[0]]
         longs = [coord_1[1], coord_2[1]]
     elif coords.any():
@@ -1043,20 +1043,21 @@ def _plot_bbox_corners(ax, coord_1=None, coord_2=None, coords=None):
     return
 
 
-def _add_cartopy_features(ax):
+def _add_cartopy_features(
+        ax, 
+        features = [
+            cfeature.BORDERS,
+            cfeature.COASTLINE,
+            cfeature.LAND,
+            cfeature.OCEAN,
+            cfeature.LAKES,
+            cfeature.RIVERS,
+        ],
+    ):
     """
     Add cartopy features to an existing matplotlib.pyplot axis.
     """
-    import cartopy.feature as cfeature
-
-    features = [
-        cfeature.BORDERS,
-        cfeature.COASTLINE,
-        cfeature.LAND,
-        cfeature.OCEAN,
-        cfeature.LAKES,
-        cfeature.RIVERS,
-    ]
+   
 
     for i in features:
         if i == cfeature.BORDERS: