Merge pull request #111 from ecmwf/feature/fix-geopotential-in-upper-level-bufr-plots

Fix geopotential plotting on upper air charts

Author: iainrussell

.gitignore

@@ -29,6 +29,7 @@ build
 .ipynb_checkpoints
 
 .vscode
+.vscode/
 *.sublime-workspace
 notebook/.ipynb_checkpoints
 compile_commands.json
@@ -42,5 +43,6 @@ tmp.*
 .env
 *.swp
 
+
 # mac
 .DS_Store

share/magics/bufr_98.xml

@@ -103,6 +103,8 @@
   <type value='2'>
 	<subtype value='91'>pilot</subtype>
 	<subtype value='101' template='temp'>temp</subtype>
+    <subtype value='109' template='temp'>temp</subtype>
+    <subtype value='111' template='temp'>temp</subtype>
   </type>
 
   <type value='4'>

share/magics/obs.xml

@@ -106,12 +106,12 @@
 	<obs_identification row='-1' column='-1'/>
  </obs_template>
 
- <obs_template type='temp' rows='4' columns='5'>:
+ <obs_template type='temp' rows='4' columns='6'>:
 	<obs_station_ring row='0' column='0'/>
 	<obs_wind/>
 	<obs_temperature row='1' column='-1'/>
 	<obs_dewpoint row='-1' column='-1'/>
-	<obs_height row='1' column='1'/>
+	<obs_height row='1' column='4'/>
 	<obs_thickness row='0' column='1' justification='left'/>
 	<obs_time_plot row='0' column='-1' />
 	<obs_identification row='-2' column='-1'/>

src/decoders/ObsDecoder.cc

@@ -27,6 +27,11 @@
 #include "TextVisitor.h"
 #include "expat.h"
 
+#include <algorithm>
+#include <set>
+#include <string>
+#include <vector>
+
 using namespace magics;
 
 map<int, bool> multilevels_;
