Ticket 817: Some parameters added (EMAX, FF_MAX and NZ). The command-line c-code...
authorulf.nordh <a002011@c21745.ad.smhi.se>
Mon, 3 Feb 2020 10:08:56 +0000 (11:08 +0100)
committerulf.nordh <a002011@c21745.ad.smhi.se>
Mon, 3 Feb 2020 10:08:56 +0000 (11:08 +0100)
13 files changed:
bin/Makefile
bin/fix_shebang.sh [new file with mode: 0755]
bin/wrwp_config.xml [new file with mode: 0644]
bin/wrwp_main [new file with mode: 0644]
bin/wrwp_main.c [deleted file]
lib/wrwp.c
lib/wrwp.h
pywrwp/VPodimVersionConverter.py [new file with mode: 0644]
pywrwp/baltrad_wrwp_pgf_plugin.py
pywrwp/pywrwp.c
test/pytest/WrwpTest.py
test/pytest/fixtures/seosd_qcvol_zdrvol_different_task.h5 [new file with mode: 0644]
test/pytest/fixtures/wrwp_config.xml [new file with mode: 0644]

index 4bfe2ab..cb1e965 100644 (file)
@@ -1,92 +1,40 @@
 ###########################################################################
-# Copyright (C) 2013 Swedish Meteorological and Hydrological Institute, SMHI,
+# Copyright (C) 2012 Swedish Meteorological and Hydrological Institute, SMHI,
 #
-# This file is part of baltrad-wrwp.
+# This file is part of beamb.
 #
-# baltrad-wrwp is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
+# beamb is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 # 
-# baltrad-wrwp is distributed in the hope that it will be useful,
+# beamb is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
+# GNU Lesser General Public License for more details.
 # 
-# You should have received a copy of the GNU General Public License
-# along with baltrad-wrwp.  If not, see <http://www.gnu.org/licenses/>.
+# You should have received a copy of the GNU Lesser General Public License
+# along with RAVE.  If not, see <http://www.gnu.org/licenses/>.
 # ------------------------------------------------------------------------
 # 
-# baltrad-wrwp make file
+# wrwp_main make file
 # @file
 # @author Anders Henja (Swedish Meteorological and Hydrological Institute, SMHI)
-# @date 2013-09-16
+# @date 2012-08-30
 ###########################################################################
 -include ../def.mk
 
-# --------------------------------------------------------------------
-# Fixed definitions
-CFLAGS=-I../lib $(LAPACKE_INCLUDE_DIR) $(CBLAS_INCLUDE_DIR) $(RAVE_MODULE_CFLAGS)
-LDFLAGS=-L../lib $(BLAS_LIB_DIR) $(CBLAS_LIB_DIR) $(LAPACK_LIB_DIR) $(LAPACKE_LIB_DIR) $(RAVE_MODULE_LDFLAGS) $(FORTRAN_CLINK_LDFLAG)
-LIBS= -lwrwp $(RAVE_MODULE_LIBRARIES) -llapacke -llapack -lcblas -lblas $(FORTRAN_CLINK_LIBS) -lm
-
-SOURCES= wrwp_main.c
-OBJECTS= $(SOURCES:.c=.o)
-TARGET=        wrwp
-
-MAKEDEPEND=gcc -MM $(CFLAGS) -o $(DF).d $<
-DEPDIR=.dep
-DF=$(DEPDIR)/$(*F)
-
-# --------------------------------------------------------------------
-# Rules
-
-# Contains dependency generation as well, so if you are not using
-# gcc, comment out everything until the $(CC) statement.
-%.o : %.c
-       @$(MAKEDEPEND); \
-       cp $(DF).d $(DF).P; \
-       sed -e 's/#.*//' -e 's/^[^:]*: *//' -e 's/ *\\$$//' \
-               -e '/^$$/ d' -e 's/$$/ :/' < $(DF).d >> $(DF).P; \
-       \rm -f $(DF).d
-       $(CC) -c $(CFLAGS) $<
-
-# Ensures that the .dep directory exists
-.PHONY=$(DEPDIR)
-$(DEPDIR):
-       +@[ -d $@ ] || mkdir -p $@
-
-.PHONY=all
-all:           $(TARGET)
-
-$(TARGET): $(DEPDIR) $(OBJECTS)
-       $(CC) $(LDFLAGS) -o $@ $(OBJECTS) $(LIBS)
-
 .PHONY=install
 install:
-       @"$(HLHDF_INSTALL_BIN)" -f -o -C $(TARGET) "${DESTDIR}$(prefix)/bin/$(TARGET)"
+       @mkdir -p "${DESTDIR}${prefix}/bin/"
+       @./fix_shebang.sh ${PYTHON_BIN} wrwp_main "${DESTDIR}${prefix}/bin/"
+       @mkdir -p "${DESTDIR}${prefix}/config/"
+       @cp -v -f *.xml "${DESTDIR}${prefix}/config/"
 
 .PHONY=clean
 clean:
-               @\rm -f *.o core *~ $(TARGET)
-               @\rm -fr $(DEPDIR)
+       @\rm -f *.o core *~
+       @\rm -fr $(DEPDIR)
 
 .PHONY=distclean                
 distclean:     clean
