2
\$\begingroup\$

I have some data that represents the area, speed(meters per second), and azimuth(rads) of a polygon. The objective is to determine the mean for each set of parameters. With that mean value, apply a threshold to remove data points that fell outside of the mean deviation threshold. Then recalculate the means for the parameters set.

I am using match statements so you will need python 3.10

data

DATA = [
    [0.045750000000000984, 0.04939999999999974, 0.05244999999999857, 0.054399999999998436, 0.056099999999999456, 0.05375000000000011, 0.047700000000000985, 0.053700000000000594, 0.053249999999999686, 16.81209831859011, 27.983008063348617, 20.301826794339316, 13.421484359287598,
        17.10557180223468, 14.686987606192902, 41.18082589372154, 13.065486605460656, 4.608155337880411, 0.7571951341718915, 0.3457777398771206, 0.46056961460740087, 0.2819901554949584, 0.7205708724903156, -0.00040991561554303056, 0.3860776453927442, 0.09629955035400072, 0.2468722764784534],
    [0.03290000000000123, 0.033550000000000045, 0.03949999999999949, 0.035550000000001615, 0.03725000000000067, 0.03515000000000022, 0.03714999999999999, 0.041350000000001316, 0.03869999999999882, 14.902912040542054, 15.452401255651447, 12.5484544315597, 10.83783196920491,
        13.246734412134916, 12.636336066844674, 19.61093585673077, 8.898327945757917, 26.237157593452878, 0.36589566404637536, 0.5140598344685876, 0.06092353057665044, 0.5433592643289422, 0.3300217799943245, 0.1295695630661193, 0.5813772760582753, 0.97888742111902, 0.2144469684077374],
    [0.15945000000000153, 0.1419000000000014, 0.1346, 0.14314999999999908, 0.14485000000000084, 0.1293000000000014, 0.14264999999999992, 0.1434499999999994, 0.08944999999999861, 39.54982055905117, 67.05598966736439, 14.170464300161013, 25.57103344352307, 22.224605807376317,
        28.901343579758898, 19.05130673155152, 25.61289901736464, 122.31441473369556, 0.5250514668415044, 0.3119136953061388, 0.3818093125132205, 0.2669557624118119, -0.4876722558321518, 0.03814067199937899, -0.07502082862542193, 0.15371372713178985, -1.1026290280269917],
    [0.00340000000000014, 0.0066000000000002515, 0.008499999999999603, 0.007999999999999518, 0.0065500000000002015, 0.008100000000000543, 0.0032500000000004834, 0.018350000000000196, 0.020150000000000542, 22.364276562983136, 20.506557739907425, 6.243036570120132, 19.3458717607394,
        10.456833536986126, 19.65879720239808, 14.307243539019602, 9.091403905263792, 18.47420699835527, 0.8485690702159555, 0.40094252533477115, 0.1157765828210181, 0.6944022345849206, 0.46766736064272035, 1.153657712332056, 0.7600461290670872, 0.4170470750740715, 0.49920935232659946],
    [None, None, 0.003549999999999955, 0.005799999999999788, 0.005799999999999788, 0.006900000000000275, 0.007049999999999682, 0.006800000000000173, 0.005200000000000045, None, None, 8.06110268044255, 13.695219329907925, 0.0,
        12.81872890572318, 6.190861120089245, 3.9861242549543703, 4.403761437912741, None, None, -0.8342338501323102, -0.35990880371032785, 0.0, -2.9708131816772876, -0.12378167714772335, -1.4646911328887278, -2.475146432921943]
]

main.py


import functools
from typing import Iterable, Dict
import pandas as pd
import numpy as np

DATA = ...
INDEX = ['89234', '89467', '89519', '89680', '90104']
#COLUMN = np.repeat(['area', 'MPS', 'azimuth'], 9)
# added MULTI_COLUMN per request
DATE_RANGE = pd.to_datetime(range(9), unit='m', origin=pd.Timestamp('1960-01-01'))

MULTI_COLUMN = pd.MultiIndex.from_tuples((
    (date, param) for param in ('area', 'MPS', 'azimuth')
    for date in DATE_RANGE), names=['validTime', 'parameters'])

