src/Environment/metarproperties.cxx

   1 // metarproperties.cxx -- Parse a METAR and write properties
   2 //
   3 // Written by David Megginson, started May 2002.
   4 // Rewritten by Torsten Dreyer, August 2010
   5 //
   6 // Copyright (C) 2002  David Megginson - david@megginson.com
   7 //
   8 // This program is free software; you can redistribute it and/or
   9 // modify it under the terms of the GNU General Public License as
  10 // published by the Free Software Foundation; either version 2 of the
  11 // License, or (at your option) any later version.
  12 //
  13 // This program is distributed in the hope that it will be useful, but
  14 // WITHOUT ANY WARRANTY; without even the implied warranty of
  15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16 // General Public License for more details.
  17 //
  18 // You should have received a copy of the GNU General Public License
  19 // along with this program; if not, write to the Free Software
  20 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  21 //
  22
  23 #ifdef HAVE_CONFIG_H
  24 #  include <config.h>
  25 #endif
  26
  27 #include <cstring> // for strlen
  28
  29 #include "metarproperties.hxx"
  30 #include "fgmetar.hxx"
  31 #include "environment.hxx"
  32 #include "atmosphere.hxx"
  33 #include "metarairportfilter.hxx"
  34 #include <simgear/scene/sky/cloud.hxx>
  35 #include <simgear/structure/exception.hxx>
  36 #include <simgear/misc/strutils.hxx>
  37 #include <simgear/magvar/magvar.hxx>
  38 #include <simgear/timing/sg_time.hxx>
  39 #include <Main/fg_props.hxx>
  40
  41 using std::string;
  42
  43 namespace Environment {
  44
  45 static vector<string> coverage_string;
  46
  47 /**
  48  * @brief Helper class to wrap SGMagVar functionality and cache the variation and dip for
  49  *        a certain position.
  50  */
  51 class MagneticVariation : public SGMagVar {
  52 public:
  53   /**
  54    * Constructor
  55    */
  56   MagneticVariation() : _lat(1), _lon(1), _alt(1) {
  57     recalc( 0.0, 0.0, 0.0 );
  58   }
  59
  60   /**
  61    * @brief get the magnetic variation for a specific position at the current time
  62    * @param lon the positions longitude in degrees
  63    * @param lat the positions latitude in degrees
  64    * @param alt the positions height above MSL (aka altitude) in feet
  65    * @return the magnetic variation in degrees
  66    */
  67   double get_variation_deg( double lon, double lat, double alt );
  68
  69   /**
  70    * @brief get the magnetic dip for a specific position at the current time
  71    * @param lon the positions longitude in degrees
  72    * @param lat the positions latitude in degrees
  73    * @param alt the positions height above MSL (aka altitude) in feet
  74    * @return the magnetic dip in degrees
  75    */
  76   double get_dip_deg( double lon, double lat, double alt );
  77 private:
  78   void recalc( double lon, double lat, double alt );
  79   SGTime _time;
  80   double _lat, _lon, _alt;
  81 };
  82
  83 inline void MagneticVariation::recalc( double lon, double lat, double alt )
  84 {
  85   // calculation of magnetic variation is expensive. Cache the position
  86   // and perform this calculation only if it has changed
  87   if( _lon != lon || _lat != lat || _alt != alt ) {
  88     SG_LOG(SG_ENVIRONMENT, SG_DEBUG, "Recalculating magvar for lon=" << lon << ", lat=" << lat << ", alt=" << alt );
  89     _lon = lon;
  90     _lat = lat;
  91     _alt = alt;
  92
  93     SGGeod location(SGGeod::fromDegFt(lon, lat, alt));
  94    _time.update( location, 0, 0 );
  95     update( lon, lat, alt, _time.getJD() );
  96   }
  97 }
  98
  99 inline double MagneticVariation::get_variation_deg( double lon, double lat, double alt )
 100 {
 101   recalc( lon, lat, alt );
 102   return get_magvar() * SGD_RADIANS_TO_DEGREES;
 103 }
 104
 105 inline double MagneticVariation::get_dip_deg( double lon, double lat, double alt )
 106 {
 107   recalc( lon, lat, alt );
 108   return get_magdip() * SGD_RADIANS_TO_DEGREES;
 109 }
 110
 111 MetarProperties::MetarProperties( SGPropertyNode_ptr rootNode ) :
 112   _rootNode(rootNode),
 113   _metarValidNode( rootNode->getNode( "valid", true ) ),
 114   _station_elevation(0.0),
 115   _station_latitude(0.0),
 116   _station_longitude(0.0),
 117   _min_visibility(16000.0),
 118   _max_visibility(16000.0),
 119   _base_wind_dir(0),
 120   _base_wind_range_from(0),
 121   _base_wind_range_to(0),
 122   _wind_speed(0.0),
 123   _wind_from_north_fps(0.0),
 124   _wind_from_east_fps(0.0),
 125   _gusts(0.0),
 126   _temperature(0.0),
 127   _dewpoint(0.0),
 128   _humidity(0.0),
 129   _pressure(0.0),
 130   _sea_level_temperature(0.0),
 131   _sea_level_dewpoint(0.0),
 132   _sea_level_pressure(29.92),
 133   _rain(0.0),
 134   _hail(0.0),
 135   _snow(0.0),
 136   _snow_cover(false),
 137   _day(0),
 138   _hour(0),
 139   _minute(0),
 140   _cavok(false),
 141   _magneticVariation(new MagneticVariation())
 142 {
 143   // Hack to avoid static initialization order problems on OSX
 144   if( coverage_string.empty() ) {
 145     coverage_string.push_back(SGCloudLayer::SG_CLOUD_CLEAR_STRING);
 146     coverage_string.push_back(SGCloudLayer::SG_CLOUD_FEW_STRING);
 147     coverage_string.push_back(SGCloudLayer::SG_CLOUD_SCATTERED_STRING);
 148     coverage_string.push_back(SGCloudLayer::SG_CLOUD_BROKEN_STRING);
 149     coverage_string.push_back(SGCloudLayer::SG_CLOUD_OVERCAST_STRING);
 150   }
 151   // don't tie metar-valid, so listeners get triggered
 152   _metarValidNode->setBoolValue( false );
 153   _tiedProperties.setRoot( _rootNode );
 154   _tiedProperties.Tie("data", this, &MetarProperties::get_metar, &MetarProperties::set_metar );
 155   _tiedProperties.Tie("station-id", this, &MetarProperties::get_station_id, &MetarProperties::set_station_id );
 156   _tiedProperties.Tie("station-elevation-ft", &_station_elevation );
 157   _tiedProperties.Tie("station-latitude-deg", &_station_latitude );
 158   _tiedProperties.Tie("station-longitude-deg", &_station_longitude );
 159   _tiedProperties.Tie("station-magnetic-variation-deg", this, &MetarProperties::get_magnetic_variation_deg );
 160   _tiedProperties.Tie("station-magnetic-dip-deg", this, &MetarProperties::get_magnetic_dip_deg );
 161   _tiedProperties.Tie("min-visibility-m", &_min_visibility );
 162   _tiedProperties.Tie("max-visibility-m", &_max_visibility );
 163   _tiedProperties.Tie("base-wind-range-from", &_base_wind_range_from );
 164   _tiedProperties.Tie("base-wind-range-to", &_base_wind_range_to );
 165   _tiedProperties.Tie("base-wind-speed-kt", this, &MetarProperties::get_wind_speed, &MetarProperties::set_wind_speed );
 166   _tiedProperties.Tie("base-wind-dir-deg", this, &MetarProperties::get_base_wind_dir, &MetarProperties::set_base_wind_dir );
 167   _tiedProperties.Tie("base-wind-from-north-fps", this, &MetarProperties::get_wind_from_north_fps, &MetarProperties::set_wind_from_north_fps );
 168   _tiedProperties.Tie("base-wind-from-east-fps",this, &MetarProperties::get_wind_from_east_fps, &MetarProperties::set_wind_from_east_fps );
 169   _tiedProperties.Tie("gust-wind-speed-kt", &_gusts );
 170   _tiedProperties.Tie("temperature-degc", &_temperature );
 171   _tiedProperties.Tie("dewpoint-degc", &_dewpoint );
 172   _tiedProperties.Tie("rel-humidity-norm", &_humidity );
 173   _tiedProperties.Tie("pressure-inhg", &_pressure );
 174   _tiedProperties.Tie("temperature-sea-level-degc", &_sea_level_temperature );
 175   _tiedProperties.Tie("dewpoint-sea-level-degc", &_sea_level_dewpoint );
 176   _tiedProperties.Tie("pressure-sea-level-inhg", &_sea_level_pressure );
 177   _tiedProperties.Tie("rain-norm", &_rain );
 178   _tiedProperties.Tie("hail-norm", &_hail );
 179   _tiedProperties.Tie("snow-norm", &_snow);
 180   _tiedProperties.Tie("snow-cover", &_snow_cover );
 181   _tiedProperties.Tie("day", &_day );
 182   _tiedProperties.Tie("hour", &_hour );
 183   _tiedProperties.Tie("minute", &_minute );
 184   _tiedProperties.Tie("decoded", this, &MetarProperties::get_decoded );
 185   _tiedProperties.Tie("cavok", &_cavok );
 186 }
 187
 188 MetarProperties::~MetarProperties()
 189 {
 190   delete _magneticVariation;
 191 }
 192
 193 void MetarProperties::invalidate()
 194 {
 195   if( _metarValidNode->getBoolValue() )
 196     _metarValidNode->setBoolValue(false);
 197 }
 198
 199 static const double thickness_value[] = { 0, 65, 600, 750, 1000 };
 200
 201 const char* MetarProperties::get_metar() const
 202 {
 203     if (!_metar) return "";
 204     return _metar->getData();
 205 }
 206
 207 void MetarProperties::set_metar( const char * metarString )
 208 {
 209     SGSharedPtr<FGMetar> m;
 210     if ((metarString == NULL) || (strlen(metarString) == 0)) {
 211         setMetar(m);
 212         return;
 213     }
 214
 215     try {
 216         m = new FGMetar( metarString );
 217     }
 218     catch( sg_io_exception ) {
 219         SG_LOG( SG_ENVIRONMENT, SG_WARN, "Can't parse metar: " << metarString );
 220         _metarValidNode->setBoolValue(false);
 221         return;
 222     }
 223
 224     setMetar(m);
 225 }
 226
 227 void MetarProperties::setMetar( SGSharedPtr<FGMetar> m )
 228 {
 229     _metar = m;
 230     _decoded.clear();
 231     if (!m) {
 232         return;
 233     }
 234
 235     const vector<string> weather = m->getWeather();
 236     for( vector<string>::const_iterator it = weather.begin(); it != weather.end(); ++it ) {
 237         if( false == _decoded.empty() ) _decoded.append(", ");
 238         _decoded.append(*it);
 239     }
 240
 241     _min_visibility = m->getMinVisibility().getVisibility_m();
 242     _max_visibility = m->getMaxVisibility().getVisibility_m();
 243
 244     const SGMetarVisibility *dirvis = m->getDirVisibility();
 245     for ( int i = 0; i < 8; i++, dirvis++) {
 246         SGPropertyNode *vis = _rootNode->getChild("visibility", i, true);
 247         double v = dirvis->getVisibility_m();
 248
 249         vis->setDoubleValue("min-m", v);
 250         vis->setDoubleValue("max-m", v);
 251     }
 252
 253     set_base_wind_dir(m->getWindDir());
 254     _base_wind_range_from = m->getWindRangeFrom();
 255     _base_wind_range_to = m->getWindRangeTo();
 256     set_wind_speed(m->getWindSpeed_kt());
 257
 258     _gusts = m->getGustSpeed_kt();
 259     _temperature = m->getTemperature_C();
 260     _dewpoint = m->getDewpoint_C();
 261     _humidity = m->getRelHumidity();
 262     _pressure = m->getPressure_inHg();
 263
 264     {
 265         // 1. check the id given in the metar
 266         FGAirport* a = FGAirport::findByIdent(m->getId());
 267
 268         // 2. if unknown, find closest airport with metar to current position
 269         if( a == NULL ) {
 270             SGGeod pos = SGGeod::fromDeg(
 271                 fgGetDouble( "/position/longitude-deg", 0.0 ),
 272                 fgGetDouble( "/position/latitude-deg", 0.0 ) );
 273             a = FGAirport::findClosest(pos, 10000.0, MetarAirportFilter::instance() );
 274         }
 275
 276         // 3. otherwise use ground elevation
 277         if( a != NULL ) {
 278             _station_elevation = a->getElevation();
 279             const SGGeod & towerPosition = a->getTowerLocation();
 280             _station_latitude = towerPosition.getLatitudeDeg();
 281             _station_longitude = towerPosition.getLongitudeDeg();
 282             _station_id = a->ident();
 283         } else {
 284             _station_elevation = fgGetDouble("/position/ground-elev-m", 0.0 ) * SG_METER_TO_FEET;
 285             _station_latitude = fgGetDouble( "/position/latitude-deg", 0.0 );
 286             _station_longitude = fgGetDouble( "/position/longitude-deg", 0.0 );
 287             _station_id = "XXXX";
 288         }
 289     }
 290
 291     {    // calculate sea level temperature, dewpoint and pressure
 292         FGEnvironment dummy; // instantiate a dummy so we can leech a method
 293         dummy.set_elevation_ft( _station_elevation );
 294         dummy.set_temperature_degc( _temperature );
 295         dummy.set_dewpoint_degc( _dewpoint );
 296         _sea_level_temperature = dummy.get_temperature_sea_level_degc();
 297         _sea_level_dewpoint = dummy.get_dewpoint_sea_level_degc();
 298
 299         double elevation_m = _station_elevation * SG_FEET_TO_METER;
 300         double fieldPressure = FGAtmo::fieldPressure( elevation_m, _pressure * atmodel::inHg );
 301         _sea_level_pressure = P_layer(0, elevation_m, fieldPressure, _temperature + atmodel::freezing, atmodel::ISA::lam0) / atmodel::inHg;
 302     }
 303
 304     bool isBC = false;
 305     bool isBR = false;
 306     bool isFG = false;
 307     bool isMI = false;
 308     bool isHZ = false;
 309
 310     {
 311         for( unsigned i = 0; i < 3; i++ ) {
 312             SGPropertyNode_ptr n = _rootNode->getChild("weather", i, true );
 313             vector<struct SGMetar::Weather> weather = m->getWeather2();
 314             struct SGMetar::Weather * w = i < weather.size() ? &weather[i] : NULL;
 315             n->getNode("intensity",true)->setIntValue( w != NULL ? w->intensity : 0 );
 316             n->getNode("vincinity",true)->setBoolValue( w != NULL ? w->vincinity : false );
 317             for( unsigned j = 0; j < 3; j++ ) {
 318
 319                 const string & phenomenon = w != NULL && j < w->phenomena.size() ? w->phenomena[j].c_str() : "";
 320                 n->getChild( "phenomenon", j, true )->setStringValue( phenomenon );
 321
 322                 const string & description = w != NULL && j < w->descriptions.size() ? w->descriptions[j].c_str() : "";
 323                 n->getChild( "description", j, true )->setStringValue( description );
 324
 325                 // need to know later,
 326                 // if its fog(FG) (might be shallow(MI) or patches(BC)) or haze (HZ) or mist(BR)
 327                 if( phenomenon == "FG" ) isFG = true;
 328                 if( phenomenon == "HZ" ) isHZ = true;
 329                 if( phenomenon == "BR" ) isBR = true;
 330                 if( description == "MI" ) isMI = true;
 331                 if( description == "BC" ) isBC = true;
 332             }
 333         }
 334     }
 335
 336     {
 337         static const char * LAYER = "layer";
 338         SGPropertyNode_ptr cloudsNode = _rootNode->getNode("clouds", true );
 339         const vector<SGMetarCloud> & metarClouds = m->getClouds();
 340         unsigned layerOffset = 0; // Oh, this is ugly!
 341
 342         // fog/mist/haze cloud layer does not work with 3d clouds yet :-(
 343         bool setGroundCloudLayer = _rootNode->getBoolValue("set-ground-cloud-layer", false ) &&
 344               false == (fgGetBool("/sim/rendering/shader-effects", false ) &&
 345                         fgGetBool("/sim/rendering/clouds3d-enable", false ) );
 346
 347         if( setGroundCloudLayer ) {
 348             // create a cloud layer #0 starting at the ground if its fog, mist or haze
 349
 350             // make sure layer actually starts at ground and set it's bottom at a constant
 351             // value below the station's elevation
 352             const double LAYER_BOTTOM_STATION_OFFSET =
 353               fgGetDouble( "/environment/params/fog-mist-haze-layer/offset-from-station-elevation-ft", -200 );
 354
 355             SGMetarCloud::Coverage coverage = SGMetarCloud::COVERAGE_NIL;
 356             double thickness = 0;
 357             double alpha = 1.0;
 358
 359             if( isFG ) { // fog
 360                 coverage = SGMetarCloud::getCoverage( isBC ?
 361                     fgGetString( "/environment/params/fog-mist-haze-layer/fog-bc-2dlayer-coverage", SGMetarCloud::COVERAGE_SCATTERED_STRING ) :
 362                     fgGetString( "/environment/params/fog-mist-haze-layer/fog-2dlayer-coverage", SGMetarCloud::COVERAGE_BROKEN_STRING )
 363                 );
 364
 365                 thickness = isMI ?
 366                    fgGetDouble("/environment/params/fog-mist-haze-layer/fog-shallow-thickness-ft",30) - LAYER_BOTTOM_STATION_OFFSET : // shallow fog, 10m/30ft
 367                    fgGetDouble("/environment/params/fog-mist-haze-layer/fog-thickness-ft",500) - LAYER_BOTTOM_STATION_OFFSET; // fog, 150m/500ft
 368                 alpha =  fgGetDouble("/environment/params/fog-mist-haze-layer/fog-2dlayer-alpha", 1.0);
 369             } else if( isBR ) { // mist
 370                 coverage = SGMetarCloud::getCoverage(fgGetString("/environment/params/fog-mist-haze-layer/mist-2dlayer-coverage", SGMetarCloud::COVERAGE_OVERCAST_STRING));
 371                 thickness =  fgGetDouble("/environment/params/fog-mist-haze-layer/mist-thickness-ft",2000) - LAYER_BOTTOM_STATION_OFFSET;
 372                 alpha =  fgGetDouble("/environment/params/fog-mist-haze-layer/mist-2dlayer-alpha",0.8);
 373             } else if( isHZ ) { // haze
 374                 coverage = SGMetarCloud::getCoverage(fgGetString("/environment/params/fog-mist-haze-layer/mist-2dlayer-coverage", SGMetarCloud::COVERAGE_OVERCAST_STRING));
 375                 thickness =  fgGetDouble("/environment/params/fog-mist-haze-layer/haze-thickness-ft",2000) - LAYER_BOTTOM_STATION_OFFSET;
 376                 alpha =  fgGetDouble("/environment/params/fog-mist-haze-layer/haze-2dlayer-alpha",0.6);
 377             }
 378
 379             if( coverage != SGMetarCloud::COVERAGE_NIL ) {
 380
 381                 // if there is a layer above the fog, limit the top to one foot below that layer's bottom
 382                 if( metarClouds.size() > 0 && metarClouds[0].getCoverage() != SGMetarCloud::COVERAGE_CLEAR )
 383                     thickness = metarClouds[0].getAltitude_ft() - LAYER_BOTTOM_STATION_OFFSET - 1;
 384
 385                 SGPropertyNode_ptr layerNode = cloudsNode->getChild(LAYER, 0, true );
 386                 layerNode->setDoubleValue( "coverage-type", SGCloudLayer::getCoverageType(coverage_string[coverage]) );
 387                 layerNode->setStringValue( "coverage", coverage_string[coverage] );
 388                 layerNode->setDoubleValue( "elevation-ft", _station_elevation + LAYER_BOTTOM_STATION_OFFSET );
 389                 layerNode->setDoubleValue( "thickness-ft", thickness );
 390                 layerNode->setDoubleValue( "visibility-m", _min_visibility );
 391                 layerNode->setDoubleValue( "alpha", alpha );
 392                 _min_visibility = _max_visibility =
 393                   fgGetDouble("/environment/params/fog-mist-haze-layer/visibility-above-layer-m",20000.0); // assume good visibility above the fog
 394                 layerOffset = 1;  // shudder
 395             }
 396         }
 397
 398         for( unsigned i = 0; i < 5-layerOffset; i++ ) {
 399             SGPropertyNode_ptr layerNode = cloudsNode->getChild(LAYER, i+layerOffset, true );
 400             SGMetarCloud::Coverage coverage = i < metarClouds.size() ? metarClouds[i].getCoverage() : SGMetarCloud::COVERAGE_CLEAR;
 401             double elevation =
 402                 i >= metarClouds.size() || coverage == SGMetarCloud::COVERAGE_CLEAR ?
 403                 -9999.0 :
 404                 metarClouds[i].getAltitude_ft() + _station_elevation;
 405
 406             layerNode->setDoubleValue( "alpha", 1.0 );
 407             layerNode->setStringValue( "coverage", coverage_string[coverage] );
 408             layerNode->setDoubleValue( "coverage-type", SGCloudLayer::getCoverageType(coverage_string[coverage]) );
 409             layerNode->setDoubleValue( "elevation-ft", elevation );
 410             layerNode->setDoubleValue( "thickness-ft", thickness_value[coverage]);
 411             layerNode->setDoubleValue( "span-m", 40000 );
 412             layerNode->setDoubleValue( "visibility-m", 50.0 );
 413         }
 414     }
 415
 416     _rain = m->getRain();
 417     _hail = m->getHail();
 418     _snow = m->getSnow();
 419     _snow_cover = m->getSnowCover();
 420     _day = m->getDay();
 421     _hour = m->getHour();
 422     _minute = m->getMinute();
 423     _cavok = m->getCAVOK();
 424     _metarValidNode->setBoolValue(true);
 425 }
 426
 427 void MetarProperties::setStationId( const std::string & value )
 428 {
 429     set_station_id(simgear::strutils::strip(value).c_str());
 430 }
 431
 432 double MetarProperties::get_magnetic_variation_deg() const
 433 {
 434   return _magneticVariation->get_variation_deg( _station_longitude, _station_latitude, _station_elevation );
 435 }
 436
 437 double MetarProperties::get_magnetic_dip_deg() const
 438 {
 439   return _magneticVariation->get_dip_deg( _station_longitude, _station_latitude, _station_elevation );
 440 }
 441
 442 static inline void calc_wind_hs( double north_fps, double east_fps, int & heading_deg, double & speed_kt )
 443 {
 444     speed_kt = sqrt((north_fps)*(north_fps)+(east_fps)*(east_fps)) * 3600.0 / (SG_NM_TO_METER * SG_METER_TO_FEET);
 445     heading_deg = SGMiscd::roundToInt(
 446         SGMiscd::normalizeAngle2( atan2( east_fps, north_fps ) ) * SGD_RADIANS_TO_DEGREES );
 447 }
 448
 449 void MetarProperties::set_wind_from_north_fps( double value )
 450 {
 451     _wind_from_north_fps = value;
 452     calc_wind_hs( _wind_from_north_fps, _wind_from_east_fps, _base_wind_dir, _wind_speed );
 453 }
 454
 455 void MetarProperties::set_wind_from_east_fps( double value )
 456 {
 457     _wind_from_east_fps = value;
 458     calc_wind_hs( _wind_from_north_fps, _wind_from_east_fps, _base_wind_dir, _wind_speed );
 459 }
 460
 461 static inline void calc_wind_ne( double heading_deg, double speed_kt, double & north_fps, double & east_fps )
 462 {
 463     double speed_fps = speed_kt * SG_NM_TO_METER * SG_METER_TO_FEET / 3600.0;
 464     north_fps = speed_fps * cos(heading_deg * SGD_DEGREES_TO_RADIANS);
 465     east_fps = speed_fps * sin(heading_deg * SGD_DEGREES_TO_RADIANS);
 466 }
 467
 468 void MetarProperties::set_base_wind_dir( double value )
 469 {
 470     _base_wind_dir = value;
 471     calc_wind_ne( (double)_base_wind_dir, _wind_speed, _wind_from_north_fps, _wind_from_east_fps );
 472 }
 473
 474 void MetarProperties::set_wind_speed( double value )
 475 {
 476     _wind_speed = value;
 477     calc_wind_ne( (double)_base_wind_dir, _wind_speed, _wind_from_north_fps, _wind_from_east_fps );
 478 }
 479
 480
 481 } // namespace Environment