N0 calibration (weighting procedure)

[danpower101]

This is part of the review of the openjournals/joss-reviews#6025

The weighting procedure currently implemented has since been revised (Schrön et al., 2017).

It is a refinement of the Köhli 2015 weighting already implemented in this code (it now takes into account the influence of pressure and humidity, as well as allowing the weighting of samples <0.5m from the sensor.

I feel this should be added as an option to allow users to match the weighting strategy of the networks currently implementing it (e.g., https://cosmoz.csiro.au/about). It involves some different weights and a convergence approach (described in section 2.3 from the paper).

Schrön, M., Köhli, M., Scheiffele, L., Iwema, J., Bogena, H. R., Lv, L., Martini, E., Baroni, G., Rosolem, R., Weimar, J., Mai, J., Cuntz, M., Rebmann, C., Oswald, S. E., Dietrich, P., Schmidt, U., and Zacharias, S.: Improving calibration and validation of cosmic-ray neutron sensors in the light of spatial sensitivity, Hydrol. Earth Syst. Sci., 21, 5009–5030, https://doi.org/10.5194/hess-21-5009-2017, 2017.

[joaquinperaza]

@danpower101 Please feel free to check if the last changes in the software addressed your suggestions.

[danpower101]

I couldn't find changes in the main branch so I presume this is it:
https://github.com/soilwater/crnpy/blob/0.6.0/src/crnpy/crnpy.py

This part needs a few more adjustments.

Line 943:945

This step shouldn't be optional for the Schrön 2017 method.

The impact of pressure on the weighting steps is introduced through the rscaled function, so allowing this to be skipped removes this impact. All sites will record pressure as standard, so it should always be available.

Hveg you might want to default to 0 for when data is not available, as it is less likely to be recorded during calibration campaigns (at least with older sites this is the case). It would be better to be included if it is available though.

Weighting

• When bulk density used in the following calculations it should be the site average, whereas I think currently the calculations are using the specific bulk density of each sample. This is because neutrons are impacted by the bulk density across the whole area and not just the sample area.

The steps here should be:

1. Estimate field scale SM (start with equal average - y_average)
2. Re-scale radius - using pressure and (if available) Hveg
3. Calculate D86 of each vertical profile using the average field SM estimate (y_average) and the rescaled r parameter.
4. Weight each profile using the vertical weighting functions to get an average SM value for each profile
5. Horizontally weight the profiles, using the horizontal weighting scheme, the theta (or y) should be the average SM value (y_average). The influence of theta/y on neutrons is thought of as being from the whole field which is why we have to estimate it here and then keep running the weighting until it converges on a value.

e.g., this:

Wr[r0_idx] = WrX(r[r0_idx], x[r0_idx], y[r0_idx])

becomes this:

Wr[r0_idx] = WrX(r[r0_idx], x_average, y_average)

...absolute humidity (x) and theta (y) are from average values, and the radius includes the sensitivity to pressure/Hveg as calculated in the r-scale step.

6. Use the horizontal weights to average the SM profiles, and get a revised estimate (y_average_new)
7. Compare the initial estimate (y_average) with the new estimate (y_average_new), to see if they are converging (Schrön et al., 2017 suggest 1% difference).
8. Repeat as needed, until the values converge, using y_average_new as your new estimate of SM (y_average).

[danpower101]

I was able to work with the .ipynb notebook for the calibration (which is working).

I noticed however that the profile is being created using the distance from the sensor, which it shouldn't. Doing it this way will mean, for example, that a profile taken at North 25m will be combined with a profile taken at South 25m and treated as a single profile. Better to use individual profile assignments to make sure they're treated separately.

[joaquinperaza]

@danpower101, I fixed the method to require a unique profile indetifier, in the example I created a new column ID to compensate for the lack of that information in the sample dataset, also the unit tests were updated to account for this additional input.

[danpower101]

Works for me, happy to close this issue.