def normalize(
        frame_cache: pd.DataFrame, thresh) -> pd.DataFrame:
    """normalize values in the dataframe"""

    frame_cache[['area', 'MPS', 'azimuth']] = np.block(
        list(outliers2nan(frame_cache, thresh=thresh)))

    return frame_cache


def outliers2nan(
        frame_cache: pd.DataFrame, thresh: Dict[str, float]) -> Iterable[np.ndarray]:
    """
    where a parameter is greater than the minimum mean threshold set the value to NaN
    """

    for param, frame in frame_cache.groupby('parameters', axis=1, sort=False):
        mean = mean_stack(frame, kind=param)

        match param:
            case 'azimuth':

                deg_diff = (
                    np.rad2deg(frame_cache['azimuth']) - np.rad2deg(mean))

                yield frame.where(abs((deg_diff + 180) % 360 - 180) < thresh[param])

            case _:

                yield frame.where(abs(frame-mean) < mean*thresh[param])


def mean_stack(
        frame: pd.DataFrame, axis=1, kind=None) -> np.ndarray:
    """DataFrame.mean(1) stacked as a np.ndarray"""

    match kind:
        case 'azimuth':
            return (
                functools.reduce(_mean_azimuth_reduction, frame.to_numpy().T)[:, np.newaxis])

        case _:
            return (
                frame.mean(axis=axis).to_numpy()[:, np.newaxis])


def _mean_azimuth_reduction(col1: np.ndarray, col2: np.ndarray) -> np.ndarray:

    np.nan_to_num(x=col1, nan=col2, copy=False)

    mean_rads = (
        np.arctan2(np.sin(col1)+np.sin(col2), np.cos(col1)+np.cos(col2)))

    return mean_rads


def get_means(df: pd.DataFrame) -> pd.DataFrame:
    """start subroutine"""

    return pd.DataFrame({
        param: mean_stack(frame, kind=param).squeeze()
        for param, frame in df.groupby('parameters', axis=1, sort=False)
    }, index=INDEX)


def start():
    """construct df and display means"""

    bad_df = pd.DataFrame(
        DATA, index=INDEX, columns=MULTI_COLUMN)

    bad_means = get_means(bad_df)

    good_df = normalize(
        bad_df.droplevel('validTime', axis=1).copy(), thresh={'area': .65, 'MPS': .65, 'azimuth': 60.0})

    good_means = get_means(good_df)
    good_df.columns = MULTI_COLUMN
    print(f"""
BAD
{bad_df.stack()}

GOOD
{good_df.stack()}

BAD
{bad_means}

GOOD
{good_means}

DIF
{abs(bad_means-good_means)}
    """)


if __name__ == '__main__':
    start()

result

A good example of the use case is on the last MPS column of index 89519 which reported a speed of 122 meters per second. That value gets omitted from the final return value.

BAD
validTime         1960-01-01 00:00:00  1960-01-01 00:01:00  1960-01-01 00:02:00  1960-01-01 00:03:00  1960-01-01 00:04:00  1960-01-01 00:05:00  1960-01-01 00:06:00  1960-01-01 00:07:00  1960-01-01 00:08:00
      parameters                                                                                                                                                                                             
89234 MPS                   16.812098            27.983008            20.301827            13.421484            17.105572            14.686988            41.180826            13.065487             4.608155
      area                   0.045750             0.049400             0.052450             0.054400             0.056100             0.053750             0.047700             0.053700             0.053250
      azimuth                0.757195             0.345778             0.460570             0.281990             0.720571            -0.000410             0.386078             0.096300             0.246872
89467 MPS                   14.902912            15.452401            12.548454            10.837832            13.246734            12.636336            19.610936             8.898328            26.237158
      area                   0.032900             0.033550             0.039500             0.035550             0.037250             0.035150             0.037150             0.041350             0.038700
      azimuth                0.365896             0.514060             0.060924             0.543359             0.330022             0.129570             0.581377             0.978887             0.214447
89519 MPS                   39.549821            67.055990            14.170464            25.571033            22.224606            28.901344            19.051307            25.612899           122.314415
      area                   0.159450             0.141900             0.134600             0.143150             0.144850             0.129300             0.142650             0.143450             0.089450
      azimuth                0.525051             0.311914             0.381809             0.266956            -0.487672             0.038141            -0.075021             0.153714            -1.102629
