radixsnd2arl.cc

/*
 * Example utility to convert a vertical profile of meteorological
 * data to the ARL format
 */

#include <algorithm>
#include <ctime>
#include <iostream>
#include <string>
#include <vector>

#include "radixcommand/commandline.hh"
#include "radixcore/stringfunctions.hh"
#include "radixmath/constants.hh"
#include "radixmath/util.hh"

#include "radixio/arldatastream.hh"
#include "radixio/csvfile.hh"

using namespace radix;

void addHour(struct tm *time, int hours)
{
  int seconds = hours * 60 * 60;

  time_t date_seconds = mktime(time) + seconds;

  *time = *localtime(&date_seconds);
}

/**
 * @brief interpolateValues Interpolate/extrapolate missing values from a data
 * vector Basic capability for now - linear interpolation, extrapolation assumes
 * constant from ends.
 * @param pressureValues Pressure vector for meteorology - requires all values
 * to be present. Used for determining 'interpolation bounds' (there's
 * definitely a better phrase for this)
 * @param valuesToInterpolate
 */
bool interpolateValues(const std::vector<float> &pressureValues,
                       std::vector<float> &valuesToInterpolate,
                       bool circular = false)
{
  float missingValue = -9999.f;

  // Ensure pressure and other vector are same length for interpolation
  if (pressureValues.size() != valuesToInterpolate.size())
  {
    std::cout << "Error! Pressure vector is not same size ("
              << pressureValues.size() << ") as vector to interpolate ("
              << valuesToInterpolate.size() << ") - can't use to interpolate"
              << std::endl;
    return false;
  }

  // Ensure we have no missing values in the pressure dataset
  for (float pressure : pressureValues)
  {
    if (pressure == -9999.f)
    {
      std::cout << "Error! Pressure vector contains " << missingValue
                << " values - can't use to interpolate" << std::endl;
      return false;
    }
  }

  // Write out initial values
  std::cout << "Interpolating " << valuesToInterpolate.size()
            << " values; initial:" << std::endl
            << "  ";
  for (float f : valuesToInterpolate)
  {
    std::cout << f << " ";
  }
  std::cout << std::endl;

  // Loop through vector first to find first/last non-missing indices
  size_t firstIndex = valuesToInterpolate.size(), lastIndex = 0;
  bool foundFirst = false;
  for (size_t i = 0; i < valuesToInterpolate.size(); ++i)
  {
    if (valuesToInterpolate[i] != missingValue)
    {
      if (!foundFirst)
      {
        firstIndex = i;
        foundFirst = true;
      }
      lastIndex = i;
    }
  }
  // Fill in values before first and after last indices
  std::cout << "  Found first (" << firstIndex << ") and last (" << lastIndex
            << ") non-missing indices; filling in values before & after..."
            << std::endl;
  for (size_t i = 0; i < firstIndex; ++i)
  {
    valuesToInterpolate[i] = valuesToInterpolate[firstIndex];
  }
  for (size_t i = lastIndex + 1; i < valuesToInterpolate.size(); ++i)
  {
    valuesToInterpolate[i] = valuesToInterpolate[lastIndex];
  }

  // Fill in missing values in central parts of vector
  std::cout << "  Filling in missing data in central part of vector..."
            << std::endl;
  for (size_t i = firstIndex + 1; i < lastIndex; ++i)
  {
    // Search for a missing value
    size_t lastGood = i - 1;
    if (valuesToInterpolate[i] == missingValue)
    {
      // Get the next good value
      while (valuesToInterpolate[i] == missingValue)
      {
        i++;
      }
      size_t nextGood = i;
      // Interpolate between the two good values
      for (size_t j = lastGood + 1; j < nextGood; ++j)
      {
        float lastGoodValue = valuesToInterpolate[lastGood],
              nextGoodValue = valuesToInterpolate[nextGood];

        if (circular && fabs(lastGoodValue - nextGoodValue) > 180.0)
        {
          std::cout << "    Circular interpolation with distance > 180 "
                       "degrees: performing correction"
                    << std::endl;
          if (lastGoodValue < nextGoodValue)
          {
            lastGoodValue += 360.0;
          }
          else
          {
            nextGoodValue += 360.0;
          }
        }

        valuesToInterpolate[j] =
            lastGoodValue +
            ((nextGoodValue - lastGoodValue) *
             ((pressureValues[j] - pressureValues[lastGood]) /
              (pressureValues[nextGood] - pressureValues[lastGood])));

        if (circular)
        {
          valuesToInterpolate[j] = fmod(valuesToInterpolate[j], 360.0);
        }
      }
    }
  }

  // Write out final values
  std::cout << "Interpolation complete; final:" << std::endl << "  ";
  for (float f : valuesToInterpolate)
  {
    std::cout << f << " ";
  }
  std::cout << std::endl;

  return true;
}