@@ -52,7 +57,7 @@ class BufrSubType : public map<string, string> {
 class BufrFamily : public map<string, BufrSubType> {
 public:
     BufrFamily(const string&);
-    ~BufrFamily(){};
+    ~BufrFamily() {};
 
     void currentType(const map<string, string>& def) { type_ = def.find("value")->second; }
 
@@ -176,7 +181,7 @@ static void XMLCALL character(void* userData, const char* s, int len) {
 
 BufrIdentifiers::BufrIdentifiers(int centre) : centre_(centre) {
     ostringstream file, deffile;
-    file    << buildSharePath("bufr_" + to_string(centre) + ".xml");
+    file << buildSharePath("bufr_" + to_string(centre) + ".xml");
     deffile << buildSharePath("bufr_98.xml");
     char buf[BUFSIZ];
     XML_Parser parser = XML_ParserCreate(NULL);
@@ -293,9 +298,9 @@ class BufrAccessor {
     BufrAccessor(const string& descriptor) : descriptor_(descriptor) {
         init();
         auto token = translator_.find(descriptor_);
-        eccodes_ = (token == translator_.end()) ? descriptor_ : token->second;
+        eccodes_   = (token == translator_.end()) ? descriptor_ : token->second;
     }
-    virtual ~BufrAccessor(){};
+    virtual ~BufrAccessor() {};
     virtual void operator()(const ObsDecoder&, MvObs&, string&) const {}
     virtual void operator()(const ObsDecoder&, MvObs& obs, double& val) const { val = obs.value(eccodes_); }
     virtual void print() {}
@@ -311,7 +316,7 @@ map<string, string> BufrAccessor::translator_;
 class BufrMultiValueAccessor : public BufrAccessor {
 public:
     BufrMultiValueAccessor(const string& descriptor, int index) : BufrAccessor(descriptor), index_(index) {}
-    virtual ~BufrMultiValueAccessor(){};
+    virtual ~BufrMultiValueAccessor() {};
     virtual void operator()(const ObsDecoder&, MvObs& obs, double& val) const {
         const BufrIdentifiers& table = BufrIdentTable::get(obs.originatingCentre());
         val                          = obs.valueByOccurrence(index_, descriptor_);
@@ -381,7 +386,7 @@ class BufrHighCloudHeightAccessor : public BufrMultiValueAccessor {
 class BufrTypeAccessor : public BufrAccessor {
 public:
     BufrTypeAccessor() { descriptor_ = "type"; };
-    virtual ~BufrTypeAccessor(){};
+    virtual ~BufrTypeAccessor() {};
     void operator()(const ObsDecoder&, MvObs& obs, string& val) const {
         val = obs.messageType();
 
@@ -459,7 +464,7 @@ class BufrMultiLevelAccessor : public BufrAccessor {
     void operator()(const ObsDecoder& decoder, MvObs& obs, double& val) const {
         const BufrIdentifiers& table = BufrIdentTable::get(obs.originatingCentre());
         int type                     = obs.messageType();
-        
+
         map<int, bool>::const_iterator multilevel = multilevels_.find(type);
         if (multilevel == multilevels_.end()) {
             MagLog::warning() << "BufrMultiLevelAccessor> Unknown observation type [" << val << "]\n";
@@ -468,7 +473,7 @@ class BufrMultiLevelAccessor : public BufrAccessor {
         if (type == 0 || type == 1) {
             // surface data
             val = obs.value(surface_);
-            if ( val == kBufrMissingValue ) {
+            if (val == kBufrMissingValue) {
                 // Trying the descriptor for multilevel: some centres are using it for surface.
                 val = obs.value(altitude_);
             }
@@ -492,8 +497,43 @@ class BufrMultiLevelAccessor : public BufrAccessor {
 
 class BufrGeopotentialAccessor : public BufrMultiLevelAccessor {
 public:
-    BufrGeopotentialAccessor() : BufrMultiLevelAccessor("geopotential") {}
+    BufrGeopotentialAccessor() : BufrMultiLevelAccessor("geopotential") {
+        keys_       = {"geopotential", "nonCoordinateGeopotential", "nonCoordinateGeopotentialHeight"};
+        heightKeys_ = {"nonCoordinateGeopotentialHeight"};
+    }
+
     virtual ~BufrGeopotentialAccessor() {}
+
+    void operator()(const ObsDecoder& decoder, MvObs& obs, double& val) const {
+        const BufrIdentifiers& table = BufrIdentTable::get(obs.originatingCentre());
+        int type                     = obs.messageType();
+
+        map<int, bool>::const_iterator multilevel = multilevels_.find(type);
+        if (multilevel == multilevels_.end()) {
+            MagLog::warning() << "BufrMultiLevelAccessor> Unknown observation type [" << val << "]\n";
+            val = kBufrMissingValue;
+        }
+
+        if (type != 0 && type != 1 && type != 4 && type != 5) {
+            for (auto key : keys_) {
+                val = obs.valueByPressureLevel(decoder.level_, key);
+                // std::cout << "level: " << decoder.level_ << " key: " << key << " val: " << val << std::endl;
+                if (val != kBufrMissingValue) {
+                    if (heightKeys_.find(key) != heightKeys_.end()) {
+                        val *= 9.81;
+                    }
+                    return;
+                }
+            };
+        }
+
+        BufrMultiLevelAccessor::operator()(decoder, obs, val);
+    }
+
+
+protected:
+    std::vector<std::string> keys_;
+    std::set<std::string> heightKeys_;
 };
 
 class BufrTemperatureAccessor : public BufrMultiLevelAccessor {
@@ -670,7 +710,6 @@ void ObsDecoder::customisedPoints(const Transformation& transformation, const st
 
     map<string, BufrAccessor*> accessors;
     for (std::set<string>::const_iterator token = tokens.begin(); token != tokens.end(); ++token) {
-
         try {
             BufrAccessor* accessor = SimpleObjectMaker<BufrAccessor>::create(*token);
             accessors.insert(make_pair(*token, accessor));
@@ -679,7 +718,6 @@ void ObsDecoder::customisedPoints(const Transformation& transformation, const st
             BufrAccessor* accessor = new BufrAccessor(*token);
             accessors.insert(make_pair(*token, accessor));
         }
-
     }
 
     // Create the type accessor!