89680 MPS                   22.364277            20.506558             6.243037            19.345872            10.456834            19.658797            14.307244             9.091404            18.474207
      area                   0.003400             0.006600             0.008500             0.008000             0.006550             0.008100             0.003250             0.018350             0.020150
      azimuth                0.848569             0.400943             0.115777             0.694402             0.467667             1.153658             0.760046             0.417047             0.499209
90104 MPS                         NaN                  NaN             8.061103            13.695219             0.000000            12.818729             6.190861             3.986124             4.403761
      area                        NaN                  NaN             0.003550             0.005800             0.005800             0.006900             0.007050             0.006800             0.005200
      azimuth                     NaN                  NaN            -0.834234            -0.359909             0.000000            -2.970813            -0.123782            -1.464691            -2.475146

GOOD
               1960-01-01 00:00:00  1960-01-01 00:01:00  1960-01-01 00:02:00  1960-01-01 00:03:00  1960-01-01 00:04:00  1960-01-01 00:05:00  1960-01-01 00:06:00  1960-01-01 00:07:00  1960-01-01 00:08:00
89234 MPS                16.812098            27.983008            20.301827            13.421484            17.105572            14.686988                  NaN            13.065487                  NaN
      area                0.045750             0.049400             0.052450             0.054400             0.056100             0.053750             0.047700             0.053700             0.053250
      azimuth             0.757195             0.345778             0.460570             0.281990             0.720571            -0.000410             0.386078             0.096300             0.246872
89467 MPS                14.902912            15.452401            12.548454            10.837832            13.246734            12.636336            19.610936             8.898328                  NaN
      area                0.032900             0.033550             0.039500             0.035550             0.037250             0.035150             0.037150             0.041350             0.038700
      azimuth             0.365896             0.514060             0.060924             0.543359             0.330022             0.129570             0.581377             0.978887             0.214447
89519 MPS                39.549821                  NaN                  NaN            25.571033            22.224606            28.901344            19.051307            25.612899                  NaN
      area                0.159450             0.141900             0.134600             0.143150             0.144850             0.129300             0.142650             0.143450             0.089450
      azimuth                  NaN             0.311914             0.381809             0.266956            -0.487672             0.038141            -0.075021             0.153714            -1.102629
89680 MPS                22.364277            20.506558             6.243037            19.345872            10.456834            19.658797            14.307244             9.091404            18.474207
      area                0.003400             0.006600             0.008500             0.008000             0.006550             0.008100             0.003250                  NaN                  NaN
      azimuth             0.848569             0.400943             0.115777             0.694402             0.467667             1.153658             0.760046             0.417047             0.499209
90104 MPS                      NaN                  NaN             8.061103                  NaN                  NaN                  NaN             6.190861             3.986124             4.403761
      area                     NaN                  NaN             0.003550             0.005800             0.005800             0.006900             0.007050             0.006800             0.005200
      azimuth                  NaN                  NaN            -0.834234                  NaN                  NaN                  NaN                  NaN            -1.464691            -2.475146

BAD
           area        MPS   azimuth
89234  0.051833  18.796161  0.230576
89467  0.036789  14.930121  0.455432
89519  0.136533  40.494653 -0.524696
89680  0.009211  15.605359  0.552226
90104  0.005871   7.022257 -1.823553

GOOD
           area        MPS   azimuth
89234  0.051833  17.625209  0.230576
89467  0.036789  13.516742  0.455432
89519  0.136533  26.818502 -0.525529
89680  0.006343  15.605359  0.552226
90104  0.005871   5.660462 -1.969919

DIF
           area        MPS   azimuth
89234  0.000000   1.170951  0.000000
89467  0.000000   1.413380  0.000000
89519  0.000000  13.676152  0.000833
89680  0.002868   0.000000  0.000000
90104  0.000000   1.361794  0.146366
\$\endgroup\$
7
  • \$\begingroup\$ Your columns don't make much sense. You have multiple area and multiple azimuth. What is the actual semantic difference from column to column for those areas and azimuths? \$\endgroup\$
    – Reinderien
    Commented Mar 25, 2022 at 15:10
  • \$\begingroup\$ A Timestamp that was dropped from a MultiIndex \$\endgroup\$ Commented Mar 25, 2022 at 15:19
  • \$\begingroup\$ Each param[N] column is associated to a Timestamp. \$\endgroup\$ Commented Mar 25, 2022 at 15:25
  • 1
    \$\begingroup\$ OK. I'm going to suggest that you show those original data with full indices, timestamp included, prior to them receiving your repeat()ed columns. \$\endgroup\$
    – Reinderien
    Commented Mar 25, 2022 at 15:31
  • \$\begingroup\$ The columns have been updated with a MultiIndex. That being said Is there a preferred/more common axis to use timeseries data? \$\endgroup\$ Commented Mar 25, 2022 at 16:46

