add plot_stations including test

examples/KWK_getcheckdata.py

@@ -135,41 +135,11 @@
 
 ### PLOT SELECTION OF AVAILABLE STATIONS ON MAP
 if plot_stations:
-    crs = 28992
-    locs_meas_ts, locs_meas_ext, _ = kw.data_retrieve.retrieve_catalog(crs=crs)
-    locs_ts_sel = locs_meas_ts.loc[station_list]
-    locs_ext_sel = locs_meas_ext.loc[locs_meas_ext.index.isin(station_list)]
+    station_list_map = station_list.copy()
+    if "NORTHCMRT" in station_list_map:
+        northcmrt_idx = station_list_map.index("NORTHCMRT")
+        station_list_map.pop(northcmrt_idx)
 
-    fig_map,ax_map = plt.subplots(figsize=(8,8))
-    ax_map.plot(locs_ts_sel['X'], locs_ts_sel['Y'],'xk', label="timeseries")
-    ax_map.plot(locs_ext_sel['X'], locs_ext_sel['Y'],'xr', label="extremes")
-    ax_map.legend()
-
-    """
-    for iR, row in locs_ts_sel.iterrows():
-        ax_map.text(row['X'],row['Y'],row.name)
-    """
-    ax_map.set_xlim(-50000,300000) # RD
-    ax_map.set_ylim(350000,850000) # RD
-    ax_map.set_title('stations with timeseries/extremes data')
-    ax_map.set_aspect('equal')
-    ax_map.set_xlabel(f'X (EPSG:{crs})')
-    ax_map.set_ylabel(f'Y (EPSG:{crs})')
-    ax_map.grid(alpha=0.5)
-
-    # optionally add basemap/coastlines
-    try:
-        import dfm_tools as dfmt # pip install dfm_tools
-        dfmt.plot_coastlines(ax=ax_map, crs=crs)
-        dfmt.plot_borders(ax=ax_map, crs=crs)
-    except ModuleNotFoundError:
-        try:
-            import contextily as ctx # pip install contextily
-            ctx.add_basemap(ax_map, source=ctx.providers.Esri.WorldImagery, crs=crs, attribution=False)
-        except ModuleNotFoundError:
-            pass
-
-    fig_map.tight_layout()
-    fig_map.savefig(os.path.join(dir_base,'stations_map.png'), dpi=200)
-
+    fig, ax = kw.plot_stations(station_list=station_list_map, crs=28992, add_labels=False)
+    fig.savefig(os.path.join(dir_base,'stations_map.png'), dpi=200)
 

kenmerkendewaarden/data_analysis.py

@@ -6,13 +6,16 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-from kenmerkendewaarden.data_retrieve import read_measurements, xarray_to_hatyan
+from kenmerkendewaarden.data_retrieve import (read_measurements,
+                                              xarray_to_hatyan,
+                                              retrieve_catalog)
 import hatyan
 import logging
 
 __all__ = [
     "df_amount_pcolormesh",
     "plot_measurements",
+    "plot_stations",
     "derive_statistics",
     ]
 
@@ -80,6 +83,45 @@ def plot_measurements(df, df_ext=None):
     return fig, (ax1,ax2)
 
 
+def plot_stations(station_list, crs=None, add_labels=False):
+    locs_meas_ts_all, locs_meas_ext_all, _ = retrieve_catalog(crs=crs)
+    locs_ts = locs_meas_ts_all.loc[locs_meas_ts_all.index.isin(station_list)]
+    locs_ext = locs_meas_ext_all.loc[locs_meas_ext_all.index.isin(station_list)]
+    if crs is None:
+        crs = int(locs_ts["Coordinatenstelsel"].iloc[0])
+
+    fig, ax = plt.subplots(figsize=(8,8))
+    ax.plot(locs_ts['X'], locs_ts['Y'],'xk', label="timeseries")
+    ax.plot(locs_ext['X'], locs_ext['Y'],'+r', label="extremes")
+    ax.legend()
+
+    ax.set_title('stations with timeseries/extremes data')
+    ax.set_aspect('equal')
+    ax.set_xlabel(f'X (EPSG:{crs})')
+    ax.set_ylabel(f'Y (EPSG:{crs})')
+    ax.grid(alpha=0.5)
+
+    # optionally add basemap/coastlines
+    try:
+        import dfm_tools as dfmt # pip install dfm_tools
+        dfmt.plot_coastlines(ax=ax, crs=crs)
+        dfmt.plot_borders(ax=ax, crs=crs)
+    except ModuleNotFoundError:
+        try:
+            import contextily as ctx # pip install contextily
+            ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=crs, attribution=False)
+        except ModuleNotFoundError:
+            pass
+
+    fig.tight_layout()
+
+    if add_labels:
+        for iR, row in locs_ts.iterrows():
+            ax.text(row['X'],row['Y'],row.name)
+
+    return fig, ax
+
+
 def get_flat_meta_from_dataset(ds):
     list_relevantmetadata = ['WaarnemingMetadata.StatuswaardeLijst', 
                              'WaarnemingMetadata.KwaliteitswaardecodeLijst', 

tests/test_data_analysis.py

@@ -60,3 +60,8 @@ def test_df_amount_pcolormesh(dir_meas_amount, extremes):
 @pytest.mark.unittest
 def test_plot_measurements(df_meas_2010, df_ext_2010):
     kw.plot_measurements(df=df_meas_2010, df_ext=df_ext_2010)
+
+
+@pytest.mark.unittest
+def test_plot_stations():
+    kw.plot_stations(station_list=["HOEKVHLD"], crs=None, add_labels=True)