From 3ad351a228f8e2c6b83f9694ca8a03743e688284 Mon Sep 17 00:00:00 2001 From: Christian Chwala Date: Thu, 12 Aug 2021 10:26:07 +0200 Subject: [PATCH] Cleanup for v0.3.3 (#90) * bump version to v0.3.3. and minor changes * remove unused imports * Reformatted ITU table for better readability and added `# fmt: off` so that black formatting keeps this * Black formatting * update whatsnew with PR since last version --- docs/whats-new.rst | 9 +- .../anomaly_detection/cnn_ano_detection.py | 44 +++--- pycomlink/processing/baseline.py | 6 +- pycomlink/processing/k_R_relation.py | 125 ++---------------- pycomlink/processing/wet_antenna.py | 3 +- pycomlink/processing/xarray_wrapper.py | 1 - pycomlink/spatial/interpolator.py | 18 +-- pycomlink/tests/test_anomaly_detection.py | 47 ++++--- pycomlink/tests/test_cmlh5_to_xarray.py | 7 +- setup.py | 8 +- 10 files changed, 81 insertions(+), 187 deletions(-) diff --git a/docs/whats-new.rst b/docs/whats-new.rst index 7bb5fe2..b6c5bce 100644 --- a/docs/whats-new.rst +++ b/docs/whats-new.rst @@ -2,13 +2,14 @@ What's New ********************** -v0.3.x +v0.3.3 ------ Enhancements ~~~~~~~~~~~~ * Added xarray-wrapper for WAA Leijnse and updated WAA example notebook (by cchwala in PR #82) +* Add CNN-based anomaly detection for CML data (by Glawion in PR#87) * xarray wrapper now uses `xr.apply_ufunc` to apply processing functions along time dimension, instead of looping over the `channel_id` dimension. This should be a lot more flexible. (by cchwala in PR #89) @@ -17,6 +18,12 @@ Bug fixes ~~~~~~~~~ * Fixed problem with xarray_wrapper for calc_R_from_A (by cchwala in PR #89) +Maintenance +~~~~~~~~~~~~ +* Move CI from Travis to Github Actions (by maxmargraf in PR #85) +* Add readthedocs and zenodo badge to README (by maxmargraaf in PR #85) + + v0.3.2 ------ diff --git a/pycomlink/processing/anomaly_detection/cnn_ano_detection.py b/pycomlink/processing/anomaly_detection/cnn_ano_detection.py index 6629db5..a66242b 100644 --- a/pycomlink/processing/anomaly_detection/cnn_ano_detection.py +++ b/pycomlink/processing/anomaly_detection/cnn_ano_detection.py @@ -6,7 +6,6 @@ import pkg_resources - # Limit GPU memory usage to avoid processes to run out of memory. # For a list of processes blocking GPU memory on an nvidia GPU type 'nvidia-smi' in the terminal. config = tf.compat.v1.ConfigProto() @@ -15,7 +14,6 @@ set_session(tf.compat.v1.Session(config=config)) - def get_model_file_path(): return pkg_resources.resource_filename( "pycomlink", "/processing/anomaly_detection/cnn_model" @@ -24,11 +22,10 @@ def get_model_file_path(): # load model modelh5_fn = str(get_model_file_path() + "/model_anomaly_detection_60.h5") -#modelh5_fn = '/pd/home/glawion-l/pycomlink-1/pycomlink/processing/anomaly_detection/cnn_model/model_anomaly_detection_60.h5' +# modelh5_fn = '/pd/home/glawion-l/pycomlink-1/pycomlink/processing/anomaly_detection/cnn_model/model_anomaly_detection_60.h5' model = tf.keras.models.load_model(modelh5_fn, compile=False) - def _rolling_window(a, window): shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) strides = a.strides + (a.strides[-1],) @@ -36,12 +33,13 @@ def _rolling_window(a, window): a, shape=shape, strides=strides, writeable=False ) + def cnn_anomaly_detection( trsl_channel_1, trsl_channel_2, batch_size=100, verbose=0, - ): +): """ Anomaly detection using the CNN based on channel 1 and channel 2 of a CML @@ -68,10 +66,10 @@ def cnn_anomaly_detection( References ---------- - .. [1] Polz, J., Schmidt, L., Glawion, L., Graf, M., Werner, C., Chwala, C., Mollenhauer, H., Rebmann, C., Kunstmann, H., and Bumberger, J.: - Supervised and unsupervised machine-learning for automated quality control of environmental sensor data, EGU General Assembly 2021, online, - 19–30 Apr 2021, EGU21-14485, - https://doi.org/10.5194/egusphere-egu21-14485, 2021. + .. [1] Polz, J., Schmidt, L., Glawion, L., Graf, M., Werner, C., Chwala, C., Mollenhauer, H., Rebmann, C., Kunstmann, H., and Bumberger, J.: + Supervised and unsupervised machine-learning for automated quality control of environmental sensor data, EGU General Assembly 2021, online, + 19–30 Apr 2021, EGU21-14485, + https://doi.org/10.5194/egusphere-egu21-14485, 2021. """ df = pd.DataFrame() @@ -84,18 +82,16 @@ def cnn_anomaly_detection( df = df.fillna(value=-9999) - x_fts =np.moveaxis( - np.array( - [ - _rolling_window(df["trsl1"].values, 60), - _rolling_window(df["trsl2"].values, 60), - ] - ), - source=0, - destination=-1, - ) - - + x_fts = np.moveaxis( + np.array( + [ + _rolling_window(df["trsl1"].values, 60), + _rolling_window(df["trsl2"].values, 60), + ] + ), + source=0, + destination=-1, + ) cnn_pred = np.ravel(model.predict(x_fts, batch_size=batch_size, verbose=verbose)) @@ -106,8 +102,6 @@ def cnn_anomaly_detection( ): cnn_pred[i] = np.nan - df["prediction"] = np.concatenate( - (np.repeat(np.nan, 59), cnn_pred), axis=0 - ) + df["prediction"] = np.concatenate((np.repeat(np.nan, 59), cnn_pred), axis=0) - return df.prediction.values.reshape(len(trsl_channel_1)) \ No newline at end of file + return df.prediction.values.reshape(len(trsl_channel_1)) diff --git a/pycomlink/processing/baseline.py b/pycomlink/processing/baseline.py index 3da5e98..0a2e5f1 100644 --- a/pycomlink/processing/baseline.py +++ b/pycomlink/processing/baseline.py @@ -55,9 +55,9 @@ def _numba_baseline_constant(trsl, wet, n_average_last_dry): for i in range(n_average_last_dry, len(trsl)): if np.isnan(wet[i]): baseline[i] = np.NaN - elif wet[i] & ~wet[i-1]: - baseline[i] = np.mean(baseline[(i-n_average_last_dry) : i]) - elif wet[i] & wet[i-1]: + elif wet[i] & ~wet[i - 1]: + baseline[i] = np.mean(baseline[(i - n_average_last_dry) : i]) + elif wet[i] & wet[i - 1]: baseline[i] = baseline[i - 1] else: baseline[i] = trsl[i] diff --git a/pycomlink/processing/k_R_relation.py b/pycomlink/processing/k_R_relation.py index 54c0572..b589d09 100644 --- a/pycomlink/processing/k_R_relation.py +++ b/pycomlink/processing/k_R_relation.py @@ -8,6 +8,7 @@ # Functions for k-R power law calculations # ############################################ + @xarray_apply_along_time_dim() def calc_R_from_A(A, L_km, f_GHz=None, a=None, b=None, pol="H", R_min=0.1): """Calculate rain rate from attenuation using the A-R Relationship @@ -86,7 +87,7 @@ def calc_R_from_A_min_max( ------- float or iterable of float Rain rate - + Note ---- Based on: "Empirical Study of the Quantization Bias Effects in @@ -170,122 +171,14 @@ def a_b(f_GHz, pol, approx_type="ITU"): return a, b +# fmt: off ITU_table = np.array( [ - [ - 1.000e0, - 2.000e0, - 4.000e0, - 6.000e0, - 7.000e0, - 8.000e0, - 1.000e1, - 1.200e1, - 1.500e1, - 2.000e1, - 2.500e1, - 3.000e1, - 3.500e1, - 4.000e1, - 4.500e1, - 5.000e1, - 6.000e1, - 7.000e1, - 8.000e1, - 9.000e1, - 1.000e2, - ], - [ - 3.870e-5, - 2.000e-4, - 6.000e-4, - 1.800e-3, - 3.000e-3, - 4.500e-3, - 1.010e-2, - 1.880e-2, - 3.670e-2, - 7.510e-2, - 1.240e-1, - 1.870e-1, - 2.630e-1, - 3.500e-1, - 4.420e-1, - 5.360e-1, - 7.070e-1, - 8.510e-1, - 9.750e-1, - 1.060e0, - 1.120e0, - ], - [ - 3.520e-5, - 1.000e-4, - 6.000e-4, - 1.600e-3, - 2.600e-3, - 4.000e-3, - 8.900e-3, - 1.680e-2, - 3.350e-2, - 6.910e-2, - 1.130e-1, - 1.670e-1, - 2.330e-1, - 3.100e-1, - 3.930e-1, - 4.790e-1, - 6.420e-1, - 7.840e-1, - 9.060e-1, - 9.990e-1, - 1.060e0, - ], - [ - 9.120e-1, - 9.630e-1, - 1.121e0, - 1.308e0, - 1.332e0, - 1.327e0, - 1.276e0, - 1.217e0, - 1.154e0, - 1.099e0, - 1.061e0, - 1.021e0, - 9.790e-1, - 9.390e-1, - 9.030e-1, - 8.730e-1, - 8.260e-1, - 7.930e-1, - 7.690e-1, - 7.530e-1, - 7.430e-1, - ], - [ - 8.800e-1, - 9.230e-1, - 1.075e0, - 1.265e0, - 1.312e0, - 1.310e0, - 1.264e0, - 1.200e0, - 1.128e0, - 1.065e0, - 1.030e0, - 1.000e0, - 9.630e-1, - 9.290e-1, - 8.970e-1, - 8.680e-1, - 8.240e-1, - 7.930e-1, - 7.690e-1, - 7.540e-1, - 7.440e-1, - ], + [1.000e0, 2.000e0, 4.000e0, 6.000e0, 7.000e0, 8.000e0, 1.000e1, 1.200e1, 1.500e1, 2.000e1, 2.500e1, 3.000e1, 3.500e1, 4.000e1, 4.500e1, 5.000e1, 6.000e1, 7.000e1, 8.000e1, 9.000e1, 1.000e2], + [3.870e-5, 2.000e-4, 6.000e-4, 1.800e-3, 3.000e-3, 4.500e-3, 1.010e-2, 1.880e-2, 3.670e-2, 7.510e-2, 1.240e-1, 1.870e-1, 2.630e-1, 3.500e-1, 4.420e-1, 5.360e-1, 7.070e-1, 8.510e-1, 9.750e-1, 1.060e0, 1.120e0], + [3.520e-5, 1.000e-4, 6.000e-4, 1.600e-3, 2.600e-3, 4.000e-3, 8.900e-3, 1.680e-2, 3.350e-2, 6.910e-2, 1.130e-1, 1.670e-1, 2.330e-1, 3.100e-1, 3.930e-1, 4.790e-1, 6.420e-1, 7.840e-1, 9.060e-1, 9.990e-1, 1.060e0], + [9.120e-1, 9.630e-1, 1.121e0, 1.308e0, 1.332e0, 1.327e0, 1.276e0, 1.217e0, 1.154e0, 1.099e0, 1.061e0, 1.021e0, 9.790e-1, 9.390e-1, 9.030e-1, 8.730e-1, 8.260e-1, 7.930e-1, 7.690e-1, 7.530e-1, 7.430e-1], + [8.800e-1, 9.230e-1, 1.075e0, 1.265e0, 1.312e0, 1.310e0, 1.264e0, 1.200e0, 1.128e0, 1.065e0, 1.030e0, 1.000e0, 9.630e-1, 9.290e-1, 8.970e-1, 8.680e-1, 8.240e-1, 7.930e-1, 7.690e-1, 7.540e-1, 7.440e-1], ] ) +# fmt: on diff --git a/pycomlink/processing/wet_antenna.py b/pycomlink/processing/wet_antenna.py index aa933eb..5dc30d3 100644 --- a/pycomlink/processing/wet_antenna.py +++ b/pycomlink/processing/wet_antenna.py @@ -109,6 +109,7 @@ def waa_schleiss_2013(rsl, baseline, wet, waa_max, delta_t, tau): return waa + @xarray_apply_along_time_dim() def waa_leijnse_2008_from_A_obs( A_obs, @@ -168,7 +169,7 @@ def waa_leijnse_2008_from_A_obs( """ if np.any(A_obs < 0): - raise ValueError('Negative values for `A_obs` are not allowed') + raise ValueError("Negative values for `A_obs` are not allowed") # Make sure that L_km is not an array or xarray.Dataarray with a size greater # than 1, i.e. it has to be a single scalar values. This is required so that diff --git a/pycomlink/processing/xarray_wrapper.py b/pycomlink/processing/xarray_wrapper.py index 292a4e3..77b9da7 100644 --- a/pycomlink/processing/xarray_wrapper.py +++ b/pycomlink/processing/xarray_wrapper.py @@ -1,7 +1,6 @@ from functools import wraps from collections import OrderedDict import inspect -import numpy as np import xarray as xr diff --git a/pycomlink/spatial/interpolator.py b/pycomlink/spatial/interpolator.py index cbc3bad..290d0c1 100644 --- a/pycomlink/spatial/interpolator.py +++ b/pycomlink/spatial/interpolator.py @@ -17,11 +17,11 @@ class PointsToGridInterpolator(with_metaclass(abc.ABCMeta, object)): - """ PointsToGridInterpolator class docstring """ + """PointsToGridInterpolator class docstring""" @abc.abstractmethod def __init__(self): - """ Pass all configuration parameters for the interpolator here """ + """Pass all configuration parameters for the interpolator here""" return def __call__(self, x, y, z, xgrid=None, ygrid=None, resolution=None): @@ -65,13 +65,13 @@ def __call__(self, x, y, z, xgrid=None, ygrid=None, resolution=None): @abc.abstractmethod def _interpol_func(self, x, y, z, xi, yi): - """ The actual interpolation code goes here """ + """The actual interpolation code goes here""" return class IdwKdtreeInterpolator(PointsToGridInterpolator): def __init__(self, nnear=8, p=2, exclude_nan=True, max_distance=None): - """ A k-d tree based IDW interpolator for points to grid """ + """A k-d tree based IDW interpolator for points to grid""" self.nnear = nnear self.p = p self.exclude_nan = exclude_nan @@ -80,7 +80,7 @@ def __init__(self, nnear=8, p=2, exclude_nan=True, max_distance=None): self.y = None def _interpol_func(self, x, y, z, xi, yi): - """ Do IDW interpolation """ + """Do IDW interpolation""" x = np.asarray(x) y = np.asarray(y) @@ -94,10 +94,10 @@ def _interpol_func(self, x, y, z, xi, yi): self.z = z if np.array_equal(x, self.x) and np.array_equal(y, self.y): - #print('Reusing old `Invdisttree`') + # print('Reusing old `Invdisttree`') idw = self.idw else: - #print('Building new `Invdistree`') + # print('Building new `Invdistree`') idw = Invdisttree(X=list(zip(x, y))) self.idw = idw self.x = x @@ -123,7 +123,7 @@ def __init__( # coordinates_type='euclidean', # Not supported in v1.3.1 backend="C", ): - """ A ordinary kriging interpolator for points to grid""" + """A ordinary kriging interpolator for points to grid""" self.nlags = nlags self.variogram_model = variogram_model @@ -191,7 +191,7 @@ def _parse_grid_kwargs(x_list, y_list, xgrid, ygrid, resolution): def get_lon_lat_list_from_cml_list(cml_list): - """ Extract lats and lons from all CMLs """ + """Extract lats and lons from all CMLs""" lons = np.array([cml.get_center_lon_lat()[0] for cml in cml_list]) lats = np.array([cml.get_center_lon_lat()[1] for cml in cml_list]) diff --git a/pycomlink/tests/test_anomaly_detection.py b/pycomlink/tests/test_anomaly_detection.py index b522c76..90a0998 100644 --- a/pycomlink/tests/test_anomaly_detection.py +++ b/pycomlink/tests/test_anomaly_detection.py @@ -1,39 +1,42 @@ import unittest import numpy as np -from pycomlink.processing.anomaly_detection.cnn_ano_detection import cnn_anomaly_detection +from pycomlink.processing.anomaly_detection.cnn_ano_detection import ( + cnn_anomaly_detection, +) class Testcnnpred(unittest.TestCase): def test_cnnpred(self): - trsl_channel_1 = np.arange(0, 5,0.01).astype(float) - trsl_channel_2 = np.arange(0, 5,0.01).astype(float) + trsl_channel_1 = np.arange(0, 5, 0.01).astype(float) + trsl_channel_2 = np.arange(0, 5, 0.01).astype(float) trsl_channel_1[320] = np.nan pred_raw = cnn_anomaly_detection( - trsl_channel_1, - trsl_channel_2, - batch_size=1, - verbose=0, - ) + trsl_channel_1, + trsl_channel_2, + batch_size=1, + verbose=0, + ) - - assert len(pred_raw) == len(np.arange(0, 5,0.01)) + assert len(pred_raw) == len(np.arange(0, 5, 0.01)) truth_raw = np.array( - [np.nan, - np.nan, - np.nan, - 0.5839 , - 0.5869, - 0.5899, - 0.5931 , - 0.5964, - 0.5997, - 0.6031 ] - ) + [ + np.nan, + np.nan, + np.nan, + 0.5839, + 0.5869, + 0.5899, + 0.5931, + 0.5964, + 0.5997, + 0.6031, + ] + ) np.testing.assert_almost_equal( np.round(pred_raw, decimals=4)[235:245], truth_raw - ) \ No newline at end of file + ) diff --git a/pycomlink/tests/test_cmlh5_to_xarray.py b/pycomlink/tests/test_cmlh5_to_xarray.py index 01a499d..ad5f942 100644 --- a/pycomlink/tests/test_cmlh5_to_xarray.py +++ b/pycomlink/tests/test_cmlh5_to_xarray.py @@ -76,9 +76,7 @@ def test_cmlh5_to_xarray_fun(self): ) np.testing.assert_array_equal( - np.array( - [51.0909, 51.1271, 50.8026, 50.738299999999995] - ), + np.array([51.0909, 51.1271, 50.8026, 50.738299999999995]), np.array( [ cml_xr_file[12].site_a_longitude.values, @@ -86,6 +84,5 @@ def test_cmlh5_to_xarray_fun(self): cml_xr_file[12].site_a_latitude.values, cml_xr_file[12].site_b_latitude.values, ], - ) + ), ) - diff --git a/setup.py b/setup.py index 9d22fc6..26b751c 100644 --- a/setup.py +++ b/setup.py @@ -20,15 +20,15 @@ def read(fname): setup( name = "pycomlink", - version = "0.3.2", + version = "0.3.3", author = "Christian Chwala", author_email = "christian.chwala@kit.edu", - description = ("Python tools for MW link data processing"), + description = ("Python tools for CML (commercial microwave link) data processing"), license = "BSD", - keywords = "microwave links precipitation radar", + keywords = "microwave links precipitation radar cml", url = "https://github.com/pycomlink/pycomlink", download_url = ( - "https://github.com/pycomlink/pycomlink/archive/0.3.2.tar.gz"), + "https://github.com/pycomlink/pycomlink/archive/0.3.3.tar.gz"), packages=find_packages(exclude=['test']), include_package_data=True, long_description=read('README.md'),