1 Answer 1

1
\$\begingroup\$

I'm going to assume that you have real timestamps instead of the synthetic 1960 values.

I think that all of the following should go away:

  • np.block()
  • normalize()
  • for param, frame in
  • yield and any kind of generation
  • match and any kind of column-conditional branch
  • mean_stack()
  • squeeze()
  • droplevel()

I don't think that the current approach is as vectorised as it could be.

Also, the data shape coming in is insane. The index is not a "parameter" but is in fact (by your description) a "polygon_id". The second index level should be validTime. That leaves three columns - area, MPS and azimuth - all dependent variables.

Suggested

Not strictly identical because it has a different output shape, but numerically equivalent:

import functools
import pandas as pd
import numpy as np


DATA = [
    [0.045750000000000984, 0.04939999999999974, 0.05244999999999857, 0.054399999999998436, 0.056099999999999456, 0.05375000000000011, 0.047700000000000985, 0.053700000000000594, 0.053249999999999686, 16.81209831859011, 27.983008063348617, 20.301826794339316, 13.421484359287598,
        17.10557180223468, 14.686987606192902, 41.18082589372154, 13.065486605460656, 4.608155337880411, 0.7571951341718915, 0.3457777398771206, 0.46056961460740087, 0.2819901554949584, 0.7205708724903156, -0.00040991561554303056, 0.3860776453927442, 0.09629955035400072, 0.2468722764784534],
    [0.03290000000000123, 0.033550000000000045, 0.03949999999999949, 0.035550000000001615, 0.03725000000000067, 0.03515000000000022, 0.03714999999999999, 0.041350000000001316, 0.03869999999999882, 14.902912040542054, 15.452401255651447, 12.5484544315597, 10.83783196920491,
        13.246734412134916, 12.636336066844674, 19.61093585673077, 8.898327945757917, 26.237157593452878, 0.36589566404637536, 0.5140598344685876, 0.06092353057665044, 0.5433592643289422, 0.3300217799943245, 0.1295695630661193, 0.5813772760582753, 0.97888742111902, 0.2144469684077374],
    [0.15945000000000153, 0.1419000000000014, 0.1346, 0.14314999999999908, 0.14485000000000084, 0.1293000000000014, 0.14264999999999992, 0.1434499999999994, 0.08944999999999861, 39.54982055905117, 67.05598966736439, 14.170464300161013, 25.57103344352307, 22.224605807376317,
        28.901343579758898, 19.05130673155152, 25.61289901736464, 122.31441473369556, 0.5250514668415044, 0.3119136953061388, 0.3818093125132205, 0.2669557624118119, -0.4876722558321518, 0.03814067199937899, -0.07502082862542193, 0.15371372713178985, -1.1026290280269917],
    [0.00340000000000014, 0.0066000000000002515, 0.008499999999999603, 0.007999999999999518, 0.0065500000000002015, 0.008100000000000543, 0.0032500000000004834, 0.018350000000000196, 0.020150000000000542, 22.364276562983136, 20.506557739907425, 6.243036570120132, 19.3458717607394,
        10.456833536986126, 19.65879720239808, 14.307243539019602, 9.091403905263792, 18.47420699835527, 0.8485690702159555, 0.40094252533477115, 0.1157765828210181, 0.6944022345849206, 0.46766736064272035, 1.153657712332056, 0.7600461290670872, 0.4170470750740715, 0.49920935232659946],
    [None, None, 0.003549999999999955, 0.005799999999999788, 0.005799999999999788, 0.006900000000000275, 0.007049999999999682, 0.006800000000000173, 0.005200000000000045, None, None, 8.06110268044255, 13.695219329907925, 0.0,
        12.81872890572318, 6.190861120089245, 3.9861242549543703, 4.403761437912741, None, None, -0.8342338501323102, -0.35990880371032785, 0.0, -2.9708131816772876, -0.12378167714772335, -1.4646911328887278, -2.475146432921943]
]
INDEX = ['89234', '89467', '89519', '89680', '90104']

