euroargodev / User-Acceptance-Test-Python-version-of-the-OWC-tool

BODC Argo team with consultation with other Argo partners took initiative to convert the OWC software from Matlab to Python. This initiative has been undertaken based on the output from the international survey about the methods and tools used in DMQC core data. The following repository is dedicated for User Acceptance Testing.
1 stars 0 forks source link

Out of bounds for ptmp in find_10thetas #19

Closed kamwal closed 4 years ago

kamwal commented 4 years ago

In running the Python code I found the following error in find_10theta:

L:\users\argo\kamwal\OWC_Python_Test_2\argodmqc_owc\ow_calibration\find_10thetas\find_10thetas.py:175: RuntimeWarning: invalid value encountered in greater if np.nanmax(ptmp[:, depth]) > theta_levels[level] > np.nanmin(ptmp[:, depth]): L:\users\argo\kamwal\OWC_Python_Test_2\argodmqc_owc\ow_calibration\find_10thetas\find_10thetas.py:203: RuntimeWarning: invalid value encountered in greater pos_diff = np.argwhere(ptmp_diff > 0) L:\users\argo\kamwal\OWC_Python_Test_2\argodmqc_owc\ow_calibration\find_10thetas\find_10thetas.py:213: RuntimeWarning: invalid value encountered in less neg_diff = np.argwhere(ptmp_diff < 0) Traceback (most recent call last): File "owc_calibration.py", line 81, in <module> calc_piecewisefit("/", FLOAT_NAME, USER_CONFIG) File "L:\users\argo\kamwal\OWC_Python_Test_2\argodmqc_owc\ow_calibration\calc_piecewisefit\calc_piecewisefit.py", line 230, in calc_piecewisefit theta, p, index, var_s_th, th = find_10thetas(copy.deepcopy(unique_sal), File "L:\users\argo\kamwal\OWC_Python_Test_2\argodmqc_owc\ow_calibration\find_10thetas\find_10thetas.py", line 200, in find_10thetas ptmp_diff = ptmp[theta_index, depth] - ptmp[interval, depth] IndexError: index 57 is out of bounds for axis 0 with size 57

WMO: 1901847 Reference data: CTD_for_DMQC_2019V01 WMO boxes: wmo_boxes_ctd.mat Canfiguration: `# Objective Mapping Parameters #

    # max number of historical casts used in objective mapping
    'CONFIG_MAX_CASTS': 300,

    # 1=use PV constraint, 0=don't use PV constraint, in objective mapping
    'MAP_USE_PV': 0,

    # 1=use SAF separation criteria, 0=don't use SAF separation criteria, in objective mapping
    'MAP_USE_SAF': 0,

    # spatial decorrelation scales, in degrees
    'MAPSCALE_LONGITUDE_LARGE': 8,
    'MAPSCALE_LONGITUDE_SMALL': 4,
    'MAPSCALE_LATITUDE_LARGE': 4,
    'MAPSCALE_LATITUDE_SMALL': 2,

    # cross-isobath scales, dimensionless, see BS(2005)
    'MAPSCALE_PHI_LARGE': 0.1,
    'MAPSCALE_PHI_SMALL': 0.02,

    # temporal decorrelation scale, in years
    'MAPSCALE_AGE_LARGE': 20,
    'MAPSCALE_AGE_SMALL': 5,

    # exclude the top xxx dbar of the water column
    'MAP_P_EXCLUDE': 100,

    # only use historical data that are within +/- yyy dbar from float data
    'MAP_P_DELTA': 200
}`
kamwal commented 4 years ago

In the cal_piecewisefit.py function we found a bug where "unique_ptmp" has been called twice.

theta,p, index, var_s_th, th = find_10thetas(copy.deepcopy(unique_sal),
                                                      copy.deepcopy(unique_ptmp),
                                                      copy.deepcopy(unique_pres),
                                                      copy.deepcopy(unique_ptmp),
                                                      use_theta_lt[0, 0], use_theta_gt[0, 0],
                                                      use_pres_lt[0, 0], use_pres_gt[0, 0],
                                                      use_percent_gt[0, 0])

The code is not including ptmp data from mapped file. This needs to be changed to:

theta, p, index, var_s_th, th = find_10thetas(copy.deepcopy(unique_sal),
                                                      copy.deepcopy(unique_ptmp),
                                                      copy.deepcopy(unique_pres),
                                                      copy.deepcopy(unique_mapped_ptmp),
                                                      use_theta_lt[0, 0], use_theta_gt[0, 0],
                                                      use_pres_lt[0, 0], use_pres_gt[0, 0],
                                                      use_percent_gt[0, 0])
edsmall-bodc commented 4 years ago

During the algorithm to find 10 theta levels we need to find all the indices in each column that does not contain a nan. We then take the potential temperature at this point, and find out how much the potential temperature in the previous and next water column differs from this point, like so:

ptmp_diff =potential temperature here - [last water column at same depth, the current water column, the next water column at the same depth]

However, there is a problem here. What if we find a good ptmp for the very last water column? Then there is no next column

The way we get around this is that, if we are at the very last column, just use the last value again when doing the diff. However, I have an "off by one" error here:

interval = np.arange(np.max([theta_index - 1, 0]), np.min([theta_index + 1, profile_no]) + 1, dtype=int)

This is error is because the matlab indexing starts at 1, but python starts at zero. profile_no is the number of profiles we have (which in this case is 57), but that means the indexing goes from 0 to 56! Therefore, what we want is:

interval = np.arange(np.max([theta_index - 1, 0]), np.min([theta_index + 1, profile_no - 1]) + 1, dtype=int)

edsmall-bodc commented 4 years ago

@kamwal A fix for this problem can be found HERE

Checkout the branch and let me know if the problem is resolved. The branch does not contain your config file or your mapped output and float source.

edsmall-bodc commented 4 years ago

Outputs seem to match to me. Let me know if the branch resolves this.

image

image

edsmall-bodc commented 4 years ago

Issue has been resolved and pulled in here. Closing.