-               @\rm -f *.so
-
-# NOTE! This ensures that the dependencies are setup at the right time so this should not be moved
--include $(SOURCES:%.c=$(DEPDIR)/%.P)
-
-#.PHONY=install
-#install:
-#      @mkdir -p ${prefix}/bin/
-#      @cp -v -f beamb ${prefix}/bin/
-
-#.PHONY=clean
-#clean:
-#      @\rm -f *.o core *~
-#      @\rm -fr $(DEPDIR)
-
-#.PHONY=distclean               
-#distclean:    clean
diff --git a/bin/fix_shebang.sh b/bin/fix_shebang.sh
new file mode 100755 (executable)
index 0000000..7e72ce1
--- /dev/null
@@ -0,0 +1,35 @@
+#!/bin/sh
+###########################################################################
+# Copyright (C) 2012 Swedish Meteorological and Hydrological Institute, SMHI,
+#
+# This file is part of RAVE.
+#
+# RAVE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# 
+# RAVE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+# 
+# You should have received a copy of the GNU Lesser General Public License
+# along with RAVE.  If not, see <http://www.gnu.org/licenses/>.
+# ------------------------------------------------------------------------
+# 
+# Fixes the shebang line in the python scripts
+# @file
+# @author Anders Henja (Swedish Meteorological and Hydrological Institute, SMHI)
+# @date 2018-04-01
+############################################################################
+
+if [ $# -ne 3 ]; then
+  echo "Must specify <python binary> <file to be modifed> <full path to installation>"
+  exit 127
+fi
+
+cat $2 | sed -e '1s/.*python.*/#!\/usr\/bin\/env '$1'/' > "$3/$2"
+chmod 755 "$3/$2"
+echo "Copied $2 => $3/$2"
+
diff --git a/bin/wrwp_config.xml b/bin/wrwp_config.xml
new file mode 100644 (file)
index 0000000..aa69cd7
--- /dev/null
@@ -0,0 +1,63 @@
+<?xml version="1.0"?>
+<wrwp-params>
+    <param name="QUANTITIES">
+        <description>"Currently supported fields"</description>
+        <value>NV,HGHT,UWND,VWND,ff,ff_dev,dd,DBZH,DBZH_dev,NZ</value>
+    </param>
+    <param name="DMIN">
+        <description>Min distance for deriving a profile [m]</description>
+        <value>5000</value>
+    </param>
+    <param name="DMAX">
+        <description>Max distance for deriving a profile [m]</description>
+        <value>25000</value>
+    </param>
+    <param name="NMIN_WND">
+        <description>Minimum sample size for wind</description>
+        <value>40</value>
+    </param>
+    <param name="NMIN_REF">
+        <description>Minimum sample size for wind</description>
+        <value>40</value>
+    </param>
+    <param name="EMIN">
+        <description>Minimum elevation angle [deg]</description>
+        <value>0.5</value>
+    </param>
+    <param name="EMAX">
+        <description>Minimum elevation angle [deg]</description>
+        <value>45.0</value>
+    </param>
+    <param name="VMIN">
+        <description>Radial velocity threshold [m/s]</description>
+        <value>2.0</value>
+    </param>
+    <param name="FF_MAX">
+        <description>Maximum allowed layer velocity [m/s][</description>
+        <value>60.0</value>
+    </param>
+    <param name="DZ">
+        <description>Height interval for deriving a profile [m]</description>
+        <value>200</value>
+    </param>
+    <param name="HMAX">
+        <description>Maximum height of the profile [m]</description>
+        <value>12000</value>
+    </param>
+    <param name="NODATA_VP">
+        <description>Nodata value for profile</description>
+        <value>-9999</value>
+    </param>
+    <param name="UNDETECT_VP">
+        <description>Undetect value for profile</description>
+        <value>-9999</value>
+    </param>
+    <param name="GAIN_VP">
+        <description>Gain value for profile</description>
+        <value>1.0</value>
+    </param>
+    <param name="OFFSET_VP">
+        <description>Offset value for profile</description>
+        <value>0.0</value>
+    </param>
+</wrwp-params>
diff --git a/bin/wrwp_main b/bin/wrwp_main
new file mode 100644 (file)
index 0000000..2d37207
--- /dev/null
@@ -0,0 +1,381 @@
+#!/usr/bin/env python3
+'''
+Copyright (C) 2012- Swedish Meteorological and Hydrological Institute (SMHI)
+
+This file is part of RAVE.
+
+RAVE is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+RAVE is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+See the GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with RAVE.  If not, see <http://www.gnu.org/licenses/>.'''
+
+#Makes weather radar vertical profiles directly from polar volumes. Input volumes
+#can be either ver2.1 or ver2.2. In the case with input volumes of ver2.1 the
+#produced vertical profiles will be ver2.1. In the case of input volumes of
+#ver2.2 the output vertical profiles can be either v2.1 or v2.2.
+#
+#Adjustable parameters is placed in the file wrwp_config.xml which is placed under .../baltrad-wrwp/bin/ as default. 
+#The necessary parameters in the wrwp_config.xml can of course also be adjusted from the command-line when executing
+#the code.
+## @file
+## @author Ulf E. Nordh, SMHI
+## @date 2019-11-05
+
+import _wrwp
+import _rave
+import _raveio
+import _polarvolume
+import string
+import rave_tempfile
+import odim_source
+import math
+import xml.etree.cElementTree as ET    
+from rave_defines import CENTER_ID, GAIN, OFFSET, H5RAD_VERSIONS
+import os
+import time
+import _pyhl
+import VPodimVersionConverter
+import sys
+import fnmatch
+from optparse import OptionParser
+import logging, logging.handlers
+
+# The default profile output directory if nothing else is selected
+DEFAULTOUT = os.path.dirname(os.path.realpath(__file__)) + '/'
+
+def find(pattern, path):
+  # Locates a file (pattern) in a directory tree (path)
+  result = []
+  for root, dirs, files in os.walk(path):
+    for name in files:
+      if fnmatch.fnmatch(name, pattern):
+        result.append(os.path.join(root, name))
+  return result
+
+def strToNumber(strval):
+  # Converts a string into a number, either int or float
+  # strval: the string to translate
+  # return:the translated value
+  # throws ValueError if value not could be translated
+  if type(strval) is not str: # Avoid doing anything if it is not a string as input
+    return strval
+  else:
+    try:
+      return int(strval)
+    except ValueError:
+      return float(strval)
+
+def write_file(out_file, dst,exestart):
+
+  # Writing of the generated converted file + check if the writing succeded or failed    
+  dst.write(out_file, 6)
+  exetime2 = time.time() - exestart
+  if os.path.isfile(out_file):
+    logger.debug("Succeded in writing ver2.1 vertical profile to disk: " + out_file)
+    logger.debug("Total generation time for ver 2.1 vertical profile: %f s"%exetime2)
+    logger.info("Finished generating ver 2.1 vertical profile, file: %s\n"%out_file)
+  else:
+    logger.info("Failed to write ver2.1 vertical profile to disk: \n" + out_file)
+  return
+
+def change2Odim21(filename, quantities, options, exestart):
+
+  #Initiate the converter
+  converter = VPodimVersionConverter.VPodimVersionConverter(filename, quantities)
+  
+  if converter.isSupported():
+    src, prods = converter.convert(quantities, filename)
+    for dst,out_file in prods:
+      write_file(out_file, dst, exestart)
+  return
+
+def main(options):
+  
+  ## Creates a vertical profile
+
+  files = options.infiles.split(",")
+
+  if options.outfname != None and len(files) > 1:
+    raise AttributeError("Several infiles given, but only one outfile. Run again with only one infile and one outfile or avoid specifying a name so the code can define names.")
+    logger.error("Several infiles given, but only one outfile. Run again with only one infile and one outfile or avoid specifying a name so that the code can define names.")
+
+  fileCounter = 0
+  for fileItem in files:
+
+    logger.info("Starting generation of vertical profile from polar volume: %s"%fileItem)
+
+    # Load the file item 
+    obj = None
+    rio = _raveio.open(fileItem)
+    obj = rio.object
+
+    # The input must be a polar volume, if it is not raise an exception!
+    if not _polarvolume.isPolarVolume(obj):
+      raise AttributeError("Must call wrwp_main with a polar volume as input, check polar file: %s"%fileItem)
+      logger.error("Must call wrwp_main with polar volume as input, check polar file: %s"%fileItem)
+
+    RaveIO_ODIM_Version = rio.version
+    RaveIO_ODIM_H5rad_Version = rio.h5radversion
+    rio.close()
+
+    if RaveIO_ODIM_Version == 1 and RaveIO_ODIM_H5rad_Version == 1:
+      options.setOdim21 = True
+      logger.info("Input volume according to ODIM-H5 ver2.1 implies an output vertical profile with ver2.1")
+
+    if fileCounter == 0:
+      if options.setOdim21:
+        print("\nGenerated vertical profile/s has ODIM version 2.1")
+      else:
+        print("\nGenerated vertical profile/s has ODIM version 2.2")
+
+      if options.outpath != None:
+        print("Vertical profile/s will be placed in directory: %s"%options.outpath)
+      else:
+        print("Vertical profiles will be placed in the wrwp_main.py install directory")
+      fileCounter = fileCounter + 1
+
+    # Define and load the wrwp object
+    wrwp = _wrwp.new()
+    fields = None
+
+    wrwp.dmin = options.dmin
+    wrwp.dmax = options.dmax
+    wrwp.nmin_wnd = options.nmin_wnd
+    wrwp.nmin_ref = options.nmin_ref
+    wrwp.emin = options.emin
+    wrwp.emax = options.emax
+    wrwp.vmin = options.vmin
+    wrwp.ff_max = options.ff_max
+    wrwp.dz = options.dz
+    wrwp.hmax = options.hmax
+    wrwp.nodata_VP = options.nodata_VP
+    wrwp.undetect_VP = options.undetect_VP
+    wrwp.gain_VP = options.gain_VP
+    wrwp.offset_VP = options.offset_VP
+
+    if options.quantities != None:
+      fields = options.quantities
+    exestart = time.time()
+
+    try:
+      profile = wrwp.generate(obj, fields)
+
+      # Extract parameters for use in filename construction if the user has not set a filename explicitly.
+      product_vp = profile.product
+      date_vp = profile.date
+      time_vp = profile.time
+      source_vp = profile.source
+
+      source_split = [str(x) for x in source_vp.split(",")]
+      for i in range(len(source_split)):
+        if source_split[i].startswith("NOD"):
+          NOD = [str(x) for x in source_split[i].split(":")]
+  
+      rio = _raveio.new()
+      rio.object = profile
+
+      # Filename based on data from the generated profile
+      fname = str(NOD[1]) + "_" + str(product_vp).lower() + "_" + str(date_vp) + "T" + str(time_vp) + "Z.h5"
+  
+      if options.outfname == None:
+        rio.filename = options.outpath + fname
+      else:
+        rio.filename = options.outpath + options.outfname
+    except:
+      logger.info("No vertical profile could be generated from polar volume: %s, check input volume and wrwp parameter settings\n"%fileItem)
+      print("No vertical profile could be generated from polar volume: %s, check input volume and wrwp parameter settings\n"%fileItem)
+      return None
+    
+    filenameV22 = rio.filename
+    rio.save()
+
+    exetime1 = time.time() - exestart
+     
+    if not options.setOdim21: 
+      if os.path.isfile(filenameV22):
+        logger.debug("Succeded in writing generated ver 2.2 vertical profile to disk: " + filenameV22)
+        logger.debug("Total generation time for ver 2.2 vertical profile: %f s"%exetime1)
+        logger.info("Finished generating ver 2.2 vertical profile, file: %s\n"%filenameV22)
+      else:
+        logger.info("Failed to write ver2.2 vertical profile to disk: \n" + filenameV22)
+    else:
+      # Converts the vertical profile from ODIM-H5 v2.2 to ODIM-H5 2.1 if wanted
+      change2Odim21(filenameV22, options.quantities, options, exestart)
+
+
+if __name__ == "__main__":
+
+  # Finds the wrwp_config.xml under /local_disk/baltrad/blt_sys/ and sets parameters necessary for wrwp generation.
+  # Usually the installation paths in the node-installer is not changed relative to default, the use of find locates the wrwp_config.xml
+  # should another path be set in the node-installer, the only requirement is that the wrwp_config.xml is somewhere under /local_disk/baltrad/blt_sys.
+
+  path2config = find('wrwp_config.xml', '/local_disk/baltrad/blt_sys')
+
+  root = ET.parse(str(path2config[0])).getroot()
+
+  for param in root.findall('param'):
+    if param.get('name') == 'DMIN':
+      DMIN = param.find('value').text
+    if param.get('name') == 'DMAX':
+      DMAX = param.find('value').text
+    if param.get('name') == 'NMIN_WND':
+      NMIN_WND = param.find('value').text
+    if param.get('name') == 'NMIN_REF':
+      NMIN_REF = param.find('value').text
+    if param.get('name') == 'EMIN':
+      EMIN = param.find('value').text
+    if param.get('name') == 'EMAX':
+      EMAX = param.find('value').text
+    if param.get('name') == 'VMIN':
+      VMIN = param.find('value').text
+    if param.get('name') == 'FF_MAX':
+      FF_MAX = param.find('value').text
+    if param.get('name') == 'DZ':
+      DZ = param.find('value').text
+    if param.get('name') == 'HMAX':
+      HMAX = param.find('value').text
+    if param.get('name') == 'NODATA_VP':
+      NODATA_VP = param.find('value').text
+    if param.get('name') == 'UNDETECT_VP':
+      UNDETECT_VP = param.find('value').text
+    if param.get('name') == 'GAIN_VP':
+      GAIN_VP = param.find('value').text
+    if param.get('name') == 'OFFSET_VP':
+      OFFSET_VP = param.find('value').text
+    if param.get('name') == 'QUANTITIES':
+      QUANTITIES = param.find('value').text
+
+    QUANTITIES_DEF = 'NV,HGHT,UWND,VWND,ff,ff_dev,dd,DBZH,DBZH_dev,NZ'
+
+  usage = "usage: %wrwp_main --infile <infile/s> --outpath <path to profiles> [args] [h]"
+  usage += "\nGenerates weather radar wind profiles directly from polar volumes."
+  usage += "\nIf a ver2.1 input volume is used, the output vertical profile is ver2.1."
+  usage += "\nIf a ver2.2 input volume is used, the output vertical profile can be either ver2.1 and ver2.2."
+  usage += "\nAdjustable parameters are stored in wrwp_config.xml but can of course also be changed from the command line."
+  usage += "\nThe default install directory for the wrwp_config.xml is under .../baltrad-wrwp/config/"
+  usage += "\nThe script is equipped with a non-rotating log, which ends up in the same directory as wrwp_main"
+
+  parser = OptionParser(usage=usage)
+
+  parser.add_option("--infiles", dest = "infiles", help = "Name of polar volume input files (including path) " + \
+                    "from which to create a vertical profile, must be in ODIM-H5 format. Note that if only one " + \
+                    "file is given it is possible to select a filename, if no filename is specified, the code " + \
+                    "will generate a name based on NOD, product, date and time. If multiple files are given no filenames" + \
+                    "can be specified and the code generates names based on earlier mentioned attributes. " + \
+                    "If several files are given, they must be given comma-separated without space in between.")
+  parser.add_option("--outpath", dest = "outpath", default = DEFAULTOUT, help = "The path to the directory " + \
+                    "where the outfile/s is/are written, default installdir")
+  parser.add_option("--outfname", dest = "outfname", default = None, help = "Selected filename for vertical profile, in ODIM-H5 format, " + \
+                    "in the case that only one input file is given, default None.")
+  parser.add_option("--quantities", dest = "quantities", type = "string", default = QUANTITIES_DEF, help = "Comma separated list of " + \
+                    "quantities. Currently supported quantities are NV,HGHT,UWND,VWND,ff,ff_dev,dd,DBZH,DBZH_dev,NZ and are given " + \
+                    "in the wrwp_config.xml file. Default if no quantities are given is UWND, VWND, HGHT, ff, ff_dev, dd, NV, DBZH, DBZH_dev and NZ.")
+  parser.add_option("--dmin", dest = "dmin", type = "int", default = DMIN, help="Minimum distance for deriving a profile [m], default 5000. Note: integer.")
+  parser.add_option("--dmax", dest = "dmax", type = "int", default = DMAX, help="Maximum distance for deriving a profile [m], default 25000. Note: integer.")
+  parser.add_option("--nmin_wnd", dest = "nmin_wnd", type = "int", default = NMIN_WND, help="Minimum sample size for wind, default 40. Note: integer.")
+  parser.add_option("--nmin_ref", dest = "nmin_ref", type = "int", default = NMIN_REF, help="Minimum sample size for reflectivity, default 40. Note: integer.")
+  parser.add_option("--emin", dest = "emin", type = "float", default = EMIN, help="Minimum elevation angle [deg], default 4.0.")
+  parser.add_option("--emax", dest = "emax", type = "float", default = EMAX, help="Maximum elevation angle [deg], default 45.0.")
+  parser.add_option("--vmin", dest = "vmin", type = "float", default = VMIN, help="Radial velocity threshold [m/s], default 2.0")
+  parser.add_option("--ff_max", dest = "ff_max", type = "float", default = FF_MAX, help="Maximum allowed calculated " + \
+                    "layer velocity [m/s], default 60.0. Layer velocity greater than ff_max will be set as nodata. " + \
+                    "Note that this implies also setting the remaining wind related parameters for this layer as nodata")
+  parser.add_option("--dz", dest = "dz", type = "int", default = DZ, help="Height interval for the generated vertical profile [m], default 200. Note: integer.")
+  parser.add_option("--hmax", dest = "hmax", type = "int", default = HMAX, help="Maximum height of the generated vertical profile [m], default 12000. Note:integer.")
+  parser.add_option("--nodata_VP", dest = "nodata_VP", type = "int", default = NODATA_VP, help="Nodata value for vertical profile, default -9999. Note: integer.")
+  parser.add_option("--undetect_VP", dest = "undetect_VP", type = "int", default = UNDETECT_VP, help="Undetect value for vertical profile, default -9999. Note:integer.")
+  parser.add_option("--gain_VP", dest = "gain_VP", type = "float", default = GAIN_VP, help="Gain value for vertical profile, default 1.0.")
+  parser.add_option("--offset_VP", dest = "offset_VP", type = "float", default = OFFSET_VP, help="Offset value for vertical profile, default 0.0.")
+  parser.add_option("--setOdim21", dest = "setOdim21", action="store_true", default = False, help="Converts the VP to ver2.1 if set, default False.")
+  parser.add_option("--verbose", dest = "verbose", action="store_true", default = False, help="Enables verbose logging and verbose printing of some info to the terminal, default False.")
+   
+  (options, args) = parser.parse_args()
+
+  # Putting a file extension to the selected filename just in case...
+  if options.outfname != None and not options.outfname.endswith(".h5"):
+    options.outfname = options.outfname + ".h5"
+  if options.outpath != DEFAULTOUT:
+    # Create main output directory if it does not exist
+    if not os.path.exists(options.outpath):
+      os.makedirs(options.outpath)
+
+    # Put a training slash if does not exist already
+    if not options.outpath.endswith('/'):
+      options.outpath = options.outpath + '/'
+
+  # Set up the logging system
+  logger = logging.getLogger('wrwp_log')
+  handler = logging.FileHandler(DEFAULTOUT + '/wrwp_log.log')
+  formatter = logging.Formatter('%(asctime)s %(levelname)s %(message)s')
+  handler.setFormatter(formatter)
+  logger.addHandler(handler)
+
+  if options.verbose:
+    log_level = logging.DEBUG
+  else:
+    log_level = logging.INFO
+
+  logger.setLevel(log_level)
+
+
+  # Setting two "flags" responsible for filtering some printed parameter output
+  refl_flag = False
+  wnd_flag = False
+
+  if options.infiles != None:
+    filelist = [str(x) for x in options.infiles.split(",")] 
+    print("\nGenerating vertical profile/s from volume/s:")
+    for i in range(0, len(filelist)):
+      print(filelist[i])
+  
+  if options.verbose:
+    print("\nSelected vertical profile parameter settings")
+    print("----------------------------------------------")
+
+    if options.quantities != None:
+      print("Quantities in the generated vertical profile: %s"%options.quantities)
+      if "DBZH" in options.quantities or "DBZH_dev" in options.quantities or "NZ" in options.quantities:
+        refl_flag = True
+      if "ff" in options.quantities or "ff_dev" in options.quantities or "dd" in options.quantities or \
+        "UWND" in options.quantities or "VWND" in options.quantities or "NV" in options.quantities:
+        wnd_flag = True
+    else:
+      print("Quantities in the generated vertical profile: ff, ff_dev, dd, NV, DBZH, DBZH_dev, NZ")
+      refl_flag = True
+      wnd_flag = True
+
+    print("Minimum distance for vertical profile, dmin: %s [m]"%options.dmin)
+    print("Maximum distance for vertical profile, dmax: %s [m]"%options.dmax)
+    print("Minimum elevation angle used in vertical profile, emin: %s [deg]"%options.emin)
+    print("Maximum elevation angle used in the vertical profile, emax: %s [deg]"%options.emax)
+    print("Height interval for vertical profile, dz: %s [m]"%options.dz)
+    print("Maximum height for vertical profile, hmax: %s [m]"%options.hmax)
+
+    if wnd_flag:
+      print("Minimum sample size required for wind in the vertical profile, nmin_wnd: %s"%options.nmin_wnd)
+      print("Radial velocity threshold for vertical profile, vmin: %s [m/s]"%options.vmin)
+      print("Maximum allowed layer velocity for vertical profile, ff_max: %s [m/s]"%options.ff_max)
+    if refl_flag:
+      print("Minimum sample size required for reflectivity in the vertical profile, nmin_ref: %s"%options.nmin_ref)
+
+    print("Nodata value, nodata_VP: %s"%options.nodata_VP)
+    print("Undetect value, undetect_VP: %s"%options.undetect_VP)
+    print("Gain value, gain_VP: %s"%options.gain_VP)
+    print("Offset value, offset_VP: %s"%options.offset_VP)
+
+  if options.infiles != None:
+    main(options)
+  else:
+    parser.print_help()
+  
diff --git a/bin/wrwp_main.c b/bin/wrwp_main.c
deleted file mode 100644 (file)
index 66cf64e..0000000
+++ /dev/null
@@ -1,397 +0,0 @@
-/* --------------------------------------------------------------------
-Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
-
-This is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-This software is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with baltrad-wrwp.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------*/
-#include "wrwp.h"
-#include "rave_debug.h"
-#include <getopt.h>
-#include <libgen.h>
-#include "rave_utilities.h"
-
-static void PrintUsage(char* name, int full)
-{
-  char* namecopy = RAVE_STRDUP(name);
-
-  printf("Usage: %s [options] <input volume.h5> <output verticalprofile.h5>\n", basename(namecopy));
-
-  if (full) {
-    printf("--help             - Prints this output\n");
-    printf("--verbose          - Produces some information about the generated product\n");
-    printf("--debug            - Produces some debug information during the generation\n");
-    printf("--dz=<value>       - Height interval for deriving a profile [m] (default: %d)\n", DZ);
-    printf("--nodata=<value>   - Nodata value (default: %d)\n", NODATA_VP);
-    printf("--undetect=<value> - Undetect value (default: %d)\n", UNDETECT_VP);
-    printf("--gain=<value>     - Gain value (default: %f)\n", GAIN_VP);
-    printf("--offset=<value>   - Offset value (default: %f)\n", OFFSET_VP);
-    printf("--hmax=<value>     - Maximum height of the profile [m] (default: %d)\n", HMAX);
-    printf("--dmin=<value>     - Minimum distance for deriving a profile [m] (default: %d)\n", DMIN);
-    printf("--dmax=<value>     - Maximum distance for deriving a profile [m] (default: %d)\n", DMAX);
-    printf("--emin=<value>     - Minimum elevation angle [deg] (default: %f)\n", EMIN);
-    printf("--vmin=<value>     - Radial velocity threshold [m/s] (default: %f)\n", VMIN);
-    printf("--quantities=<value> A comma separated list of quanities (default: ff,ff_dev,dd,dbzh,dbzh_dev,nz)\n");
-    printf("                     Currently supported quantities are: NV,HGHT,UWND,VWND,ff,ff_dev,dd,dbzh,dbzh_dev,nz\n");
-    printf("\n");
-    printf("<input volume.h>  must be a polar volume in ODIM H5 format\n");
-    printf("<output verticalprofile.h5> will be a vertical profile in ODIM H5 format\n\n");
-  }
-
-  RAVE_FREE(namecopy);
-}
-
-/**
- * Tries to translate an string into an int value.
- * @param[in] arg - the argument
- * @param[in,out] pInt - the outvalue
- * @return 1 on success otherwise 0
- */
-int ParseInt(char* arg, int* pInt)
-{
-  if (arg != NULL) {
-    int i = 0;
-    while (arg[i] != '\0') {
-      if (arg[i] < '0' || arg[i] > '9') {
-        return 0;
-      }
-      i++;
-    }
-  }
-  if (arg != NULL && pInt != NULL && sscanf(arg, "%d",pInt) == 1) {
-    return 1;
-  }
-  return 0;
-}
-
-/**
- * Tries to translate an string into an int value.
- * @param[in] arg - the argument
- * @param[in,out] pDouble - the outvalue
- * @return 1 on success otherwise 0
- */
-int ParseDouble(char* arg, double* pDouble)
-{
-  if (arg != NULL) {
-    int i = 0;
-    int nrDots = 0;
-    while (arg[i] != '\0') {
-      if (arg[i] < '0' || arg[i] > '9') {
-        if (arg[i] == '.' && nrDots==0) {
-          nrDots++;
-        } else {
-          return 0;
-        }
-      }
-      i++;
-    }
-  }
-  if (arg != NULL && pDouble != NULL && sscanf(arg, "%lf",pDouble) == 1) {
-    return 1;
-  }
-  return 0;
-}
-
-/**
- * Extracts the quantity list from a string. Verifies that it's a comma separated list
- * @param[in] arg - the argument
- * @param[in,out] quantities - the list of quantities
- * @return 1 on success otherwise 0
- */
-int ParseQuantities(char* arg, char** quantities)
-{
-  int result = 0;
-  RaveList_t* quantityList = NULL;
-  if (arg != NULL) {
-    quantityList = RaveUtilities_getTrimmedTokens(arg, ',');
-    if (quantityList == NULL || RaveList_size(quantityList) == 0) { /* When 0 length, we let wrwp:s default mechanism decide */
-      *quantities = NULL;
-    } else {
-      int i = 0, sz = 0;
-      sz = RaveList_size(quantityList);
-      for (i = 0; i < sz; i++) {
-        char* q = (char*)RaveList_get(quantityList, i);
-        if (q != NULL) {
-          if (strcmp(q, "NV") == 0 ||
-              strcmp(q, "HGHT") == 0 ||
-              strcmp(q, "UWND") == 0 ||
-              strcmp(q, "VWND") == 0 ||
-              strcmp(q, "ff") == 0 ||
-              strcmp(q, "ff_dev") == 0 ||
-              strcmp(q, "dd") == 0 ||
-              strcmp(q, "dbzh") == 0 ||
-              strcmp(q, "dbzh_dev") == 0 ||
-              strcmp(q, "nz") == 0) {
-            continue;
-          } else {
-            fprintf(stderr, "quantity = %s is not supported\n", q);
-            goto done;
-          }
-        }
-      }
-      *quantities = RAVE_STRDUP(arg);
-      if (*quantities == NULL) {
-        fprintf(stderr, "Failed to allocate memory for quantities\n");
-        goto done;
-      }
-    }
-  }
-  result = 1;
-done:
-  if (quantityList != NULL) {
-    RaveList_freeAndDestroy(&quantityList);
-  }
-  return result;
-}
-
-
-/** Main function for deriving weather radar wind and reflectivity profiles
- * @file
- * @author G´┐Żnther Haase, SMHI
- * @date 2011-11-29
- */
-int main (int argc,char *argv[]) {
-       RaveIO_t* raveio = NULL;
-       PolarVolume_t* inobj = NULL;
-       VerticalProfile_t* result = NULL;
-       Wrwp_t* wrwp = NULL;
-       int exitcode = 127;
-       char* inputfile = NULL; /* do not delete */
-       char* outputfile = NULL; /* do not delete */
-
-       int verbose_flag=0;
-       int debug_flag=0;
-       int help_flag=0;
-       int dz_value = DZ;
-       int nodata_vp = NODATA_VP;
-       int undetect_vp = UNDETECT_VP;
-       double gain_vp = GAIN_VP;
-       double offset_vp = OFFSET_VP;
-       int hmax = HMAX;
-       int dmin = DMIN;
-       int dmax = DMAX;
-       double emin = EMIN;
-       double vmin = VMIN;
-       char* quantities = NULL;
-
-       int getopt_ret, option_index;
-
-       struct option long_options[] = {
-           {"help", no_argument, &help_flag, 1},
-           {"verbose", no_argument, &verbose_flag, 1},
-      {"debug", no_argument, &debug_flag, 1},
-           {"dz", optional_argument, 0, 'z'},
-           {"nodata", optional_argument, 0, 'n'},
-           {"undetect", optional_argument, 0, 'u'},
-      {"gain", optional_argument, 0, 'g'},
-      {"offset", optional_argument, 0, 'o'},
-      {"hmax", optional_argument, 0, 'h'},
-      {"dmin", optional_argument, 0, 'd'},
-      {"dmax", optional_argument, 0, 'D'},
-      {"emin", optional_argument, 0, 'e'},
-      {"vmin", optional_argument, 0, 'v'},
-      {"quantities", optional_argument, 0, 'q'},
-           {0, 0, 0, 0}
-       };
-
-
-       Rave_initializeDebugger();
-       Rave_setDebugLevel(RAVE_INFO);
-
-  while (1) {
-      getopt_ret = getopt_long( argc, argv, "",
-                                long_options,  &option_index);
-      if (getopt_ret == -1) break;
-      if (help_flag) {
-        PrintUsage(argv[0], 1);
-        exitcode=1;
-        goto done;
-      }
-      if (verbose_flag) {
-        Rave_setDebugLevel(RAVE_DEBUG);
-      }
-
-      switch(getopt_ret) {
-      case 0:
-        /* Here all non named options will arrive */
-        /*printf("option %s\n", long_options[option_index].name);*/
-        break;
-      case 'z':
-        if (!ParseInt(optarg, &dz_value)) {
-          fprintf(stderr, "--dz=<value> must be an integer value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'n':
-        if (!ParseInt(optarg, &nodata_vp)) {
-          fprintf(stderr, "--nodata=<value> must be an integer value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'u':
-        if (!ParseInt(optarg, &undetect_vp)) {
-          fprintf(stderr, "--undetect=<value> must be an integer value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'g':
-        if (!ParseDouble(optarg, &gain_vp)) {
-          fprintf(stderr, "--gain=<value> must be a double value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'o':
-        if (!ParseDouble(optarg, &offset_vp)) {
-          fprintf(stderr, "--offset=<value> must be a double value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'h':
-        if (!ParseInt(optarg, &hmax)) {
-          fprintf(stderr, "--hmax=<value> must be an integer value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'd':
-        if (!ParseInt(optarg, &dmin)) {
-          fprintf(stderr, "--dmin=<value> must be an integer value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'D':
-        if (!ParseInt(optarg, &dmax)) {
-          fprintf(stderr, "--dmax=<value> must be an integer value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'e':
-        if (!ParseDouble(optarg, &emin)) {
-          fprintf(stderr, "--emin=<value> must be a double value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'v':
-        if (!ParseDouble(optarg, &vmin)) {
-          fprintf(stderr, "--vmin=<value> must be a double value\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case 'q':
-        if (!ParseQuantities(optarg, &quantities)) {
-          fprintf(stderr, "--quantities=<value> must be a comma separated list of valid quantities.\n");
-          PrintUsage(argv[0], 0);
-          goto done;
-        }
-        break;
-      case '?':
-      default:
-        fprintf(stderr, "Unknown argument\n");
-        PrintUsage(argv[0], 0);
-        break;
-      }
-  }
-
-  if (argc - optind != 2) {
-    PrintUsage(argv[0], 0);
-    goto done;
-  }
-
-  inputfile=argv[optind++];
-  outputfile=argv[optind++];
-
-  if (argc<3) {
-               printf ("Usage: %s <input ODIM_H5 polar volume> <output ODIM_H5 polar volume> \n",argv[0]);
-               exit (1);
-       }
-
-       raveio = RaveIO_open(inputfile);
-       if (raveio == NULL) {
-         fprintf(stderr, "Failed to open file = %s\n", inputfile);
-         goto done;
-       }
-
-       wrwp = RAVE_OBJECT_NEW(&Wrwp_TYPE);
-       if (wrwp == NULL) {
-         fprintf(stderr, "Failed to create wrwp object\n");
-         goto done;
-       }
-
-       Wrwp_setDZ(wrwp, dz_value);
-       Wrwp_setNODATA_VP(wrwp, nodata_vp);
-       Wrwp_setUNDETECT_VP(wrwp, undetect_vp);
-       Wrwp_setOFFSET_VP(wrwp, offset_vp);
-       Wrwp_setGAIN_VP(wrwp, gain_vp);
-       Wrwp_setHMAX(wrwp, hmax);
-       Wrwp_setDMIN(wrwp, dmin);
-       Wrwp_setDMAX(wrwp, dmax);
-       Wrwp_setEMIN(wrwp, emin);
-       Wrwp_setVMIN(wrwp, vmin);
-
-       /*Opening of HDF5 radar input file.*/
-       if (RaveIO_getObjectType(raveio)== Rave_ObjectType_PVOL) {
-               inobj = (PolarVolume_t*)RaveIO_getObject(raveio);
-       }
-       else {
-               printf ("Input file is not a polar volume. Giving up ...\n");
-               goto done;
-       }
-       RaveIO_close (raveio);
-
-       result = Wrwp_generate(wrwp, inobj, quantities);
-       if (inobj == NULL) {
-               printf ("Could not derive wind profile %s, exiting ...\n", argv[1]);
-               goto done;
-       }
-
-       RaveIO_setObject(raveio, (RaveCoreObject*)result);
-
-       if (!RaveIO_save(raveio, outputfile))
-         goto done;
-
-  if (debug_flag) {
-    printf("Generated vertical profile...\n");
-    printf("Input file: %s\n", inputfile);
-    printf("Output file: %s\n", outputfile);
-    printf("DZ         = %d\n", Wrwp_getDZ(wrwp));
-    printf("NODATA     = %d\n", Wrwp_getNODATA_VP(wrwp));
-    printf("UNDETECT   = %d\n", Wrwp_getUNDETECT_VP(wrwp));
-    printf("GAIN       = %lf\n", Wrwp_getGAIN_VP(wrwp));
-    printf("OFFSET     = %lf\n", Wrwp_getOFFSET_VP(wrwp));
-    printf("HMAX       = %d\n", Wrwp_getHMAX(wrwp));
-    printf("DMIN       = %d\n", Wrwp_getDMIN(wrwp));
-    printf("DMAX       = %d\n", Wrwp_getDMAX(wrwp));
-    printf("EMIN       = %lf\n", Wrwp_getEMIN(wrwp));
-    printf("VMIN       = %lf\n", Wrwp_getVMIN(wrwp));
-    if (quantities != NULL) {
-      printf("QUANTITIES = %s\n", quantities);
-    }
-  }
-
-       exitcode = 0;
-done:
-       RAVE_OBJECT_RELEASE(raveio);
-       RAVE_OBJECT_RELEASE(inobj);
-       RAVE_OBJECT_RELEASE(result);
-       RAVE_FREE(quantities);
-
-       return exitcode;
-}
index f0d621e..62c06fb 100644 (file)
@@ -29,9 +29,16 @@ along with baltrad-wrwp.  If not, see <http://www.gnu.org/licenses/>.
  * incoming polar volume file only contained scans with elevation angle
  * smaller than the one stated in the call to the pgf. These kind of files
  * are rather common when using the central software EDGE coming with the
- * modernized Swedish radars. If suck a file is encountered, the return is just
+ * modernized Swedish radars. If such a file is encountered, the return is just
  * NULL and the exception following is treated in the pgf.
- */
+ * 
+ * @date 2020-01-15, corrected an bug in the wind direction calculation. In addition, added functionality
+ * for the new parameters EMAX and FF_MAX plus activated the parameter NZ.
+ *
+ * The majority of the parameters put in the file wrwp.h have been moved into a separate file (wrwp_config.xml) 
+ * so that they can be changed without compiling the code. The plugin has been modified so it reads the xml-file 
+ * the relevant parameters not set from the web-GUI and a newly created python command-line tool (wrwp_main.py 
+ * instead of wrwp_main.c) also reads the xml-file and sets the relevant parameters. */
 
 #include "wrwp.h"
 #include "vertical_profile.h"
@@ -54,7 +61,11 @@ struct _Wrwp_t
   int hmax; /**< Maximum height of the profile [m] */
   int dmin; /**< Minimum distance for deriving a profile [m] */
   int dmax; /**< Maximum distance for deriving a profile [m]*/
+  int nmin_wnd; /**< Minimum sample size for wind */
+  int nmin_ref; /**< Minimum sample size for refl. */
   double emin; /**< Minimum elevation angle [deg] */
+  double emax; /**< Maximum elevation angle [deg] */
+  double ff_max; /**<Maximum accepted value for ff (calculated layer velocity) [m/s]*/
   double vmin; /**< Radial velocity threshold [m/s] */
   double nodata_VP; /**< Nodata value for vertical profile */
   double gain_VP; /**< Gain for VP fields */
@@ -69,16 +80,20 @@ struct _Wrwp_t
 static int Wrwp_constructor(RaveCoreObject* obj)
 {
   Wrwp_t* wrwp = (Wrwp_t*)obj;
-  wrwp->hmax = HMAX;
-  wrwp->dz = DZ;
-  wrwp->dmin = DMIN;
-  wrwp->dmax = DMAX;
-  wrwp->emin = EMIN;
-  wrwp->vmin = VMIN;
-  wrwp->nodata_VP = NODATA_VP;
-  wrwp->gain_VP = GAIN_VP;
-  wrwp->offset_VP = OFFSET_VP;
-  wrwp->undetect_VP = UNDETECT_VP;
+  wrwp->dz = 0;
+  wrwp->hmax = 0; 
+  wrwp->dmin = 0;
+  wrwp->dmax = 0;
+  wrwp->nmin_wnd = 0;
+  wrwp->nmin_ref = 0;
+  wrwp->emin = 0.0;
+  wrwp->emax = 0.0;
+  wrwp->ff_max = 0.0;
+  wrwp->vmin = 0.0;
+  wrwp->nodata_VP = 0;
+  wrwp->gain_VP = 1.0; /* The gain cannot be initialized to 0.0! */
+  wrwp->offset_VP = 0.0;
+  wrwp->undetect_VP = 0;
   return 1;
 }
 
@@ -89,7 +104,7 @@ static void Wrwp_destructor(RaveCoreObject* obj)
 {
 }
 
-static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t* pvol, const char* name, double minSelAng)
+static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t* pvol, const char* name, double minSelAng, double maxSelAng)
 {
   int nscans = PolarVolume_getNumberOfScans(pvol);
   int i = 0;
@@ -101,9 +116,8 @@ static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t
     if (scan != NULL && PolarScan_hasAttribute(scan, name)) {
         elangleForThisScan = PolarScan_getElangle(scan);
         
-        /* Filter with respect to the selected min elangle
-           that is given in the web GUI route for WRWP generation */
-        if ((elangleForThisScan * RAD2DEG) >= minSelAng) {
+        /* Filter with respect to the selected min elangle and max elangle*/
+        if ((elangleForThisScan * RAD2DEG) >= minSelAng && (elangleForThisScan * RAD2DEG) <= maxSelAng) {
           RaveAttribute_t* attr = PolarScan_getAttribute(scan, name);
           VerticalProfile_addAttribute(vp, attr);
           found = 1;
@@ -126,6 +140,27 @@ static int WrwpInternal_addDoubleAttribute(VerticalProfile_t* vp, const char* na
   return result;
 }
 
+static int WrwpInternal_addStringAttribute(VerticalProfile_t* vp, const char* name, const char* value)
+{
+  RaveAttribute_t* attr = RaveAttributeHelp_createString(name, value);
+  int result = 0;
+  if (attr != NULL) {
+    result = VerticalProfile_addAttribute(vp, attr);
+  }
+  RAVE_OBJECT_RELEASE(attr);
+  return result;
+}
+static char* WrwpInternal_LastcharDel(char* name)
+{
+  int i = 0;
+  while(name[i] != '\0')
+  {
+  i++;       
+  }
+  name[i-1] = '\0';
+  return name;
+}
 
 /* Function that adds various quantities under a field's what in order to
    better resemble the function existing in vertical profiles from N2 */
@@ -155,9 +190,10 @@ static RaveList_t* WrwpInternal_createFieldsList(const char* fieldsToGenerate)
     RaveList_add(result, RAVE_STRDUP("ff"));
     RaveList_add(result, RAVE_STRDUP("ff_dev"));
     RaveList_add(result, RAVE_STRDUP("dd"));
-    RaveList_add(result, RAVE_STRDUP("dbzh"));
-    RaveList_add(result, RAVE_STRDUP("dbzh_dev"));
-    RaveList_add(result, RAVE_STRDUP("nz"));
+    RaveList_add(result, RAVE_STRDUP("NV"));
+    RaveList_add(result, RAVE_STRDUP("DBZH"));
+    RaveList_add(result, RAVE_STRDUP("DBZH_dev"));
+    RaveList_add(result, RAVE_STRDUP("NZ"));
   }
   return result;
 }
@@ -179,7 +215,7 @@ static int WrwpInternal_containsField(RaveList_t* fieldIds, const char* id)
   }
   return 0;
 }
-
+/* Adds attributes under a fileds what */
 static int WrwpInternal_addNodataUndetectGainOffset(RaveField_t* field, double nodata, double undetect, double gain, double offset)
 {
   if (!WrwpInternal_addDoubleAttr2Field(field, "what/nodata", nodata) ||
@@ -191,6 +227,7 @@ static int WrwpInternal_addNodataUndetectGainOffset(RaveField_t* field, double n
   return 1;
 }
 
+/* Gathering of startdatetime and enddatetime */
 static RaveDateTime_t* WrwpInternal_getStartDateTimeFromScan(PolarScan_t* scan)
 {
   RaveDateTime_t* dt = RAVE_OBJECT_NEW(&RaveDateTime_TYPE);
@@ -328,6 +365,29 @@ int Wrwp_getDMAX(Wrwp_t* self)
   return self->dmax;
 }
 
+void Wrwp_setNMIN_WND(Wrwp_t* self, int nmin_wnd)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  self->nmin_wnd = nmin_wnd;
+}
+
+int Wrwp_getNMIN_WND(Wrwp_t* self)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  return self->nmin_wnd;
+}
+
+void Wrwp_setNMIN_REF(Wrwp_t* self, int nmin_ref)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  self->nmin_ref = nmin_ref;
+}
+
+int Wrwp_getNMIN_REF(Wrwp_t* self)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  return self->nmin_ref;
+}
 void Wrwp_setEMIN(Wrwp_t* self, double emin)
 {
   RAVE_ASSERT((self != NULL), "self == NULL");
@@ -340,6 +400,30 @@ double Wrwp_getEMIN(Wrwp_t* self)
   return self->emin;
 }
 
+void Wrwp_setEMAX(Wrwp_t* self, double emax)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  self->emax = emax;
+}
+
+double Wrwp_getEMAX(Wrwp_t* self)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  return self->emax;
+}
+
+void Wrwp_setFF_MAX(Wrwp_t* self, double ff_max)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  self->ff_max = ff_max;
+}
+
+double Wrwp_getFF_MAX(Wrwp_t* self)
+{
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  return self->ff_max;
+}
+
 void Wrwp_setVMIN(Wrwp_t* self, double vmin)
 {
   RAVE_ASSERT((self != NULL), "self == NULL");
@@ -352,6 +436,7 @@ double Wrwp_getVMIN(Wrwp_t* self)
   return self->vmin;
 }
 
+/* Main code for vertical profile generation */
 VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char* fieldsToGenerate)
 {
   VerticalProfile_t* result = NULL;
@@ -366,12 +451,13 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
   double centerOfLayer=0.0, u_wnd_comp=0.0, v_wnd_comp=0.0, vdir_rad=0.0;
   int ysize = 0, yindex = 0;
   int countAcceptedScans = 0; /* counter for accepted scans i.e. scans with elangle >= selected
-                                 minimum elevatiuon angle and not being set as malfunc */
+                                 minimum elevatiuon angle and <= selected maximum elevation angle and not being set as malfunc */
 
   const char* product = "VP";
+
   RaveDateTime_t *firstStartDT = NULL, *lastEndDT = NULL;
 
-  /* Field defs */
+  /* Field definitions */
   RaveField_t *nv_field = NULL, *hght_field = NULL;
   RaveField_t *uwnd_field = NULL, *vwnd_field = NULL;
   RaveField_t *ff_field = NULL, *ff_dev_field = NULL, *dd_field = NULL;
@@ -391,9 +477,9 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
   if (WrwpInternal_containsField(wantedFields, "ff")) ff_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
   if (WrwpInternal_containsField(wantedFields, "ff_dev")) ff_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
   if (WrwpInternal_containsField(wantedFields, "dd")) dd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
-  if (WrwpInternal_containsField(wantedFields, "dbzh")) dbzh_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
-  if (WrwpInternal_containsField(wantedFields, "dbzh_dev")) dbzh_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
-  if (WrwpInternal_containsField(wantedFields, "nz")) nz_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
+  if (WrwpInternal_containsField(wantedFields, "DBZH")) dbzh_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
+  if (WrwpInternal_containsField(wantedFields, "DBZH_dev")) dbzh_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
+  if (WrwpInternal_containsField(wantedFields, "NZ")) nz_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 
   if ((WrwpInternal_containsField(wantedFields, "NV") && nv_field == NULL) ||
       (WrwpInternal_containsField(wantedFields, "HGHT") && hght_field == NULL) ||
@@ -402,14 +488,15 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
       (WrwpInternal_containsField(wantedFields, "ff") && ff_field == NULL) ||
       (WrwpInternal_containsField(wantedFields, "ff_dev") && ff_dev_field == NULL) ||
       (WrwpInternal_containsField(wantedFields, "dd") && dd_field == NULL) ||
-      (WrwpInternal_containsField(wantedFields, "dbzh") && dbzh_field == NULL) ||
-      (WrwpInternal_containsField(wantedFields, "dbzh_dev") && dbzh_dev_field == NULL) ||
-      (WrwpInternal_containsField(wantedFields, "nz") && nz_field == NULL))
+      (WrwpInternal_containsField(wantedFields, "DBZH") && dbzh_field == NULL) ||
+      (WrwpInternal_containsField(wantedFields, "DBZH_dev") && dbzh_dev_field == NULL) ||
+      (WrwpInternal_containsField(wantedFields, "NZ") && nz_field == NULL))
   {
     RAVE_ERROR0("Failed to allocate memory for the resulting vp fields");
     goto done;
   }
 
+  /* Set the spacing */
   ysize = self->hmax / self->dz;
 
   if ((nv_field != NULL && !RaveField_createData(nv_field, 1, ysize, RaveDataType_INT)) ||
@@ -436,7 +523,22 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
   // We use yindex for filling in the arrays even though we loop to hmax...
   yindex = 0;
 
-  // loop over atmospheric layers
+  // Allocate memory, and initialize with zeros, for the char array holding the accepted elevation angles.
+  // Define and initialize also the accessory strings and the counter.
+  char *theUsedElevationAngles = RAVE_CALLOC((size_t)(100), sizeof (char));
+  char angle[6] = {'\0'};
+  char comma[2] = ",";
+  int firstAngleCounter = 0;
+  double acceptedAngle = 0.0;
+  
+  // Some definitions  and initializations of things used for putting a unique set of tasks into an array
+  #define len(A) (sizeof(A) / sizeof(A[0]))
+  char *taskArgs[] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
+  int scanNumber = 0;
+  int foundTask = 0;
+  int ntask = 0;
+
+  // Loop over the atmospheric layers
   for (iz = 0; iz < self->hmax; iz += self->dz) {
     /* allocate memory and initialize with zeros */
     double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
@@ -454,30 +556,69 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
     nv = 0;
     nz = 0;
     
-    /* Define the center height of each verical layer, this will later
+    /* Define the center height of each vertical layer, this will later
        become the HGHT array */
     centerOfLayer = iz + (self->dz / 2.0);
     firstInit = 0;
 
-    // loop over scans
+    // Loop over scans in the polar volume
     for (is = 0; is < nscans; is++) {
       char* malfuncString = NULL;
+      char* taskString = NULL;
       PolarScan_t* scan = PolarVolume_getScan(inobj, is);      
       long nbins = PolarScan_getNbins(scan);
       long nrays = PolarScan_getNrays(scan);
       double rscale = PolarScan_getRscale(scan);
       double elangleForThisScan = PolarScan_getElangle(scan);
       
-      if (elangleForThisScan * RAD2DEG >= self->emin) { /* We only do the calculation for scans with elangle >= the minimum one */
+      if (elangleForThisScan * RAD2DEG >= self->emin && elangleForThisScan * RAD2DEG <= self->emax) { /* We only do the calculation for scans with elangle >= the minimum one AND elangle <= the maximum one */
         RaveDateTime_t* startDTofThisScan = WrwpInternal_getStartDateTimeFromScan(scan);
         RaveDateTime_t* endDTofThisScan = WrwpInternal_getEndDateTimeFromScan(scan);        
         RaveAttribute_t* malfuncattr = PolarScan_getAttribute(scan, "how/malfunc");
+        RaveAttribute_t* taskattr = PolarScan_getAttribute(scan, "how/task");
+
         if (malfuncattr != NULL) {
           RaveAttribute_getString(malfuncattr, &malfuncString); /* Set the malfuncString if attr is not NULL */
           RAVE_OBJECT_RELEASE(malfuncattr);
         }
+
+        if (taskattr != NULL) {
+          RaveAttribute_getString(taskattr, &taskString); /* Set the taskString if attr is not NULL */
+          RAVE_OBJECT_RELEASE(taskattr);
+        }
+
         if (malfuncString == NULL || strcmp(malfuncString, "False") == 0) { /* Assuming malfuncString = NULL means no malfunc */
           countAcceptedScans = countAcceptedScans + 1;
+
+          if (iz == 0) { /* Collect the elevation angles only for the first atmospheric layer */
+            acceptedAngle = elangleForThisScan * RAD2DEG;
+            sprintf(angle, "%2.1f", acceptedAngle);
+
+            if (firstAngleCounter == 0) {
+              strcat(theUsedElevationAngles, angle);
+              firstAngleCounter = 1;
+              taskArgs[scanNumber] = taskString;
+              scanNumber = scanNumber + 1; 
+            }
+            else {
+              strcat(theUsedElevationAngles, comma);
+              strcat(theUsedElevationAngles, angle);
+              for (ntask=0; ntask < len(taskArgs); ntask++) {
+                if (taskArgs[ntask]) {
+                  if (strcmp(taskArgs[ntask], taskString) == 0) {
+                    foundTask = 1;
+                    break;
+                  }
+                }
+              }
+              if (!foundTask) {
+                taskArgs[scanNumber] = taskString;
+                foundTask = 0; 
+              } 
+              scanNumber = scanNumber + 1;
+            }
+          }
+   
           if (firstInit == 0 && iz==0) { /* We only perform the date time calc at first iz-iteration*/
             /* Initialize using the first accepted scan and define 2 strings for the combined datetime */
             firstStartDT = RAVE_OBJECT_COPY(startDTofThisScan);
@@ -520,13 +661,17 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
                     (val != nodata) &&
                     (val != undetect) &&
                     (abs(offset + gain * val) >= self->vmin)) {
-                  *(v+nv) = offset+gain*val;
-                  *(az+nv) = 360./nrays*ir*DEG2RAD;
-                  *(A+nv*NOC) = sin(*(az+nv));
-                  *(A+nv*NOC+1) = cos(*(az+nv));
-                  *(A+nv*NOC+2) = 1;
-                  *(b+nv) = *(v+nv);
-                  nv = nv+1;
+                  if (nv < NOR) {
+                    *(v+nv) = offset+gain*val;
+                    *(az+nv) = 360./nrays*ir*DEG2RAD;
+                    *(A+nv*NOC) = sin(*(az+nv));
+                    *(A+nv*NOC+1) = cos(*(az+nv));
+                    *(A+nv*NOC+2) = 1;
+                    *(b+nv) = *(v+nv);
+                    nv = nv+1;
+                  } else {
+                    RAVE_ERROR0("NV to great, ignoring value");
+                  }
                 }
               }
             }
@@ -551,9 +696,13 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
                     (d <= self->dmax) &&
                     (val != nodata) &&
                     (val != undetect)) {
-                 *(z+nz) = dBZ2Z(offset+gain*val);
-                 zsum = zsum + *(z+nz);
-                nz = nz+1;
+                  if (nz < NOR) {
+                    *(z+nz) = dBZ2Z(offset+gain*val);
+                    zsum = zsum + *(z+nz);
+                    nz = nz+1;
+                  } else {
+                    RAVE_ERROR0("NZ to great, ignoring value");
+                  }
                 }
               }
             }
@@ -573,7 +722,7 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
       RAVE_FREE(v);
       RAVE_FREE(z);
       RAVE_FREE(az);
-      RAVE_INFO0("Could not find any acceptable scans, dropping out");
+      RAVE_INFO0("Could not find any acceptable scans, dropping out...");
       goto done;
     } 
 
@@ -612,35 +761,36 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
         vdir = vdir - 360;
       }
 
-      /* standard deviation of the wind velocity */
+      /* RMSE of the wind velocity */
       for (i=0; i<nv; i++) {
         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin(*(az+i)+beta)),2);
       }
-      vstd = sqrt (vstd/(nv-1));
+      vstd = sqrt (vstd/nv);
       
       /* Calculate the x-component (East) and y-component (North) of the wind 
          velocity using the wind direction and the magnitude of the wind velocity */
       vdir_rad = vdir * DEG2RAD;
-      u_wnd_comp = vvel * cos(vdir_rad);
-      v_wnd_comp = vvel * sin(vdir_rad);
+      u_wnd_comp = vvel * sin(vdir_rad - M_PI);
+      v_wnd_comp = vvel * cos(vdir_rad - M_PI);
     }
 
     // reflectivity calculations
-    if (nz > 1) {
-      /* standard deviation of reflectivity */
+    if (nz > 0) {
+      /* RMSE of the reflectivity */
       for (i = 0; i < nz; i++) {
         zstd = zstd + pow(*(z+i) - (zsum/nz),2);
       }
       zmean = Z2dBZ(zsum/nz);
-      zstd = sqrt(zstd/(nz-1));
+      zstd = sqrt(zstd/nz);
       zstd = Z2dBZ(zstd);
     }
 
     /* Set the hght_field values */
     if (hght_field != NULL) RaveField_setValue(hght_field, 0, yindex, centerOfLayer / 1000.0); /* in km */
 
-    /* If the number of points for wind is larger tahn threshold, set values, else set nodata */
-    if (nv < NMIN) {
+    /* If the number of points for wind is smaller than the threshold nmin_wnd or the calculated wind velocity is larger than */
+    /* threshold ff_max, set nodata, otherwise set values. */
+    if ((nv < self->nmin_wnd) || (vvel > self->ff_max)) {
       if (nv_field != NULL) RaveField_setValue(nv_field, 0, yindex, -1.0); /* nodata for counter */
       if (uwnd_field != NULL) RaveField_setValue(uwnd_field, 0, yindex, self->nodata_VP);
       if (vwnd_field != NULL) RaveField_setValue(vwnd_field, 0, yindex, self->nodata_VP);
@@ -657,7 +807,7 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
     }
 
     /* If the number of points for reflectivity is larger than threshold, set values, else set nodata */
-    if (nz < NMIN) {
+    if (nz < self->nmin_ref) {
       if (nz_field != NULL) RaveField_setValue(nz_field, 0, yindex, -1.0);
       if (dbzh_field != NULL) RaveField_setValue(dbzh_field, 0, yindex, self->nodata_VP);
       if (dbzh_dev_field != NULL) RaveField_setValue(dbzh_dev_field, 0, yindex, self->nodata_VP);
@@ -692,14 +842,15 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
     VerticalProfile_setLevels(result, ysize);
 
     /* Below are the fields that are possible to inject into the profile, note */
-    /* that the nz_field is not present, this requires that the ODIM spec */
-    /* is adjusted to allow TWO sample size arrays e.g. nv and nz, current spec */
+    /* that the presence of the nz_field violates the ODIM spec which defines only one sample size. */
+    /* We allow anyway TWO sample size arrays i.e. nv (named n in the output profile) and nz, current spec */
     /* contains only one sample size array n, with obvious consequences IF a */
     /* combined (wind + refl), or a pure refl, profile is created */
 
     if ((uwnd_field != NULL && !VerticalProfile_setUWND(result, uwnd_field)) ||
         (vwnd_field != NULL && !VerticalProfile_setVWND(result, vwnd_field)) ||
         (nv_field != NULL && !VerticalProfile_setNV(result, nv_field)) ||
+        (nz_field != NULL && !VerticalProfile_setNZ(result, nz_field)) ||
         (hght_field != NULL && !VerticalProfile_setHGHT(result, hght_field)) ||
         (ff_field != NULL && !VerticalProfile_setFF(result, ff_field)) ||
         (ff_dev_field != NULL && !VerticalProfile_setFFDev(result, ff_dev_field)) ||
@@ -728,28 +879,50 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char*
   VerticalProfile_setEndTime(result, RaveDateTime_getTime(lastEndDT));
   VerticalProfile_setProduct(result, product);
    
-  /* Need to filter the attribute data with respect to selected min elangle
-     Below attributes satisfy reguirements from E-profile, add other when needed */
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/highprf", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/lowprf", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/pulsewidth", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/wavelength", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/RXbandwidth", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/RXlossH", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/TXlossH", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/antgainH", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/azmethod", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/binmethod", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/malfunc", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/nomTXpower", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/radar_msg", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/radconstH", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/radomelossH", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/rpm", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/software", self->emin);
-  WrwpInternal_findAndAddAttribute(result, inobj, "how/system", self->emin);
+  /* Supported but not included how attributes, add when needed*/
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/highprf", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/lowprf", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/pulsewidth", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/wavelength", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/RXbandwidth", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/RXlossH", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/TXlossH", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/antgainH", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/azmethod", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/binmethod", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/malfunc", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/nomTXpower", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/radar_msg", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/radconstH", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/radomelossH", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/rpm", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/software", self->emin, self->emax);
+  //WrwpInternal_findAndAddAttribute(result, inobj, "how/system", self->emin, self->emax);
+  
+  /*Dealing with the array of strings containing uniques how/task attribute
+    This is converted to a string that will represent the unique tasks for scans in the volume
+    Note that for radar data not having this attribute in their volumes, it will not be written 
+    Initialization of some useful variables */
+  ntask = len(taskArgs);
+  char finalTasks[300] = {'\0'};
+  int taskCount = 0;
+
+  for (i = 0; i < ntask; i++) {
+    if (taskArgs[i]) {
+      taskCount = taskCount + 1;
+      strcat(finalTasks, taskArgs[i]);
+      strcat(finalTasks, comma);
+    }
+  }
+  if (taskCount != 0) {
+    WrwpInternal_LastcharDel(finalTasks); 
+    WrwpInternal_addStringAttribute(result, "how/task", finalTasks);
+  }
+  /* how attributes requested by Eprofile */
+  WrwpInternal_addStringAttribute(result, "how/angles", theUsedElevationAngles);
   WrwpInternal_addDoubleAttribute(result, "how/minrange", (double)Wrwp_getDMIN(self) / 1000.0); /* km */
   WrwpInternal_addDoubleAttribute(result, "how/maxrange", (double)Wrwp_getDMAX(self) / 1000.0); /* km */
+  RAVE_FREE(theUsedElevationAngles);
 
 done:
   RAVE_OBJECT_RELEASE(polnav);
index a3f9cbc..ccb4e17 100644 (file)
@@ -56,22 +56,11 @@ along with baltrad-wrwp.  If not, see <http://www.gnu.org/licenses/>.
 
 #define DEG2RAD     DEG_TO_RAD      /* Degrees to radians. From PROJ.4 */
 #define RAD2DEG     RAD_TO_DEG      /* Radians to degrees. From PROJ.4 */
-#define NOR         20000           /* Number of rows in matrix A used in the computation */
+#define NOR         40000           /* Number of rows in matrix A used in the computation */
 #define NOC         3               /* Number of columns in matrix A used in the computation */
 #define NRHS        1               /* Number of right-hand sides; that is, the number of columns in matrix B used in the computation */
 #define LDA         NOC             /* Leading dimension of the array specified for a */
 #define LDB         NRHS            /* Leading dimension of the array specified for b */
-#define DMIN        4000            /* Minimum distance for deriving a profile [m] */
-#define DMAX        40000           /* Maximum distance for deriving a profile [m] */
-#define NMIN        36              /* Minimum sample size */
-#define EMIN        2.5             /* Minimum elevation angle [deg] */
-#define VMIN        2.0             /* Radial velocity threshold [m/s] */
-#define DZ          200             /* Height interval for deriving a profile [m] */
-#define HMAX        12000           /* Maximum height of the profile [m] */
-#define NODATA_VP   -9999           /* Nodata value used in the vertical profile */
-#define UNDETECT_VP -9999           /* Undetect value used in the vertical profile */         
-#define GAIN_VP     1.0             /* Gain value for the fields UWND and VWND */
-#define OFFSET_VP   0.0             /* Offset value for the fields UWND and VWND */
 
 /**
  * Defines a weather radar wind product generator
@@ -84,6 +73,13 @@ typedef struct _Wrwp_t Wrwp_t;
 extern RaveCoreObjectType Wrwp_TYPE;
 
 /**
+ * Returns the height interval for deriving a profile [m]
+ * @param[in] self - self
+ * @return the height interval (default value is DZ)
+ */
+int Wrwp_getDZ(Wrwp_t* self);
+
+/**
  * Sets the height interval for deriving a profile [m]
  * @param[in] self - self
  * @param[in] dz - the height interval
@@ -147,11 +143,18 @@ double Wrwp_getOFFSET_VP(Wrwp_t* self);
 void Wrwp_setOFFSET_VP(Wrwp_t* self, double offset_vp);
 
 /**
- * Returns the height interval for deriving a profile [m]
+ * Sets maximum allowed velocity for each layer in the profile [m/s]
  * @param[in] self - self
- * @return the height interval (default value is DZ)
+ * @param[in] ff_max - maximum allowed velocity for each layer in the profile
  */
-int Wrwp_getDZ(Wrwp_t* self);
+void Wrwp_setFF_MAX(Wrwp_t* self, double ff_max);
+
+/**
+ * Returns maximum allowed velocity for each layer in the profile [m/s]
+ * @param[in] self - self
+ * @return maximum allowed velocity for each layer in the profile [m/s] (default is FF_MAX)
+ */
+double Wrwp_getFF_MAX(Wrwp_t* self);
 
 /**
  * Sets maximum height of the profile [m]
@@ -168,6 +171,34 @@ void Wrwp_setHMAX(Wrwp_t* self, int hmax);
 int Wrwp_getHMAX(Wrwp_t* self);
 
 /**
+ * Sets minimum number of points for wind
+ * @param[in] self - self
+ * @param[in] nmin_wnd - minimum number of points for wind
+ */
+void Wrwp_setNMIN_WND(Wrwp_t* self, int nmin_wnd);
+
+/**
+ * Returns minimum number of points for wind
+ * @param[in] self - self
+ * @return minimum number of points for wind (default NMIN_WND)
+ */
+int Wrwp_getNMIN_WND(Wrwp_t* self);
+
+/**
+ * Sets minimum number of points for reflectivity
+ * @param[in] self - self
+ * @param[in] nmin_ref - minimum number of points for reflectivity
+ */
+void Wrwp_setNMIN_REF(Wrwp_t* self, int nmin_ref);
+
+/**
+ * Returns minimum number of points for reflectivity
+ * @param[in] self - self
+ * @return minimum number of points for reflectivity (default NMIN_REF)
+ */
+int Wrwp_getNMIN_REF(Wrwp_t* self);
+
+/**
  * Sets minimum distance for deriving a profile [m]
  * @param[in] self - self
  * @param[in] dmin - minimum distance for deriving a profile [m]
@@ -210,6 +241,20 @@ void Wrwp_setEMIN(Wrwp_t* self, double emin);
 double Wrwp_getEMIN(Wrwp_t* self);
 
 /**
+ * Sets maximum elevation angle [deg]
+ * @param[in] self - self
+ * @param[in] emax - maximum elevation angle [deg]
+ */
+void Wrwp_setEMAX(Wrwp_t* self, double emax);
+
+/**
+ * Returns maximum elevation angle [deg]
+ * @param[in] self - self
+ * @return Maximum elevation angle [deg] (default EMAX)
+ */
+double Wrwp_getEMAX(Wrwp_t* self);
+
+/**
  * Sets radial velocity threshold [m/s]
  * @param[in] self - self
  * @param[in] vmin - radial velocity threshold [m/s]
diff --git a/pywrwp/VPodimVersionConverter.py b/pywrwp/VPodimVersionConverter.py
new file mode 100644 (file)
index 0000000..a24cc6b
--- /dev/null
@@ -0,0 +1,231 @@
+'''
+Copyright (C) 2015- Swedish Meteorological and Hydrological Institute (SMHI)
+
+Converts WRWP files in ODIM-H5 ver2.2 to WRWP files in ODIM-H5 ver 2.1 files.
+
+@file
+@author Ulf E. Nordh, SMHI
+@date 2019-10-23'''
+
+import _pyhl
+import numpy
+import re, sys, traceback, time
+
+DATATYPES = {"char":"char",
+             "int8":"char",
+             "uint8":"uchar",
+             "int16":"short",
+             "uint16":"ushort",
+             "int32":"int",
+             "int64":"long",
+             "float32":"float",
+             "float64":"double"}
+
+# Main class for the version converter.
+class VPodimVersionConverter(object):
+  def __init__(self, filename, quantities):
+    self._filename = filename
+    if not _pyhl.is_file_hdf5(filename):
+      raise Exception("Not a HDF5 file")
+    self._nodelist = _pyhl.read_nodelist(self._filename)
+    self._nodelist.selectAll()
+    self._nodelist.fetch()    
+    self._nodenames = self._nodelist.getNodeNames().keys()
+    self._converted_files = []  # Contains tuples of (nodelist, suggested name)
+
+  def convert(self,quantities, filename):
+    if self.isVerticalProfile():
+      self._convertVP(quantities, filename)
+    else:
+      raise Exception("No Support for file type: %s"%self.getWhatObject())
+    return self._nodelist, self._converted_files
+  
+  def isVerticalProfile(self):
+    return self.getWhatObject() in ["VP"]
+    
+  def isSupported(self):
+    return self.isVerticalProfile()
+    
+  def getWhatObject(self):
+    return self._nodelist.getNode("/what/object").data()
+  
+  def _getDatasetAndDataIndexForQuantity(self, quantity):
+    dsetidx = 1
+    dsetloop = True
+    while dsetloop:
+      dloop = True
+      dsetname = "/dataset%i"%dsetidx
+      if dsetname not in self._nodenames:
+        dsetloop = False
+        break
+      didx = 1
+      while dloop:
+        try:
+          name = "/dataset%i/data%i"%(dsetidx,didx)
+          if name in self._nodenames:
+            dataname = "/dataset%i/data%i/what/quantity"%(dsetidx,didx)
+            if dataname in self._nodenames:
+              if self._nodelist.getNode(dataname).data() == quantity:
+                return (dsetidx, didx)
+          else:
+            dloop = False
+        except Exception:
+          traceback.print_exc(file=sys.stdout)
+          dloop = False
+        didx = didx + 1
+      dsetidx = dsetidx + 1
+    return None
+    
+  def _getDatasetByQuantity(self, quantity):
+    indices = self._getDatasetAndDataIndexForQuantity(quantity)
+    if indices != None:
+      name = "/dataset%i/data%i/data"%(indices[0], indices[1])
+      return name
+    return None
+    
+  def _convertValueTypeToStr(self, value):
+    if type(value) is str:
+      result = "string"
+    elif type(value) is float:
+      result = "float"
+    elif type(value) is int:
+      result = "int"
+    else:
+      raise Exception("Unsupported data type")
+    return result
+
+  def _copyData(self, nodelist, name, oname=None):
+    d = self._nodelist.getNode(name).data()
+    datatype = str(d.dtype)
+    if datatype in DATATYPES:
+      translDatatype= DATATYPES[datatype]
+    else:
+      raise Exception("Unsupported datatype %s"%datatype)
+
+    c = _pyhl.compression(_pyhl.COMPRESSION_ZLIB)
+    c.level = 6
+    if oname:
+      n = _pyhl.node(_pyhl.DATASET_ID, oname, c)
+    else:
+      n = _pyhl.node(_pyhl.DATASET_ID, name, c)
+    n.setArrayValue(-1, list(d.shape), d, translDatatype, -1)
+    nodelist.addNode(n)
+
+  def _addGroup(self, nodelist, name):
+    node = _pyhl.node(_pyhl.GROUP_ID, name)
+    nodelist.addNode(node)
+
+  def _addAttribute(self, nodelist, name, value):
+    node = _pyhl.node(_pyhl.ATTRIBUTE_ID, name)
+    node.setScalarValue(-1, value, self._convertValueTypeToStr(value), -1)
+    nodelist.addNode(node)
+
+  def _copyAttribute(self, nodelist, name, ntype=None, oname=None):
+    d = self._nodelist.getNode(name).data()
+    if oname:
+      n = _pyhl.node(_pyhl.ATTRIBUTE_ID, oname)
+    else:
+      n = _pyhl.node(_pyhl.ATTRIBUTE_ID, name)
+      
+    if ntype:
+      n.setScalarValue(-1, d, ntype, -1)
+    else:
+      n.setScalarValue(-1, d, self._convertValueTypeToStr(d), -1)
+    
+    nodelist.addNode(n)
+
+  def _copyAttributeIfExistsToGroup(self, nodelist, name, toname=None):
+    if name in self._nodenames:
+      if toname == None:
+        toname = name
+      self._copyAttribute(nodelist, name, None, toname)
+      return True
+    return False
+
+  def _copyWhatObject(self, nodelist):
+    d = self._nodelist.getNode("/what/object").data()
+    self._addAttribute(nodelist, "/what/object", d)
+    
+  def _populateNodelistWithDataAndAttributes(self, nodelist, quantity):
+
+    if quantity == "NV":
+      quantity = "n"
+    if quantity == "NZ":
+      quantity = "nz"
+
+    datasetLocation = self._getDatasetByQuantity(quantity)
+
+    # Add the datsets and their attributes to the nodelist
+    data = self._nodelist.getNode(datasetLocation).data()
+    self._copyData(nodelist, datasetLocation, datasetLocation)
+    self._copyAttributeIfExistsToGroup(nodelist, datasetLocation[:-4] + "what/gain", datasetLocation[:-4] + "what/gain")
+    self._copyAttributeIfExistsToGroup(nodelist, datasetLocation[:-4] + "what/offset", datasetLocation[:-4] + "what/offset")
+    self._copyAttributeIfExistsToGroup(nodelist, datasetLocation[:-4] + "what/nodata", datasetLocation[:-4] + "what/nodata")
+
+    if quantity == "DBZH":
+      self._addAttribute(nodelist, datasetLocation[:-4] + "what/quantity", "dbz")
+    elif quantity == "DBZH_dev":
+      self._addAttribute(nodelist, datasetLocation[:-4] + "what/quantity", "dbz_dev")
+    else:
+      self._copyAttributeIfExistsToGroup(nodelist, datasetLocation[:-4] + "what/quantity", datasetLocation[:-4] + "what/quantity")
+
+    self._copyAttributeIfExistsToGroup(nodelist, datasetLocation[:-4] + "what/undetect", datasetLocation[:-4] + "what/undetect")   
+    
+  def _addVPInformation(self, nodelist, quantities):
+    # Prepare the nodelist
+    self._addGroup(nodelist, "/how")
+    self._addGroup(nodelist, "/what")
+    self._addGroup(nodelist, "/where")
+    self._addGroup(nodelist, "/dataset1")
+    self._addGroup(nodelist, "/dataset1/what")
+
+    for i in range(1, len(quantities) + 1):
+
+      self._addGroup(nodelist, "/dataset1/data%d"%i)
+      self._addGroup(nodelist, "/dataset1/data%d/what"%i)
+    
+    # Root attribute
+    self._addAttribute(nodelist, "/Conventions", "ODIM_H5/V2_1")
+
+    # Attributes under /what
+    self._copyAttributeIfExistsToGroup(nodelist, "/what/date", "/what/date")
+    self._copyAttributeIfExistsToGroup(nodelist, "/what/object", "/what/object")
+    self._copyAttributeIfExistsToGroup(nodelist, "/what/time", "/what/time")
+    self._addAttribute(nodelist, "/what/version", "H5rad 2.1")
+    self._copyAttributeIfExistsToGroup(nodelist, "/what/source", "/what/source")
+
+    # Attributes under /where
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/height", "/where/height")
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/interval", "/where/interval")
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/lat", "/where/lat")
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/levels", "/where/levels")
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/lon", "/where/lon")
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/maxheight", "/where/maxheight")
+    self._copyAttributeIfExistsToGroup(nodelist, "/where/minheight", "/where/minheight")
+
+    # Attribute under /dataset1/what
+    self._copyAttributeIfExistsToGroup(nodelist, "/dataset1/what/product", "/dataset1/what/product")
+    self._copyAttributeIfExistsToGroup(nodelist, "/dataset1/what/enddate", "/dataset1/what/enddate")
+    self._copyAttributeIfExistsToGroup(nodelist, "/dataset1/what/endtime", "/dataset1/what/endtime")
+    self._copyAttributeIfExistsToGroup(nodelist, "/dataset1/what/startdate", "/dataset1/what/startdate")
+    self._copyAttributeIfExistsToGroup(nodelist, "/dataset1/what/starttime", "/dataset1/what/starttime")
+        
+    # Attributes under /how
+    self._copyAttributeIfExistsToGroup(nodelist, "/how/angles", "/how/angles")
+    self._copyAttributeIfExistsToGroup(nodelist, "/how/maxrange", "/how/maxrange")
+    self._copyAttributeIfExistsToGroup(nodelist, "/how/minrange", "/how/minrange")
+    self._copyAttributeIfExistsToGroup(nodelist, "/how/task", "/how/task")
+    
+    for QUANTITY in quantities:
+      self._populateNodelistWithDataAndAttributes(nodelist, QUANTITY)
+
+  def _convertVP(self, quantities, filename):
+    nodelist = _pyhl.nodelist()
+    quantities = [str(x) for x in quantities.split(",")]
+
+    # Populate the nodelist with attributes and datasets needed for DIANA
+    self._addVPInformation(nodelist, quantities)
+
+    # Create the filename and send the modified file back to b2d
+    self._converted_files.append((nodelist, filename))
+
index 77c5c88..bd49fda 100644 (file)
@@ -41,6 +41,13 @@ along with baltrad-wrwp.  If not, see <http://www.gnu.org/licenses/>.
 ## Slight adjustment done to the call to wrwp.generate in function generate.
 ## The call is done with try-except to deal with the situation that
 ## the returned result is NULL. In addition, logging is inserted.
+##
+## @co-autor Ulf E. Nordh, SMHI
+## @date 2019-10-08
+##
+## From additional requirements, some of the parameters given in wrwp.h is given in a configuration
+## file instead for greater ease of adjustments without having to re-compile the code. The selected format
+## is xml.
 
 import _wrwp
 import _rave
@@ -52,11 +59,15 @@ import rave_pgf_logger
 import rave_tempfile
 import odim_source
 import math
-
+import os
+import sys
+import fnmatch
+import xml.etree.cElementTree as ET    
 from rave_defines import CENTER_ID, GAIN, OFFSET
 
 logger = rave_pgf_logger.create_logger()
 
+
 ravebdb = None
 try:
   import rave_bdb
@@ -64,6 +75,15 @@ try:
 except:
   pass
 
+# Locates a file (pattern) in a directory tree (path)
+def find(pattern, path):
+  result = []
+  for root, dirs, files in os.walk(path):
+    for name in files:
+      if fnmatch.fnmatch(name, pattern):
+        result.append(os.path.join(root, name))
+  return result
+
 ## Creates a dictionary from a rave argument list
 #@param arglist the argument list
 #@return a dictionary
@@ -98,6 +118,8 @@ def generate(files, arguments):
   args = arglist2dict(arguments)
   wrwp = _wrwp.new()
   fields = None
+
+  # Parameter values from the web-GUI
   if "interval" in args.keys():
     wrwp.dz = strToNumber(args["interval"])
   if "maxheight" in args.keys():
@@ -113,6 +135,31 @@ def generate(files, arguments):
   if "fields" in args.keys():
     fields = args["fields"]
 
+  # Finds the wrwp_config.xml under /local_disk/baltrad/blt_sys/ and sets the auxillary parameters not defined in the web-GUI.
+  # Usually the installation paths in the node-installer is not changed relative to default, the use of find locates the wrwp_config.xml
+  # if another path is set in the node-installer, the only requirement is that the wrwp_config.xml is somewhere under /local_disk/baltrad/blt_sys.
+
+  path2config = find('wrwp_config.xml', '/local_disk/baltrad/blt_sys')
+
+  root = ET.parse(str(path2config[0])).getroot()
+  for param in root.findall('param'):
+    if param.get('name') == 'EMAX':
+      wrwp.emax = strToNumber(param.find('value').text)
+    if param.get('name') == 'NMIN_WND':
+      wrwp.nmin_wnd = strToNumber(param.find('value').text)
+    if param.get('name') == 'NMIN_REF':
+      wrwp.nmin_ref = strToNumber(param.find('value').text)
+    if param.get('name') == 'FF-MAX':
+      wrwp.ff_max = strToNumber(param.find('value').text)
+    if param.get('name') == 'NODATA_VP':
+      wrwp.nodata_VP = strToNumber(param.find('value').text)
+    if param.get('name') == 'UNDETECT_VP':
+      wrwp.undetect_VP = strToNumber(param.find('value').text)
+    if param.get('name') == 'GAIN_VP':
+      wrwp.gain_VP = strToNumber(param.find('value').text)
+    if param.get('name') == 'OFFSET_VP':
+      wrwp.offset_VP = strToNumber(param.find('value').text)
+
   if len(files) != 1:
     raise AttributeError("Must call plugin with _one_ polar volume")
   
index d6bf7cd..be81c1b 100644 (file)
@@ -203,8 +203,24 @@ static PyObject* _pywrwp_getattro(PyWrwp* self, PyObject* name)
     return PyInt_FromLong(Wrwp_getDMAX(self->wrwp));
   } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("emin", name) == 0) {
     return PyFloat_FromDouble(Wrwp_getEMIN(self->wrwp));
+  }  else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("emax", name) == 0) {
+    return PyFloat_FromDouble(Wrwp_getEMAX(self->wrwp));
   } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("vmin", name) == 0) {
     return PyFloat_FromDouble(Wrwp_getVMIN(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("nmin_wnd", name) == 0) {
+    return PyInt_FromLong(Wrwp_getNMIN_WND(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("nmin_ref", name) == 0) {
+    return PyInt_FromLong(Wrwp_getNMIN_REF(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("ff_max", name) == 0) {
+    return PyFloat_FromDouble(Wrwp_getFF_MAX(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("nodata_VP", name) == 0) {
+    return PyInt_FromLong(Wrwp_getNODATA_VP(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("undetect_VP", name) == 0) {
+    return PyInt_FromLong(Wrwp_getUNDETECT_VP(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("gain_VP", name) == 0) {
+    return PyFloat_FromDouble(Wrwp_getGAIN_VP(self->wrwp));
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("offset_VP", name) == 0) {
+    return PyFloat_FromDouble(Wrwp_getOFFSET_VP(self->wrwp));
   }
   return PyObject_GenericGetAttr((PyObject*)self, name);
 }
@@ -250,6 +266,14 @@ static int _pywrwp_setattro(PyWrwp* self, PyObject* name, PyObject* val)
     } else {
       raiseException_gotoTag(done, PyExc_TypeError, "emin must be an integer or a float");
     }
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("emax", name) == 0) {
+    if (PyFloat_Check(val)) {
+      Wrwp_setEMAX(self->wrwp, PyFloat_AsDouble(val));
+    } else if (PyInt_Check(val)) {
+      Wrwp_setEMAX(self->wrwp, (double)PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "emax must be an integer or a float");
+    }
   } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("vmin", name) == 0) {
     if (PyFloat_Check(val)) {
       Wrwp_setVMIN(self->wrwp, PyFloat_AsDouble(val));
@@ -258,7 +282,55 @@ static int _pywrwp_setattro(PyWrwp* self, PyObject* name, PyObject* val)
     } else {
       raiseException_gotoTag(done, PyExc_TypeError, "vmin must be an integer or a float");
     }
-  }
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("nmin_wnd", name) == 0) {
+    if (PyInt_Check(val)) {
+      Wrwp_setNMIN_WND(self->wrwp, PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "nmin_wnd must be an integer");
+    }
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("nmin_ref", name) == 0) {
+    if (PyInt_Check(val)) {
+      Wrwp_setNMIN_REF(self->wrwp, PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "nmin_ref must be an integer");
+    }
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("ff_max", name) == 0) {
+    if (PyFloat_Check(val)) {
+      Wrwp_setFF_MAX(self->wrwp, PyFloat_AsDouble(val));
+    } else if (PyInt_Check(val)) {
+      Wrwp_setFF_MAX(self->wrwp, (double)PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "ff_max must be an integer or a float");
+    }
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("nodata_VP", name) == 0) {
+    if (PyInt_Check(val)) {
+      Wrwp_setNODATA_VP(self->wrwp, PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "nodata_VP must be an integer");
+    }
+  } else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("undetect_VP", name) == 0) {
+    if (PyInt_Check(val)) {
+      Wrwp_setUNDETECT_VP(self->wrwp, PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "undetect_VP must be an integer");
+    }
+  }  else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("gain_VP", name) == 0) {
+    if (PyFloat_Check(val)) {
+      Wrwp_setGAIN_VP(self->wrwp, PyFloat_AsDouble(val));
+    } else if (PyInt_Check(val)) {
+      Wrwp_setGAIN_VP(self->wrwp, (double)PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "gain_VP must be an integer or a float");
+    }
+  }  else if (PY_COMPARE_STRING_WITH_ATTRO_NAME("offset_VP", name) == 0) {
+    if (PyFloat_Check(val)) {
+      Wrwp_setOFFSET_VP(self->wrwp, PyFloat_AsDouble(val));
+    } else if (PyInt_Check(val)) {
+      Wrwp_setOFFSET_VP(self->wrwp, (double)PyInt_AsLong(val));
+    } else {
+      raiseException_gotoTag(done, PyExc_TypeError, "offset_VP must be an integer or a float");
+    }
+  } 
 
   result = 0;
 done:
index 8c43ae9..681510c 100644 (file)
@@ -29,11 +29,84 @@ import string
 import _wrwp
 import _helpers
 import _raveio, _rave
+import xml.etree.cElementTree as ET
+import sys
+import os
+
+sys.path.append(os.path.realpath(__file__))
+
+def strToNumber(strval):
+  # Converts a string into a number, either int or float
+  # strval: the string to translate
+  # return:the translated value
+  # throws ValueError if value not could be translated
+  if type(strval) is not str: # Avoid doing anything if it is not a string as input
+    return strval
+  else:
+    try:
+      return int(strval)
+    except ValueError:
+      return float(strval)
+
+# Parse the xml-file to get the default values, the file is placed in the directory below i.e. fixtures
+root = ET.parse('./fixtures/wrwp_config.xml').getroot()
+for param in root.findall('param'):
+  if param.get('name') == 'DMIN':
+    DMIN = strToNumber(param.find('value').text)
+  if param.get('name') == 'DMAX':
+    DMAX = strToNumber(param.find('value').text)
+  if param.get('name') == 'NMIN_WND':
+    NMIN_WND = strToNumber(param.find('value').text)
+  if param.get('name') == 'NMIN_REF':
+    NMIN_REF = strToNumber(param.find('value').text)
+  if param.get('name') == 'EMIN':
+    EMIN = strToNumber(param.find('value').text)
+  if param.get('name') == 'EMAX':
+    EMAX = strToNumber(param.find('value').text)
+  if param.get('name') == 'VMIN':
+    VMIN = strToNumber(param.find('value').text)
+  if param.get('name') == 'FF_MAX':
+    FF_MAX = strToNumber(param.find('value').text)
+  if param.get('name') == 'DZ':
+    DZ = strToNumber(param.find('value').text)
+  if param.get('name') == 'HMAX':
+    HMAX = strToNumber(param.find('value').text)
+  if param.get('name') == 'NODATA_VP':
+    NODATA_VP = strToNumber(param.find('value').text)
+  if param.get('name') == 'UNDETECT_VP':
+    UNDETECT_VP = strToNumber(param.find('value').text)
+  if param.get('name') == 'GAIN_VP':
+    GAIN_VP = strToNumber(param.find('value').text)
+  if param.get('name') == 'OFFSET_VP':
+    OFFSET_VP = strToNumber(param.find('value').text)
+  if param.get('name') == 'QUANTITIES':
+    QUANTITIES = param.find('value').text
+
+
+def load_wrwp_defaults_to_obj():
+  wrwp = _wrwp.new()
+  wrwp.hmax = HMAX
+  wrwp.dz = DZ
+  wrwp.emin = EMIN
+  wrwp.dmin = DMIN
+  wrwp.dmax = DMAX
+  wrwp.nmin_wnd = NMIN_WND
+  wrwp.nmin_ref = NMIN_REF
+  wrwp.emax = EMAX
+  wrwp.vmin = VMIN
+  wrwp.ff_max = FF_MAX
+  wrwp.nodata_VP = NODATA_VP
+  wrwp.undetect_VP = UNDETECT_VP
+  wrwp.gain_VP = GAIN_VP
+  wrwp.offset_VP = OFFSET_VP
+
+  return wrwp
 
 class WrwpTest(unittest.TestCase):
   FIXTURE = "fixtures/pvol_seang_20090501T120000Z.h5"
   FIXTURE2 = "fixtures/selul_pvol_20151114T1615Z.h5"
   FIXTURE3 = "fixtures/seang_zdr_qcvol_one_scan_and_elangle_lower_than_wrwp_threshold.h5"
+  FIXTURE4 = "fixtures/seosd_qcvol_zdrvol_different_task.h5"
   
   def setUp(self):
     _helpers.triggerMemoryStatus()
@@ -46,7 +119,7 @@ class WrwpTest(unittest.TestCase):
     self.assertNotEqual(-1, str(type(obj)).find("WrwpCore"))
 
   def test_dz(self):
-    obj = _wrwp.new()
+    obj = load_wrwp_defaults_to_obj()
     self.assertEqual(200, obj.dz)
     obj.dz = 100
     self.assertEqual(100, obj.dz)
@@ -58,7 +131,7 @@ class WrwpTest(unittest.TestCase):
     self.assertEqual(100, obj.dz)
 
   def test_hmax(self):
-    obj = _wrwp.new()
+    obj = load_wrwp_defaults_to_obj()
     self.assertEqual(12000, obj.hmax)
     obj.hmax = 100
     self.assertEqual(100, obj.hmax)
@@ -70,8 +143,8 @@ class WrwpTest(unittest.TestCase):
     self.assertEqual(100, obj.hmax)
 
   def test_dmin(self):
-    obj = _wrwp.new()
-    self.assertEqual(4000, obj.dmin)
+    obj = load_wrwp_defaults_to_obj()
+    self.assertEqual(5000, obj.dmin)
     obj.dmin = 100
     self.assertEqual(100, obj.dmin)
     try:
@@ -82,8 +155,8 @@ class WrwpTest(unittest.TestCase):
     self.assertEqual(100, obj.dmin)
     
   def test_dmax(self):
-    obj = _wrwp.new()
-    self.assertEqual(40000, obj.dmax)
+    obj = load_wrwp_defaults_to_obj()
+    self.assertEqual(25000, obj.dmax)
     obj.dmax = 100
     self.assertEqual(100, obj.dmax)
     try:
@@ -94,28 +167,62 @@ class WrwpTest(unittest.TestCase):
     self.assertEqual(100, obj.dmax)
 
   def test_emin(self):
-    obj = _wrwp.new()
-    self.assertAlmostEqual(2.5, obj.emin, 4)
+    obj = load_wrwp_defaults_to_obj()
+    self.assertAlmostEqual(0.5, obj.emin, 4)
     obj.emin = 3.5
     self.assertAlmostEqual(3.5, obj.emin, 4)
     obj.emin = 4
     self.assertAlmostEqual(4.0, obj.emin, 4)
 
+  def test_emax(self):
+    obj = load_wrwp_defaults_to_obj()
+    self.assertAlmostEqual(45.0, obj.emax, 4)
+    obj.emax = 35.0
+    self.assertAlmostEqual(35.0, obj.emax, 4)
+    obj.emax = 4
+    self.assertAlmostEqual(4.0, obj.emax, 4)
+
   def test_vmin(self):
-    obj = _wrwp.new()
+    obj = load_wrwp_defaults_to_obj()
     self.assertAlmostEqual(2.0, obj.vmin, 4)
     obj.vmin = 3.5
     self.assertAlmostEqual(3.5, obj.vmin, 4)
     obj.vmin = 4
     self.assertAlmostEqual(4.0, obj.vmin, 4)
-  
+
+  def test_ff_max(self):
+    obj = load_wrwp_defaults_to_obj()
+    self.assertAlmostEqual(60.0, obj.ff_max, 4)
+    obj.ff_max = 3.5
+    self.assertAlmostEqual(3.5, obj.ff_max, 4)
+    obj.ff_max = 90.0
+    self.assertAlmostEqual(90.0, obj.ff_max, 4)
+
+  def test_nmin_wnd(self):
+    obj = load_wrwp_defaults_to_obj()
+    self.assertEqual(40, obj.nmin_wnd, 4)
+    obj.nmin_wnd = 70
+    self.assertEqual(70, obj.nmin_wnd, 4)
+    obj.nmin_wnd = 20
+    self.assertEqual(20, obj.nmin_wnd, 4)
+
+  def test_nmin_ref(self):
+    obj = load_wrwp_defaults_to_obj()
+    self.assertEqual(40, obj.nmin_ref, 4)
+    obj.nmin_ref = 70
+    self.assertEqual(70, obj.nmin_ref, 4)
+    obj.nmin_ref = 20
+    self.assertEqual(20, obj.nmin_ref, 4)
+
   def test_generate(self):
     pvol = _raveio.open(self.FIXTURE).object
-    generator = _wrwp.new()
-    generator.hmax = 2000
-    generator.dz = 200
+    wrwp = load_wrwp_defaults_to_obj()
+    wrwp.hmax = 2000
+    wrwp.dz = 200
+    fields = None
     
-    vp = generator.generate(pvol, "UWND,VWND,HGHT,NV")
+    fields = QUANTITIES
+    vp = wrwp.generate(pvol, fields)
     
     uwnd = vp.getUWND()
     vwnd = vp.getVWND()
@@ -145,26 +252,83 @@ class WrwpTest(unittest.TestCase):
     self.assertEqual(pvol.source, vp.source)
     self.assertEqual(pvol.date, vp.date)
     self.assertEqual(pvol.time, vp.time)
+
+  def test_generate_from_default(self):
+    pvol = _raveio.open(self.FIXTURE).object
+    wrwp = load_wrwp_defaults_to_obj()
+    fields = None
+
+    exceptionTest = False
+
+    vp = wrwp.generate(pvol, fields)
+
+    ff = vp.getFF()
+    dd = vp.getDD()
+    dbzh = vp.getDBZ()
+    dbzh_dev = vp.getDBZDev()
+    NZ = vp.getNZ()
+    uwnd = vp.getUWND()
+
+    self.assertEqual(1, ff.xsize)
+    self.assertEqual(60, ff.ysize)
+    self.assertEqual("ff", ff.getAttribute("what/quantity"))
+
+    self.assertEqual(1, dd.xsize)
+    self.assertEqual(60, dd.ysize)
+    self.assertEqual("dd", dd.getAttribute("what/quantity"))
+
+    self.assertEqual(1, dbzh.xsize)
+    self.assertEqual(60, dbzh.ysize)
+    self.assertEqual("DBZH", dbzh.getAttribute("what/quantity"))
+
+    self.assertEqual(1, dbzh_dev.xsize)
+    self.assertEqual(60, dbzh_dev.ysize)
+    self.assertEqual("DBZH_dev", dbzh_dev.getAttribute("what/quantity"))
+
+    self.assertEqual(1, NZ.xsize)
+    self.assertEqual(60, NZ.ysize)
+    self.assertEqual("nz", NZ.getAttribute("what/quantity"))
+
+    self.assertEqual(60,vp.getLevels())
+    self.assertEqual(200, vp.interval)
+    self.assertEqual(0, vp.minheight)
+    self.assertEqual(12000, vp.maxheight)
+    self.assertEqual(pvol.source, vp.source)
+    self.assertEqual(pvol.date, vp.date)
+    self.assertEqual(pvol.time, vp.time)
+
+    try:
+      self.assertEqual(1, uwnd.xsize)
+      self.assertEqual(60, uwnd.ysize)
+      self.assertEqual("UWND", uwnd.getAttribute("what/quantity"))
+    except:
+      exceptionTest = True
+
+    self.assertEqual(True, exceptionTest)
+
     
   def X_test_generate_2(self):
     pvol = _raveio.open(self.FIXTURE).object
-    generator = _wrwp.new()
-    generator.hmax = 2000
-    generator.dz = 200
+    wrwp = load_wrwp_defaults_to_obj()
+    wrwp.hmax = 2000
+    wrwp.dz = 200
+    fields = None
     
-    vp = generator.generate(pvol)
+    vp = wrwp.generate(pvol, fields)
 
     robj = _raveio.new()
     robj.object = vp
-    robj.save("slask.h5")
+    robj.save("slasktest.h5")
 
   def test_generate_with_several_howattributes(self):
     pvol = _raveio.open(self.FIXTURE2).object
-    generator = _wrwp.new()
-    generator.hmax = 2000
-    generator.dz = 200
+    wrwp = load_wrwp_defaults_to_obj()
+    wrwp.hmax = 2000
+    wrwp.dz = 200
+    fields = None
     
-    vp = generator.generate(pvol, "UWND,VWND,HGHT,NV")
+    fields = QUANTITIES
+    vp = wrwp.generate(pvol, fields)
     
     uwnd = vp.getUWND()
     vwnd = vp.getVWND()
@@ -195,7 +359,7 @@ class WrwpTest(unittest.TestCase):
     self.assertEqual(pvol.date, vp.date)
     self.assertEqual(pvol.time, vp.time)
     
-    self.assertEqual(900, vp.getAttribute("how/lowprf"))
+    '''self.assertEqual(900, vp.getAttribute("how/lowprf"))
     self.assertEqual(1200, vp.getAttribute("how/highprf"))
     self.assertAlmostEqual(0.61, vp.getAttribute("how/pulsewidth"), 4)
     self.assertAlmostEqual(5.35, vp.getAttribute("how/wavelength"), 4)
@@ -212,9 +376,9 @@ class WrwpTest(unittest.TestCase):
     self.assertAlmostEqual(0.2, vp.getAttribute("how/radomelossH"), 4)
     self.assertAlmostEqual(2.0, vp.getAttribute("how/rpm"), 4)
     self.assertEqual("PARTEC2", vp.getAttribute("how/software"))
-    self.assertEqual("ERIC", vp.getAttribute("how/system"))
-    self.assertEqual(4.0, vp.getAttribute("how/minrange"))
-    self.assertEqual(40.0, vp.getAttribute("how/maxrange"))
+    self.assertEqual("ERIC", vp.getAttribute("how/system"))'''
+    self.assertEqual(5.0, vp.getAttribute("how/minrange"))
+    self.assertEqual(25.0, vp.getAttribute("how/maxrange"))
 
     robj = _raveio.new()
     robj.object = vp
@@ -222,15 +386,15 @@ class WrwpTest(unittest.TestCase):
 
   def test_second_generate_with_several_howattributes(self):
     pvol = _raveio.open(self.FIXTURE3).object
-    generator = _wrwp.new()
-    generator.hmax = 12000
-    generator.dz = 200
-    generator.emin = 4.0
+    wrwp = load_wrwp_defaults_to_obj()
+    wrwp.emin = 4.0
+    fields = None
 
     exceptionTest = False
 
     try:
-      vp = generator.generate(pvol, "UWND,VWND,HGHT,NV")
+      fields = QUANTITIES
+      vp = wrwp.generate(pvol, fields)
 
       uwnd = vp.getUWND()
       vwnd = vp.getVWND()
@@ -261,7 +425,7 @@ class WrwpTest(unittest.TestCase):
       self.assertEqual(pvol.date, vp.date)
       self.assertEqual(pvol.time, vp.time)
     
-      self.assertEqual(450, vp.getAttribute("how/lowprf"))
+      '''self.assertEqual(450, vp.getAttribute("how/lowprf"))
       self.assertEqual(600, vp.getAttribute("how/highprf"))
       self.assertAlmostEqual(0.5, vp.getAttribute("how/pulsewidth"), 4)
       self.assertAlmostEqual(5.348660945892334, vp.getAttribute("how/wavelength"), 4)
@@ -277,9 +441,9 @@ class WrwpTest(unittest.TestCase):
       self.assertAlmostEqual(71.4227294921875, vp.getAttribute("how/radconstH"), 4)
       self.assertAlmostEqual(3.0, vp.getAttribute("how/rpm"), 4)
       self.assertEqual("EDGE", vp.getAttribute("how/software"))
-      self.assertEqual("EECDWSR-2501C-SDP", vp.getAttribute("how/system"))
-      self.assertEqual(4.0, vp.getAttribute("how/minrange"))
-      self.assertEqual(40.0, vp.getAttribute("how/maxrange"))
+      self.assertEqual("EECDWSR-2501C-SDP", vp.getAttribute("how/system"))'''
+      self.assertEqual(5.0, vp.getAttribute("how/minrange"))
+      self.assertEqual(25.0, vp.getAttribute("how/maxrange"))
 
       robj = _raveio.new()
       robj.object = vp
@@ -289,5 +453,60 @@ class WrwpTest(unittest.TestCase):
 
     self.assertEqual(True, exceptionTest)
 
+  def test_third_generate_with_several_howattributes(self):
+    pvol = _raveio.open(self.FIXTURE4).object
+    wrwp = load_wrwp_defaults_to_obj()
+    wrwp.emin = 4.0
+    fields = None
+
+    exceptionTest = False
+
+    try:
+      fields = QUANTITIES
+      vp = wrwp.generate(pvol, fields)
+
+      uwnd = vp.getUWND()
+      vwnd = vp.getVWND()
+      hght = vp.getHGHT()
+      nv = vp.getNV()
+
+      self.assertEqual(1, uwnd.xsize)
+      self.assertEqual(60, uwnd.ysize)
+      self.assertEqual("UWND", uwnd.getAttribute("what/quantity"))
+
+      self.assertEqual(1, vwnd.xsize)
+      self.assertEqual(60, vwnd.ysize)
+      self.assertEqual("VWND", vwnd.getAttribute("what/quantity"))
+
+      self.assertEqual(1, hght.xsize)
+      self.assertEqual(60, hght.ysize)
+      self.assertEqual("HGHT", hght.getAttribute("what/quantity"))
+
+      self.assertEqual(1, nv.xsize)
+      self.assertEqual(60, nv.ysize)
+      self.assertEqual("n", nv.getAttribute("what/quantity"))
+
+      self.assertEqual(60, vp.getLevels())
+      self.assertEqual(200, vp.interval)
+      self.assertEqual(0, vp.minheight)
+      self.assertEqual(12000, vp.maxheight)
+      self.assertEqual(pvol.source, vp.source)
+      self.assertEqual(pvol.date, vp.date)
+      self.assertEqual(pvol.time, vp.time)
+
+      self.assertEqual(5.0, vp.getAttribute("how/minrange"))
+      self.assertEqual(25.0, vp.getAttribute("how/maxrange"))
+      self.assertEqual("4.0,8.0,14.0,24.0,40.0", vp.getAttribute("how/angles"))
+      self.assertEqual("osd_zdr,osd_zdr_fastshort,osd_ldr_fastshort,osd_zdr_longshort", vp.getAttribute("how/task"))
+
+      robj = _raveio.new()
+      robj.object = vp
+      robj.save("slask3.h5")
+    except:
+      exceptionTest = True
+
+    self.assertEqual(False, exceptionTest)
+
+
 if __name__ == "__main__":
   unittest.main()
diff --git a/test/pytest/fixtures/seosd_qcvol_zdrvol_different_task.h5 b/test/pytest/fixtures/seosd_qcvol_zdrvol_different_task.h5
new file mode 100644 (file)
index 0000000..86fb6af
Binary files /dev/null and b/test/pytest/fixtures/seosd_qcvol_zdrvol_different_task.h5 differ
diff --git a/test/pytest/fixtures/wrwp_config.xml b/test/pytest/fixtures/wrwp_config.xml
new file mode 100644 (file)
index 0000000..aa69cd7
--- /dev/null
@@ -0,0 +1,63 @@
+<?xml version="1.0"?>
+<wrwp-params>
+    <param name="QUANTITIES">
+        <description>"Currently supported fields"</description>
+        <value>NV,HGHT,UWND,VWND,ff,ff_dev,dd,DBZH,DBZH_dev,NZ</value>
+    </param>
+    <param name="DMIN">
+        <description>Min distance for deriving a profile [m]</description>
+        <value>5000</value>
+    </param>
+    <param name="DMAX">
+        <description>Max distance for deriving a profile [m]</description>
+        <value>25000</value>
+    </param>
+    <param name="NMIN_WND">
+        <description>Minimum sample size for wind</description>
+        <value>40</value>
+    </param>
+    <param name="NMIN_REF">
+        <description>Minimum sample size for wind</description>
+        <value>40</value>
+    </param>
+    <param name="EMIN">
+        <description>Minimum elevation angle [deg]</description>
+        <value>0.5</value>
+    </param>
+    <param name="EMAX">
+        <description>Minimum elevation angle [deg]</description>
+        <value>45.0</value>
+    </param>
+    <param name="VMIN">
+        <description>Radial velocity threshold [m/s]</description>
+        <value>2.0</value>
+    </param>
+    <param name="FF_MAX">
+        <description>Maximum allowed layer velocity [m/s][</description>
+        <value>60.0</value>
+    </param>
+    <param name="DZ">
+        <description>Height interval for deriving a profile [m]</description>
+        <value>200</value>
+    </param>
+    <param name="HMAX">
+        <description>Maximum height of the profile [m]</description>
+        <value>12000</value>
+    </param>
+    <param name="NODATA_VP">
+        <description>Nodata value for profile</description>
+        <value>-9999</value>
+    </param>
+    <param name="UNDETECT_VP">
+        <description>Undetect value for profile</description>
+        <value>-9999</value>
+    </param>
+    <param name="GAIN_VP">
+        <description>Gain value for profile</description>
+        <value>1.0</value>
+    </param>
+    <param name="OFFSET_VP">
+        <description>Offset value for profile</description>
+        <value>0.0</value>
+    </param>
+</wrwp-params>