DATE_RANGE = pd.to_datetime(range(9), unit='m', origin=pd.Timestamp('1960-01-01'))


MULTI_COLUMN = pd.MultiIndex.from_tuples((
    (date, param) for param in ('area', 'MPS', 'azimuth')
    for date in DATE_RANGE), names=['validTime', 'parameters'])


def outliers2nan(bad_df: pd.DataFrame, bad_means: pd.DataFrame, thresh: pd.DataFrame) -> pd.DataFrame:
    """
    where a parameter is greater than the minimum mean threshold set the value to NaN
    """

    area_mps_means = bad_means.loc[:, ('area', 'MPS')]                      # Means (area and MPS)
    area_mps = bad_df.loc[:, ('area', 'MPS')]                               # Dataframe (area and MPS)
    lhs = (area_mps - area_mps_means).abs()                                 # Left-hand side for comparison
    rhs_unaligned = area_mps_means * thresh.loc[:, ('area', 'MPS')].values  # Right-hand side: index on polygon level only
    _, rhs = lhs.align(rhs_unaligned)  # Join-repeat the threshold side to the second level
    good_df = bad_df.copy()            # Good dataframe starts out same as bad
    good_df[lhs > rhs] = np.nan        # Anything over the threshold becomes NaN

    deg_diff = np.rad2deg(bad_df.azimuth - bad_means.azimuth)        # Difference in rad, to deg
    lhs = ((deg_diff + 180) % 360 - 180).abs()                       # Left-hand side for comparison
    good_df.loc[lhs > thresh.azimuth.values[0], 'azimuth'] = np.nan  # Anything over the threshold becomes NaN
    return good_df


def _mean_azimuth_reduction(col1: np.ndarray, col2: np.ndarray) -> np.ndarray:
    col1 = np.nan_to_num(col1, nan=col2)

    mean_rads = np.arctan2(
        np.sin(col1) + np.sin(col2),
        np.cos(col1) + np.cos(col2),
    )

    return mean_rads


def get_means(df: pd.DataFrame) -> pd.DataFrame:
    means = df.loc[:, ('MPS', 'area')].groupby(level='polygon_id').mean()

    grouped_azimuth = np.reshape(
        df.azimuth.values,
        (len(df.groupby(level='polygon_id')), -1)
    ).T
    means['azimuth'] = functools.reduce(_mean_azimuth_reduction, grouped_azimuth)

    return means


def start() -> None:
    """construct df and display means"""

    thresh = pd.DataFrame.from_records(
        [{'area': .65, 'MPS': .65, 'azimuth': 60.0}]
    )

    bad_df = pd.DataFrame(
        DATA, index=INDEX, columns=MULTI_COLUMN,
    ).stack(level=0, dropna=False).rename_axis(('polygon_id', 'validTime'))

    bad_means = get_means(bad_df)

    good_df = outliers2nan(bad_df, bad_means, thresh)
    good_means = get_means(good_df)

    print(f"""
BAD
{bad_df}

GOOD
{good_df}

BAD
{bad_means}

GOOD
{good_means}

DIF
{(bad_means-good_means).abs()}
    """)


if __name__ == '__main__':
    start()
\$\endgroup\$
3
  • \$\begingroup\$ Awesome. I have actual timestamps. They are roughly a 1-2min interval. When you refer to a dataframe as being sane or insane that relates to how it is indexed? I’ve tried googling the concept but not coming up with much. Would you be able to point me to a resource so i can look into it further. \$\endgroup\$ Commented Mar 26, 2022 at 16:07
  • 1
    \$\begingroup\$ Yes, it's about how it's indexed, as well as how columns are named and typed. My apologies: data sanity is not a technical term, and is slightly subjective, but: as a guideline, independent variables (IDs, time) should be indexes and dependent variables (area, MPS, azimuth) should be columns. \$\endgroup\$
    – Reinderien
    Commented Mar 26, 2022 at 18:15
  • \$\begingroup\$ That makes a lot of sense. Thanks. \$\endgroup\$ Commented Mar 26, 2022 at 18:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.