int main(int argc, char **argv)
{
  std::cout << "************************" << std::endl;
  std::cout << "***** radixsnd2arl *****" << std::endl;
  std::cout << "************************" << std::endl;

  // Set up command line options
  CommandLine commandLine(argc, argv);
  commandLine.addOption("i", "Input csv file containing met data", true);
  commandLine.addOption("clat", "Centre latitude of output ARL file (degrees)",
                        true);
  commandLine.addOption("clon", "Centre longitude of output ARL file (degrees)",
                        true);
  commandLine.addOption("e", "Extent of output ARL file (km) [500]", false);
  commandLine.addOption("r", "Resolution of output ARL file (km) [10]", false);
  commandLine.addOption("t", "Time of data start (YYYYMMDDHH) [1951111917]",
                        false);
  commandLine.addOption("n", "Number of one hour timesteps to output [1]",
                        false);
  commandLine.addOption("g", "Add ground level elevation to height values",
                        false);
  commandLine.addOption("o", "Output ARL file", false);

  // Ensure required options present
  std::vector<std::string> commandErrors;
  if (!commandLine.validate(commandErrors))
  {
    std::cout << "Error in arguments..." << std::endl;
    for (std::string error : commandErrors)
    {
      std::cout << "\t" << error << std::endl;
    }
    std::cout << std::endl;
    commandLine.printParsedLine(std::cout);

    return -1;
  }

  // Get command line options
  std::string inputCsvPath = commandLine.get<std::string>("i");
  std::string outputArlPath =
      commandLine.get<std::string>("o", inputCsvPath + ".bin");
  float extent          = commandLine.get<float>("e", 500.0);
  float resolution      = commandLine.get<float>("r", 10.0);
  float centreLat       = commandLine.get<float>("clat");
  float centreLon       = commandLine.get<float>("clon");
  int numberTimesteps   = commandLine.get<int>("n", 1);
  float groundElevation = commandLine.get<float>("g", 0.f);
  std::string startTime = commandLine.get<std::string>("t", "1951111917");
  // Get the grid size
  int numberGridCells = (int)(extent / resolution);
  if (!numberGridCells % 2 == 1)
  {
    numberGridCells++;
  }
  // Parse the start time
  int year  = from_string(startTime.substr(0, 4), 1951);
  int month = from_string(startTime.substr(4, 2), 11);
  int day   = from_string(startTime.substr(6, 2), 19);
  int hour  = from_string(startTime.substr(8, 2), 17);

  struct tm metTime = {0, 0, 0};
  metTime.tm_year   = year - 1900;
  metTime.tm_mon    = month - 1;
  metTime.tm_mday   = day;
  metTime.tm_hour   = hour;

  std::cout << "Creating ARL-formatted met file with parameters:" << std::endl;
  std::cout << "  " << extent << " by " << extent << " km grid centred on ("
            << centreLat << "," << centreLon << ") with resolution "
            << resolution << " (" << numberGridCells
            << " cells in each direction)" << std::endl;
  std::cout << "  " << numberTimesteps
            << " timesteps (of 1 hour each) starting from " << metTime.tm_year
            << "-" << metTime.tm_mon + 1 << "-" << metTime.tm_mday << " at "
            << metTime.tm_hour << "00" << std::endl
            << std::endl;

  // If we have correct options, parse input file
  std::cout << "Reading input csv file: " << inputCsvPath << std::endl;

  // Header strings to denote fields
  std::string pressureString = "pres", tempString = "temp", relHumString = "rh",
              dewPtString = "tdew", wDirString = "wdir", wSpdString = "wspd",
              heightString = "zhgt";
  int pressureIndex = -1, tempIndex = -1, relHumIndex = -1, dewPtIndex = -1,
      wDirIndex = -1, wSpdIndex = -1, heightIndex = -1;
  bool usingRelHum = false, usingDewPt = false;

  // Open and read the csv
  CSVFile inputCsv(inputCsvPath);
  std::vector<std::vector<std::string>> inputData;
  bool inputReadSuccess = inputCsv.data(inputData);
  if (!inputReadSuccess)
  {
    std::cout << "Failed to read input csv!" << std::endl;
    std::cout
        << "Please ensure the file exists and has the appropriate permissions."
        << std::endl;
    return -2;
  }
  // Ensure there is data in here
  if (inputData.size() < 5)
  {
    std::cout << "Not enough data in input csv!" << std::endl;
    std::cout
        << "There are only " << inputData.size()
        << " lines in this file - at least 5 are required (4 headers + data)"
        << std::endl;
    return -3;
  }

  // Work out which fields are which data points - should be 2nd line (1)
  std::cout << "Finding data fields..." << std::endl;
  for (int entry = 0; entry < inputData[1].size(); ++entry)
  {
    if (inputData[1][entry] == pressureString)
    {
      std::cout << "  Found pressure at index " << entry << std::endl;
      pressureIndex = entry;
    }
    else if (inputData[1][entry] == tempString)
    {
      std::cout << "  Found temperature at index " << entry << std::endl;
      tempIndex = entry;
    }
    else if (inputData[1][entry] == relHumString)
    {
      std::cout << "  Found relative humidity at index " << entry << std::endl;
      relHumIndex = entry;
      usingRelHum = true;
    }
    else if (inputData[1][entry] == dewPtString)
    {
      std::cout << "  Found dew point at index " << entry << std::endl;
      dewPtIndex = entry;
      usingDewPt = true;
    }
    else if (inputData[1][entry] == wDirString)
    {
      std::cout << "  Found wind direction at index " << entry << std::endl;
      wDirIndex = entry;
    }
    else if (inputData[1][entry] == wSpdString)
    {
      std::cout << "  Found wind speed at index " << entry << std::endl;
      wSpdIndex = entry;
    }
    else if (inputData[1][entry] == heightString)
    {
      std::cout << "  Found height at index " << entry << std::endl;
      heightIndex = entry;
    }
  }
  // If we are missing one, we can't continue
  if ((pressureIndex == -1) || (tempIndex == -1) ||
      (relHumIndex == -1 && dewPtIndex == -1) || (wDirIndex == -1) ||
      (wSpdIndex == -1))
  {
    std::cout << "Missing data fields in input!" << std::endl;
    if (pressureIndex == -1)
    {
      std::cout << "  Missing pressure field (denoted by " << pressureString
                << ")";
    }
    if (tempIndex == -1)
    {
      std::cout << "  Missing temperature field (denoted by " << tempString
                << ")";
    }
    if (relHumIndex == -1 && dewPtIndex == -1)
    {
      std::cout << "  Missing both relative humidity field (denoted by "
                << relHumString << ") and dew point field (denoted by "
                << dewPtString << ") - need at least one of these";
    }
    if (wDirIndex == -1)
    {
      std::cout << "  Missing wind direction field (denoted by " << wDirString
                << ")";
    }
    if (wSpdIndex == -1)
    {
      std::cout << "  Missing wind speed field (denoted by " << wSpdString
                << ")";
    }
    if (heightIndex == -1)
    {
      std::cout << "  Missing height field (denoted by " << heightString << ")";
    }
    std::cout << "Please ensure these fields are present and rerun.";
    return -4;
  }

  // Ensure we aren't using both relative humidity and dew point (duplicate
  // data) Prefer use of relative humidity
  if (usingRelHum)
  {
    usingDewPt = false;
  }

  // Read the data
  std::vector<float> inputPressures, inputTemps, inputRelHums, inputDewPts,
      inputWDirs, inputWSpds, inputHeights;
  std::cout << "Data fields found - reading data..." << std::endl;
  float missingValue = -9999.f;
  for (int row = 4; row < inputData.size(); ++row)
  {
    if (inputData[row].size() < inputData[1].size())
    {
      std::cout << "  Warning: this row (" << row + 1
                << ") has less entries than the field names row" << std::endl;
      ;
    }

    bool foundPressure = false, foundTemp = false, foundRelHum = false,
         foundDewPt = false, foundWDir = false, foundWSpd = false,
         foundHeight = false;
    for (int entry = 0; entry < inputData[row].size(); ++entry)
    {
      float thisValue = from_string(inputData[row][entry], missingValue);

      if (entry == pressureIndex)
      {
        foundPressure = true;
        inputPressures.push_back(thisValue);
      }
      else if (entry == tempIndex)
      {
        foundTemp = true;
        inputTemps.push_back(thisValue);
      }
      else if (usingRelHum && entry == relHumIndex)
      {
        foundRelHum = true;
        inputRelHums.push_back(thisValue);
      }
      else if (usingDewPt && entry == dewPtIndex)
      {
        foundDewPt = true;
        inputDewPts.push_back(thisValue);
      }
      else if (entry == wSpdIndex)
      {
        foundWSpd = true;
        inputWSpds.push_back(thisValue);
      }
      else if (entry == wDirIndex)
      {
        foundWDir = true;
        inputWDirs.push_back(thisValue);
      }
      else if (entry == heightIndex)
      {
        foundHeight = true;
        inputHeights.push_back(thisValue);
      }
    }
    // Add missing data if any of the above weren't found
    if (!foundPressure)
    {
      std::cout << "  Warning: couldn't find pressure in row " << row
                << std::endl;
      inputPressures.push_back(missingValue);
    }
    if (!foundTemp)
    {
      std::cout << "  Warning: couldn't find temperature in row " << row
                << std::endl;
      inputTemps.push_back(missingValue);
    }
    if (usingRelHum && !foundRelHum)
    {
      std::cout << "  Warning: couldn't find relative humidity in row " << row
                << std::endl;
      inputRelHums.push_back(missingValue);
    }
    if (usingDewPt && !foundDewPt)
    {
      std::cout << "  Warning: couldn't find dew point in row " << row
                << std::endl;
      inputDewPts.push_back(missingValue);
    }
    if (!foundWSpd)
    {
      std::cout << "  Warning: couldn't find wind speed in row " << row
                << std::endl;
      inputWSpds.push_back(missingValue);
    }
    if (!foundWDir)
    {
      std::cout << "  Warning: couldn't find wind direction in row " << row
                << std::endl;
      inputWDirs.push_back(missingValue);
    }
    if (!foundHeight)
    {
      std::cout << "  Warning: couldn't find height in row " << row
                << std::endl;
      inputHeights.push_back(missingValue);
    }
  }
  std::cout << "Data read complete: " << inputPressures.size()
            << " entries read." << std::endl;

  // Check if the grid is of appropriate size - with too small grid/too many
  // vertical levels, the index header can become oversized leading to errors
  int indexBase = 108, levelBase = 48;
  int indexSize  = indexBase + (inputPressures.size() * levelBase);
  int recordSize = numberGridCells * numberGridCells;
  if (indexSize > recordSize)
  {
    // The index header will be too big - instruct user to either reduce number
    // of vertical levels or increase number of grid cells
    std::cout << "Error: index header will be too large!" << std::endl
              << "  The size of the index header in the ARL-formatted output ("
              << indexBase << " + number of vertical levels * " << levelBase
              << " = " << indexSize
              << ") will be larger than the size of each record (number of "
                 "grid cells ^ 2 = "
              << recordSize << ")." << std::endl;
    std::cout << "  Please either increase the number of horizontal grid cells "
                 "(increase extent/resolution) or reduce the number of input "
                 "vertical levels."
              << std::endl;
    throw std::exception();
  }

  // Interpolate the data to remove missing values
  std::cout << "Interpolating values to remove missing data..." << std::endl;
  std::cout << "Temperature:" << std::endl;
  interpolateValues(inputPressures, inputTemps);
  if (usingRelHum)
  {
    std::cout << "Relative humidity:" << std::endl;
    interpolateValues(inputPressures, inputRelHums);
  }
  if (usingDewPt)
  {
    std::cout << "Dew point:" << std::endl;
    interpolateValues(inputPressures, inputDewPts);
  }
  std::cout << "Wind speed:" << std::endl;
  interpolateValues(inputPressures, inputWSpds);
  std::cout << "Wind direction:" << std::endl;
  interpolateValues(inputPressures, inputWDirs, true);
  std::cout << "Height:" << std::endl;
  interpolateValues(inputPressures, inputHeights);
  std::cout << "Interpolation complete." << std::endl;

  // Convert the data into ARL format
  std::cout << "Converting data to ARL format..." << std::endl;
  ARLDataStream outputStream(outputArlPath, std::ios::out);
  for (int timestep = 0; timestep < numberTimesteps; ++timestep)
  {
    // Write index header section
    std::cout << "    Writing index headers for timestep " << timestep << "..."
              << std::endl;
    ARLRecordHeader thisRecordHeader;
    thisRecordHeader.year  = metTime.tm_year;
    thisRecordHeader.month = metTime.tm_mon + 1;
    thisRecordHeader.day   = metTime.tm_mday;
    thisRecordHeader.hour  = metTime.tm_hour;
    thisRecordHeader.ic    = 0;
    thisRecordHeader.il    = 0;
    thisRecordHeader.cgrid = "99";
    thisRecordHeader.kvar  = "INDX";
    thisRecordHeader.nexp  = 0;
    thisRecordHeader.prec  = 0.f;
    thisRecordHeader.var1  = 0.f;
    ARLIndexHeader thisIndexHeader;
    thisIndexHeader.model_id = "NFDB";
    thisIndexHeader.icx      = 0;
    thisIndexHeader.mn       = 0;
    thisIndexHeader.pole_lat = centreLat;
    thisIndexHeader.pole_lon = centreLon;
    thisIndexHeader.ref_lat  = centreLat;
    thisIndexHeader.ref_lon  = centreLon;
    thisIndexHeader.size     = resolution;
    thisIndexHeader.orient   = 0.f;
    thisIndexHeader.tang_lat = centreLat;
    thisIndexHeader.sync_xp  = (numberGridCells + 1) / 2;
    thisIndexHeader.sync_yp  = (numberGridCells + 1) / 2;
    thisIndexHeader.sync_lat = centreLat;
    thisIndexHeader.sync_lon = centreLon;
    thisIndexHeader.dummy    = 0.f;
    thisIndexHeader.nx       = numberGridCells;
    thisIndexHeader.ny       = numberGridCells;
    thisIndexHeader.nz       = inputPressures.size();
    thisIndexHeader.z_flag   = 2;
    thisIndexHeader.lenh     = numberGridCells * numberGridCells;
    thisIndexHeader.levels   = inputPressures;
    thisIndexHeader.num_vars_at_levels =
        std::vector<int>(inputPressures.size(), 5);
    std::vector<std::string> surfaceVarNames = {"PRSS", "TEMP", "RELH", "UWND",
                                                "VWND"};
    std::vector<std::string> varNames        = {"HGTS", "TEMP", "RELH", "UWND",
                                         "VWND"};

    thisIndexHeader.var_names.push_back(surfaceVarNames);

    for (size_t level = 1; level < inputPressures.size(); ++level)
    {
      thisIndexHeader.var_names.push_back(varNames);
    }
    std::vector<int> checkSums = std::vector<int>(5, 0);
    for (size_t level = 0; level < inputPressures.size(); ++level)
    {
      thisIndexHeader.check_sums.push_back(checkSums);
    }

    // Write the headers
    outputStream.write_record_header(thisRecordHeader);
    outputStream.write_index_header(thisRecordHeader, thisIndexHeader);
    std::cout << "    Headers written" << std::endl;

    // Write meteorological variables for each level
    for (int level = 0; level < inputPressures.size(); ++level)
    {
      std::cout << "    Writing meteorological data for timestep " << timestep
                << ", level " << level << "..." << std::endl;

      thisRecordHeader.il = level;

      // Write pressure/height variables
      {
        if (level == 0)
        {
          thisRecordHeader.kvar = "PRSS";
          // Add ground elevation to height (default is 0)
          float thisHeight      = inputHeights[level] + groundElevation;
          thisRecordHeader.var1 = thisHeight;
          std::vector<std::vector<float>> thisData(
              numberGridCells, std::vector<float>(numberGridCells, thisHeight));
          std::cout << "      pressure...";
          outputStream.write_record_header(thisRecordHeader);
          outputStream.write_record(thisRecordHeader, thisIndexHeader,
                                    thisData);
          std::cout << "written" << std::endl;
        }
        else
        {
          thisRecordHeader.kvar = "HGTS";
          thisRecordHeader.var1 = inputHeights[level];
          std::vector<std::vector<float>> thisData(
              numberGridCells,
              std::vector<float>(numberGridCells, inputHeights[level]));
          std::cout << "      height...";
          outputStream.write_record_header(thisRecordHeader);
          outputStream.write_record(thisRecordHeader, thisIndexHeader,
                                    thisData);
          std::cout << "written" << std::endl;
        }
      }

      // Write temperature levels
      // First convert temperature from Celsius to Kelvin
      inputTemps[level] = inputTemps[level] - ABS_ZERO_CELSIUS;
      {
        thisRecordHeader.kvar = "TEMP";
        thisRecordHeader.var1 = inputTemps[level];
        std::vector<std::vector<float>> thisData(
            numberGridCells,
            std::vector<float>(numberGridCells, inputTemps[level]));
        std::cout << "      temperature...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }

      // Write relative humidity levels
      {
        double thisVar = 0.0;
        if (usingRelHum)
        {
          thisVar = inputRelHums[level];
        }
        else if (usingDewPt)
        {
          thisVar =
              dewPointToRelativeHumidity(inputDewPts[level], inputTemps[level]);
        }

        thisRecordHeader.kvar = "RELH";
        thisRecordHeader.var1 = thisVar;
        std::vector<std::vector<float>> thisData(
            numberGridCells, std::vector<float>(numberGridCells, thisVar));
        std::cout << "      relative humidity...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }

      // Write wind data
      // Need to calculate wind u and v components from direction/speed
      float thisWindU =
          (inputWSpds[level] * 0.5144444444) *
          sin(toRadians(fmod((inputWDirs[level] + 180.0), 360.0)));
      float thisWindV =
          (inputWSpds[level] * 0.5144444444) *
          cos(toRadians(fmod((inputWDirs[level] + 180.0), 360.0)));
      std::cout << "      Initial wspd = " << inputWSpds[level]
                << " knots, wdir: " << inputWDirs[level] << " degrees"
                << std::endl;
      std::cout << "      Converted wind components: u = " << thisWindU
                << " m/s, v = " << thisWindV << "m/s" << std::endl;
      {
        thisRecordHeader.kvar = "UWND";
        thisRecordHeader.var1 = thisWindU;
        std::vector<std::vector<float>> thisData(
            numberGridCells, std::vector<float>(numberGridCells, thisWindU));
        std::cout << "      wind (u component)...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }
      {
        thisRecordHeader.kvar = "VWND";
        thisRecordHeader.var1 = thisWindV;
        std::vector<std::vector<float>> thisData(
            numberGridCells, std::vector<float>(numberGridCells, thisWindV));
        std::cout << "      wind (v component)...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }
    }
    std::cout << " Written timestep " << timestep << std::endl;

    // Add an hour to time for next timestep
#ifdef _WIN32
    metTime.tm_year = metTime.tm_year + 400;
#endif  // _WIN32
    addHour(&metTime, 1);
#ifdef _WIN32
    metTime.tm_year = metTime.tm_year - 400;
#endif
  }

  std::cout << "File write complete!" << std::endl;

  return 0;
}