SuperDARN / rst

Radar Software Toolkit (RST)
https://superdarn.github.io/rst/
GNU General Public License v3.0
22 stars 17 forks source link

threshold parameter in raw files #142

Closed egthomas closed 6 years ago

egthomas commented 6 years ago

This issue originates from the dattorawacf buffer overflow previously addressed in #116.

The following circular series of dat -> rawacf -> dat conversions does not reproduce the original dat file:

dattorawacf 2005123120t.dat > 2005123120t.rawacf
rawacftodat 2005123120t.rawacf > 2005123120t.dat2
cmpraw 2005123120t.dat 2005123120t.dat2 > diff.txt

However this circular series of conversions does reproduce the original dat file (note the -t 0 option with dattorawacf:

dattorawacf -t 0 2005123120t.dat > 2005123120t.rawacf
rawacftodat 2005123120t.rawacf > 2005123120t.dat2
cmpraw 2005123120t.dat 2005123120t.dat2 > diff.txt

(setting the -t 0 option with rawacftodat does not appear to have an effect on either example)

The reason for this discrepancy is because RawEncode in rawwrite.c of the raw library compares the lag0 power at each range gate to the floor of half the threshold times the noise calculated from the clear frequency search:

https://github.com/SuperDARN/rst/blob/e6e9a1a0427436555dc84eea1e374becf0f59921/codebase/superdarn/src.lib/tk/raw.1.22/src/rawwrite.c#L60-L101

Similarly, OldRawWrite in rawwrite.c of the oldraw library also checks against the threshold and noise:

https://github.com/SuperDARN/rst/blob/e6e9a1a0427436555dc84eea1e374becf0f59921/codebase/superdarn/src.lib/tk/oldraw.1.16/src/rawwrite.c#L246-L274

I've asked Dieter and Mike Ruohoniemi about this and both seem to think the threshold parameter was used to reduce file sizes. It appears to be set at each radar site as a command line option of rawacfwrite; if no value is provided, the default is zero (although the code Dieter sent me suggests it used to default to 3 for older qnx4 versions of rawacfwrite).

Returning to my example above, the threshold value stored in the 2005123120t.dat file is already set to 3. It therefore doesn't make sense to me why passing the dat records through the rawwrite logic would remove additional records if some thresholding was already applied at the radar site as implied by the value stored in the file header.

Another troubling part is that this conversion discrepancy also appears in fit files. For example, the following set of commands illustrates how a dattorawacf conversion results in different output from make_fit:

make_fit -old 2005123120t.dat 2005123120t.fit
dattorawacf 2005123120t.dat > 2005123120t.rawacf
rawacftodat 2005123120t.rawacf > 2005123120t.dat2
make_fit -old 2005123120t.dat2 2005123120t.fit2
cmpfit 2005123120t.fit 2005123120t.fit2 > diff.txt

However, repeating these commands with the -t 0 option does not show any difference in the resulting fit files:

make_fit -old 2005123120t.dat 2005123120t.fit
dattorawacf -t 0 2005123120t.dat > 2005123120t.rawacf
rawacftodat 2005123120t.rawacf > 2005123120t.dat2
make_fit -old 2005123120t.dat2 2005123120t.fit2
cmpfit 2005123120t.fit 2005123120t.fit2 > diff.txt

To conclude, I think it's a bad thing that by default our file format conversion routines do not reproduce the original dat files. This can result in differences in the higher level fit files produced from these dat files and therefore impact data reproducibility. I can think of a couple of ways to address this issue:

1) Set the default threshold value in the dattorawacf, rawacftodat, trim_raw, and cat_raw binaries to zero. This will overwrite the threshold values stored in the files and prevent any additional records from being removed. 2) Remove the comparisons between lag0 power and the threshold/noise from the various rawwrite functions entirely. This is probably the more extreme option and may well be inconsistent with the intended behavior at the radar site or elsewhere.

This problem seems to straddle the responsibilities of the Data Analysis, Data Distribution, and Operating Software working groups, so please share this with anyone else who may have some insight.

pasha-ponomarenko commented 6 years ago

Thanks a lot , Evan! You clarified some issues that were teasing me for some time. In #116 I wrote that

some .dat files ... would have all pwr0 values but some of the respective low-power ACFs below a certain threshold would not be recorded.

Besides the reproducibility problem, there are two more side effects that are detrimental from point of view of data analysis: (1) It makes little sense to use the clear frequency search, that is performed during the duraition of a single pulse sequence, for filtering data averaged over 3.5-7 s. The statistical fluctuation level for ths estimate is 1/sqrt(N_ave) higher than that of the averaged data, so this precedure is grossly inferior to the one based on the 10 weakest lag 0 power values (although the latter also required some corrections). (2) The noise level estimated by FITACF and used for determining SNR is based on lag 0 power from raw.acfd. Due to pre-filtering, these data woud have some zeroes that are ignored while calculating the average of the of the ten weakest ACFs so that this noise estimate is generally somewhat higher (though not by too much!) than the actual noise level.

Furthermore, the use of the fixed threshold value of 3 would have very different effect on data obtained with different A/D converters. I observed before that the .raw data supposedly recorded with higher resolution have much fewer zero records than those recorded in .dat files.

So, I wholeheartedly agree that any preliminary filtering should be abandoned during both recording the data and performing old/new type conversion.

ksterne commented 6 years ago

I don't think I can really add much here. I think this is probably before me and I haven't done a lot of work with dat or fit files.

The only experience I think I can lend here is that in the dattorawacf program, I think there were cases were the operations parameters were written, but the acfd array was not. I ended up making one of the older pull requests on the VTRST 3.5 repo where a fitacf library was assuming for each "record" there was an acfd array. So, this could be apart of the problem I would imagine.

egthomas commented 6 years ago

There are multiple references to this thresholding procedure throughout the RADOPS 2000 documentation (attached below). For example, on page 35 describing the raw_write task:

By specifying the “-t” option a threshold value can be applied to the lag-zero power. Data below this threshold is not stored in the file providing a further reduction in the size of the output file. The threshold function is defined as:

threshold*NOISE/2.

Data with lag-zero power less than the result of the above sum is rejected

If this were something that was only ever applied at the radar site then it's something we'll have to live with, but the fact that additional data are removed when a dat file is converted or rewritten should not be the default behavior.

reference.pdf

egthomas commented 6 years ago

Merged in #146.