Skip to content

Commit

Permalink
Ticket 51: Update HAC quality indicator field
Browse files Browse the repository at this point in the history
  • Loading branch information
andershenja committed Oct 28, 2024
1 parent f8131d7 commit 0ed5d52
Show file tree
Hide file tree
Showing 3 changed files with 155 additions and 28 deletions.
1 change: 1 addition & 0 deletions Lib/odc_hac.py
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,7 @@ def hacFilter(self, scan, quant="DBZH", enough=100):
self.thresh = ARGS["default"].thresh
## Got site-specific threshold?


qind = _ravefield.new()
qind.setData(zeros(hac_data.shape, uint8))
qind.addAttribute("how/task", "eu.opera.odc.hac")
Expand Down
103 changes: 79 additions & 24 deletions librave/toolbox/odc_hac.c
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,11 @@ along with HLHDF. If not, see <http://www.gnu.org/licenses/>.
*/

#include "odc_hac.h"

#include "polarscan.h"
#include "rave_attribute.h"
#include "rave_field.h"
#include "rave_types.h"
#include <stdio.h>

int hacFilter(PolarScan_t* scan, RaveField_t* hac, char* quant) {
PolarScanParam_t* param = NULL;
Expand All @@ -39,34 +43,85 @@ int hacFilter(PolarScan_t* scan, RaveField_t* hac, char* quant) {
nrays = PolarScan_getNrays(scan);

if (PolarScan_hasParameter(scan, quant)) {
param = PolarScan_getParameter(scan, quant);
qind = PolarScan_getQualityFieldByHowTask(scan, "eu.opera.odc.hac");
nodata = PolarScanParam_getNodata(param);
param = PolarScan_getParameter(scan, quant);
qind = PolarScan_getQualityFieldByHowTask(scan, "eu.opera.odc.hac");
nodata = PolarScanParam_getNodata(param);

if (qind == NULL) {
qind = RAVE_OBJECT_NEW(&RaveField_TYPE);
if (qind != NULL) { // (scan.nrays, scan.nbins)
if (!RaveField_createData(qind, nrays, nbins, RaveDataType_UCHAR)) {
RAVE_ERROR0("Failed to create data field");
goto done;
}
attr = RaveAttributeHelp_createString("how/task", "eu.opera.odc.hac");
if (attr == NULL) {
RAVE_ERROR0("Failed to add how/task to quality field");
goto done;
}
if (!RaveField_addAttribute(qind, attr)) {
RAVE_ERROR0("Failed to add attribute to qfield");
goto done;
}
RAVE_OBJECT_RELEASE(attr);
if (!PolarScan_addQualityField(scan, qind)) {
RAVE_ERROR0("Failed to add hac quality field to scan");
goto done;
}
} else {
RAVE_ERROR0("Failed to create rave data field");
goto done;
}
}

attr = RaveField_getAttribute(qind, "how/task_args");
RaveAttribute_getDouble(attr, &thresh);
RAVE_OBJECT_RELEASE(attr);
attr = RaveField_getAttribute(hac, "how/count");
RaveAttribute_getLong(attr, &N);

for (ir=0; ir<nrays; ir++) {
for (ib=0; ib<nbins; ib++) {
rvt = PolarScanParam_getValue(param, ib, ir, &val);
attr = RaveField_getAttribute(qind, "how/task_args");
if (attr == NULL) {
thresh = 60.0;
attr = RaveAttributeHelp_createDouble("how/task_args", thresh);
if (attr == NULL || !RaveField_addAttribute(qind, attr)) {
RAVE_ERROR0("Failed to add how/task_args to quality field");
goto done;
}
} else {
RaveAttribute_getDouble(attr, &thresh);
}
RAVE_OBJECT_RELEASE(attr);

attr = RaveAttributeHelp_createDouble("what/offset", 0.0);
if (attr == NULL || !RaveField_addAttribute(qind, attr)) {
RAVE_ERROR0("Failed to add what/offset to quality field");
goto done;
}
RAVE_OBJECT_RELEASE(attr);

if (rvt==RaveValueType_DATA) {
RaveField_getValue(hac, ib, ir, &ni);
Pi = 100 * (ni/(double)N);
attr = RaveAttributeHelp_createDouble("what/gain", 1.0/255.0);
if (attr == NULL || !RaveField_addAttribute(qind, attr)) {
RAVE_ERROR0("Failed to add what/gain to quality field");
goto done;
}
RAVE_OBJECT_RELEASE(attr);

if (Pi > thresh) {
PolarScanParam_setValue(param, ib, ir, nodata);
RaveField_setValue(qind, ib, ir, val);
}
}
}
}
attr = RaveField_getAttribute(hac, "how/count");
RaveAttribute_getLong(attr, &N);

for (ir=0; ir<nrays; ir++) {
for (ib=0; ib<nbins; ib++) {
rvt = PolarScanParam_getValue(param, ib, ir, &val);
RaveField_setValue(qind, ib, ir, 255);
if (rvt==RaveValueType_DATA) {
RaveField_getValue(hac, ib, ir, &ni);
Pi = 100 * (ni/(double)N);
if (Pi > thresh) {
PolarScanParam_setValue(param, ib, ir, nodata);
RaveField_setValue(qind, ib, ir, 0);
}
}
}
}

retval = 1;
retval = 1;
}
done:
RAVE_OBJECT_RELEASE(param);
RAVE_OBJECT_RELEASE(qind);
RAVE_OBJECT_RELEASE(attr);
Expand Down
79 changes: 75 additions & 4 deletions test/pytest/odc_hac_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,8 @@
import math
import string
import _odc_hac, odc_hac, rave_zdiff_quality_plugin
import _raveio, _ravefield
import _polarscanparam,_polarvolume
import _raveio, _ravefield, _rave
import _polarscan,_polarscanparam,_polarvolume
import numpy

class odc_hac_test(unittest.TestCase):
Expand Down Expand Up @@ -139,7 +139,7 @@ def test_rave_zdiff_quality_plugin_process_pvol(self):
scan1 = self.create_scan(3.0, 4.0, 53.0, 40.0)
scan2 = self.create_scan(4.0, 3.0, 40.0, 53.0)
scan2.elangle = 10.0 * math.pi / 180.0

pvol = _polarvolume.new()
pvol.addScan(scan1)
pvol.addScan(scan2)
Expand All @@ -163,7 +163,78 @@ def test_rave_zdiff_quality_plugin_getQualityFields(self):
fields = qp.getQualityFields()
self.assertEqual(1, len(fields))
self.assertEqual("eu.opera.odc.zdiff", fields[0])


def test_odcFilter(self):
scan = self.create_odcFilter_scan()

hac = _ravefield.new()
hac.addAttribute("how/count", 4)
hac.setData(numpy.asarray([[1,2],[3,4]], numpy.int32))

_odc_hac.hacFilter(scan, hac, "DBZH")

qf = scan.getQualityFieldByHowTask("eu.opera.odc.hac")
self.assertTrue((numpy.asarray([[255,255],[0,0]], numpy.uint8) == scan.getQualityFieldByHowTask("eu.opera.odc.hac").getData()).all())
self.assertAlmostEqual(0.0, qf.getAttribute("what/offset"), 4)
self.assertAlmostEqual(1.0/255.0, qf.getAttribute("what/gain"), 4)
self.assertAlmostEqual(60.0, qf.getAttribute("how/task_args"), 4)

def test_odcFilter_addAttributes(self):
scan = self.create_odcFilter_scan(None)

hac = _ravefield.new()
hac.addAttribute("how/count", 4)
hac.setData(numpy.asarray([[1,2],[3,4]], numpy.int32))

_odc_hac.hacFilter(scan, hac, "DBZH")

qf = scan.getQualityFieldByHowTask("eu.opera.odc.hac")
self.assertTrue((numpy.asarray([[255,255],[0,0]], numpy.uint8) == scan.getQualityFieldByHowTask("eu.opera.odc.hac").getData()).all())
self.assertAlmostEqual(0.0, qf.getAttribute("what/offset"), 4)
self.assertAlmostEqual(1.0/255.0, qf.getAttribute("what/gain"), 4)
self.assertAlmostEqual(60.0, qf.getAttribute("how/task_args"), 4)

def test_odcFilter_addFieldAndAttributes(self):
scan = self.create_odcFilter_scan(None, False)

hac = _ravefield.new()
hac.addAttribute("how/count", 4)
hac.setData(numpy.asarray([[1,2],[3,4]], numpy.int32))

_odc_hac.hacFilter(scan, hac, "DBZH")

qf = scan.getQualityFieldByHowTask("eu.opera.odc.hac")
self.assertTrue((numpy.asarray([[255,255],[0,0]], numpy.uint8) == scan.getQualityFieldByHowTask("eu.opera.odc.hac").getData()).all())
self.assertAlmostEqual(0.0, qf.getAttribute("what/offset"), 4)
self.assertAlmostEqual(1.0/255.0, qf.getAttribute("what/gain"), 4)
self.assertAlmostEqual(60.0, qf.getAttribute("how/task_args"), 4)


def create_odcFilter_scan(self, thresh=60.0, addQualityField=True):
scan = _polarscan.new()

param = _polarscanparam.new()
param.quantity="DBZH"
param.setData(numpy.asarray([[1,2],[3,4]], numpy.uint8))
scan.addParameter(param)

param = _polarscanparam.new()
param.quantity="TH"
param.setData(numpy.asarray([[2,3],[4,5]], numpy.uint8))
scan.addParameter(param)

if addQualityField:
qind = _ravefield.new()
qind.setData(numpy.asarray([[2,2],[3,3]], numpy.uint8))
qind.addAttribute("how/task", "eu.opera.odc.hac")
if thresh is not None :
qind.addAttribute("how/task_args", thresh)

scan.addQualityField(qind)

return scan


def create_scan(self, dbzh0_0=3.0, dbzh0_1=4.0, th0_0=53.0, th0_1=40.0):
scan = _raveio.open(self.SCAN_FIXTURE).object
param_dbzh = scan.getParameter("DBZH")
Expand Down

0 comments on commit 0ed5d52

Please sign in to comment.