StationSampler Class Reference

#include <ERF_StationSampler.H>

Collaboration diagram for StationSampler:

Public Member Functions
	StationSampler (const std::string &pp_prefix)
bool	empty () const
amrex::Vector< Station > &	stations ()
const amrex::Vector< Station > &	stations () const
int	numColumns () const
void	buildColumns ()
void	appendRow (amrex::Real time, amrex::Real epoch_time, const amrex::Vector< amrex::Real > &values)
void	flush ()
void	setRestart (bool is_restart)
void	checkRestartFiles () const
void	setWriteTimestamp (bool write_timestamp, const std::string &datetime_format)
int	bufferSteps () const
bool	allowUnavailable () const
bool	lowHeightWarned () const
void	setLowHeightWarned ()
bool	levelsReported () const
void	setLevelsReported ()
bool	stencilsValidFor (const amrex::Vector< amrex::BoxArray > &ba, int finest_level) const
void	rememberStencilGrids (const amrex::Vector< amrex::BoxArray > &ba, int finest_level)

Private Member Functions
void	writeHeader (const Station &station, std::ostream &os) const
std::string	columnSignature (const Station &station) const
void	checkHeaderMatches (const Station &station, const std::string &filename) const

Private Attributes
amrex::Vector< Station >	m_stations
int	m_ncolumns = 0
int	m_buffer_steps = 100
bool	m_allow_unavailable = false
amrex::Vector< amrex::Real >	m_times
amrex::Vector< amrex::Real >	m_epoch_times
amrex::Vector< amrex::Real >	m_rows
amrex::Vector< amrex::BoxArray >	m_stencil_grids
bool	m_is_restart = false
bool	m_low_height_warned = false
bool	m_levels_reported = false
bool	m_write_timestamp = false
std::string	m_datetime_format
std::string	m_dir = "Output_Stations"
bool	m_opened = false

Constructor & Destructor Documentation

StationSampler()

StationSampler::StationSampler ( const std::string &  pp_prefix )

explicit

{
    ParmParse pp(pp_prefix);
 
    Vector<std::string> names;
    const int nstation = pp.countval("station_names");
    if (nstation > 0) {
        pp.getarr("station_names", names, 0, nstation);
    }
 
    pp.queryAdd("station_buffer_steps", m_buffer_steps);
    if (m_buffer_steps < 1) { m_buffer_steps = 1; }
 
    pp.queryAdd("station_output_dir", m_dir);
 
    pp.queryAdd("station_allow_unavailable", m_allow_unavailable);
 
    for (const auto& name : names)
    {
        check_station_name(name);
        for (const auto& earlier : m_stations) {
            if (earlier.name == name) {
                amrex::Abort("Station '" + name + "' is named twice in erf.station_names; the "
                             "two would write the same file");
            }
        }
 
        Station station;
        station.name = name;
 
        ParmParse pps(pp_prefix + "." + name);
 
        const int nfield = pps.countval("field");
        if (nfield <= 0) {
            Abort(station_key_error(name, "field", "must name at least one variable"));
        }
        Vector<std::string> fields;
        pps.getarr("field", fields, 0, nfield);
        for (const auto& field : fields) {
            StationVar var;
            var.name = field;
            station.vars.push_back(var);
        }
 
        // A station is placed either in lat/lon or in domain coordinates, never both.
        const std::string lon_key = pps.contains("long") ? "long" : "lon";
        const bool has_lat = pps.contains("lat");
        const bool has_lon = pps.contains("long") || pps.contains("lon");
        const bool has_x   = pps.contains("x");
        const bool has_y   = pps.contains("y");
 
        if ((has_lat || has_lon) && (has_x || has_y)) {
            Abort("Station '" + name + "': specify either lat/long or x/y, not both");
        }
 
        if (has_lat || has_lon) {
            if (!has_lat || !has_lon) {
                Abort("Station '" + name + "': lat and long must both be given");
            }
            const int nlat = pps.countval("lat");
            const int nlon = pps.countval(lon_key.c_str());
            if (nlat != nlon) {
                Abort("Station '" + name + "': lat has " + std::to_string(nlat) +
                      " values but " + lon_key + " has " + std::to_string(nlon) +
                      "; they are paired, so the counts must match");
            }
            Vector<Real> lat, lon;
            pps.getarr("lat", lat, 0, nlat);
            pps.getarr(lon_key.c_str(), lon, 0, nlon);
            for (int i = 0; i < nlat; ++i) {
                StationLoc loc;
                loc.use_latlon = true;
                loc.req_lat = lat[i];
                loc.req_lon = lon[i];
                station.locs.push_back(loc);
            }
        } else {
            if (!has_x || !has_y) {
                Abort("Station '" + name + "': give either lat/long, or x and y in domain coordinates");
            }
            const int nx = pps.countval("x");
            const int ny = pps.countval("y");
            if (nx != ny) {
                Abort("Station '" + name + "': x has " + std::to_string(nx) +
                      " values but y has " + std::to_string(ny) +
                      "; they are paired, so the counts must match");
            }
            Vector<Real> xv, yv;
            pps.getarr("x", xv, 0, nx);
            pps.getarr("y", yv, 0, ny);
            for (int i = 0; i < nx; ++i) {
                StationLoc loc;
                loc.use_latlon = false;
                loc.x = xv[i];
                loc.y = yv[i];
                station.locs.push_back(loc);
            }
        }
 
        // Heights come in one of two ways, and the difference is only what the
        // zero of the column is: .height_agl measures from the local terrain,
        // .height_abs from the bottom of the domain, in the same z the geometry
        // is given in.  Both are read into the one list the rest of the sampler
        // uses, with a flag saying which zero it is, so nothing downstream has
        // to ask which key the user typed.
        const int n_agl = pps.countval("height_agl");
        const int n_abs = pps.countval("height_abs");
 
        if (n_agl > 0 && n_abs > 0) {
            Abort("Station '" + name + "': give heights as either erf." + name +
                  ".height_agl (above the local terrain) or erf." + name +
                  ".height_abs (in the model's z coordinate), not both");
        }
        if (pps.contains("height")) {
            Abort("Station '" + name + "': erf." + name + ".height is not a key; say "
                  "erf." + name + ".height_agl for metres above the local terrain, or erf." +
                  name + ".height_abs for metres in the model's z coordinate");
        }
 
        station.heights_are_agl = (n_abs == 0);
        const char* height_key  = station.heights_are_agl ? "height_agl" : "height_abs";
        const int   nheight     = std::max(n_agl, n_abs);
        if (nheight > 0) {
            pps.getarr(height_key, station.heights, 0, nheight);
        }
 
        m_stations.push_back(std::move(station));
    }
}

Here is the call graph for this function:

Member Function Documentation

allowUnavailable()

bool StationSampler::allowUnavailable ( ) const

inline

{ return m_allow_unavailable; }

appendRow()

void StationSampler::appendRow	(	amrex::Real	time,
		amrex::Real	epoch_time,
		const amrex::Vector< amrex::Real > &	values
	)

{
    AMREX_ALWAYS_ASSERT(static_cast<int>(values.size()) == m_ncolumns);
 
    if (ParallelDescriptor::IOProcessor()) {
        m_times.push_back(time);
        m_epoch_times.push_back(epoch_time);
        m_rows.insert(m_rows.end(), values.begin(), values.end());
 
        if (static_cast<int>(m_times.size()) >= m_buffer_steps) { flush(); }
    }
}

Here is the call graph for this function:

bufferSteps()

int StationSampler::bufferSteps ( ) const

inline

{ return m_buffer_steps; }

buildColumns()

void StationSampler::buildColumns

(

)

{
    m_ncolumns = 0;
    for (auto& station : m_stations) {
        station.columns.clear();
        for (int il = 0; il < static_cast<int>(station.locs.size()); ++il) {
            station.locs[il].col_begin = static_cast<int>(station.columns.size());
            for (int iv = 0; iv < static_cast<int>(station.vars.size()); ++iv) {
                if (station.vars[iv].is_2d) { station.columns.push_back(m_ncolumns++); }
            }
            for (int ih = 0; ih < static_cast<int>(station.heights.size()); ++ih) {
                for (int iv = 0; iv < static_cast<int>(station.vars.size()); ++iv) {
                    if (!station.vars[iv].is_2d) { station.columns.push_back(m_ncolumns++); }
                }
            }
        }
    }
}

checkHeaderMatches()

void StationSampler::checkHeaderMatches ( const Station &  station, const std::string &  filename  ) const

private

{
    const std::vector<std::string> found = read_station_header(filename);
    const std::string found_format = find_format_line(found);
    const std::string want_format  = columnSignature(station);
 
    if (found_format == want_format) { return; }
 
    if (found_format.empty()) {
        Abort("Station '" + station.name + "': " + filename + " has no '# format:' line, so it "
              "was not written by this version of ERF's station output and cannot be appended "
              "to.  Move it aside and start a new series.");
    }
 
    // A tag mismatch is a different thing from a hash mismatch and deserves a
    // different message: nothing about the run is wrong, the file is just old.
    const std::string found_tag = found_format.substr(0, found_format.find(' '));
    if (found_tag != std::string(station_format_tag)) {
        Abort("Station '" + station.name + "': " + filename + " is in station file format '" +
              found_tag + "' and this ERF writes '" + station_format_tag + "'.  Move it aside "
              "and start a new series.");
    }
 
    // Same format, different columns.  Point at the first column line that
    // differs, which is the useful thing to say; the signature is what decided.
    std::ostringstream expected;
    writeHeader(station, expected);
    std::vector<std::string> want;
    {
        std::istringstream is(expected.str());
        std::string line;
        while (std::getline(is, line)) { want.push_back(line); }
    }
 
    // The signature lines differ by construction, so they say nothing a reader
    // does not already know; what is worth reporting is the first column that
    // differs.
    auto without_signature = [](const std::vector<std::string>& lines) {
        std::vector<std::string> out;
        for (const auto& line : lines) {
            if (line.compare(0, 10, "# format: ") != 0) { out.push_back(line); }
        }
        return out;
    };
    const std::vector<std::string> want_cmp  = without_signature(want);
    const std::vector<std::string> found_cmp = without_signature(found);
 
    std::string detail = "\n  the difference is not one the header spells out";
    for (std::size_t i = 0; i < want_cmp.size() || i < found_cmp.size(); ++i) {
        const bool have_want  = (i < want_cmp.size());
        const bool have_found = (i < found_cmp.size());
        if (have_want && have_found && want_cmp[i] == found_cmp[i]) { continue; }
        if (have_want && have_found) {
            detail = "\n  the file says:        " + found_cmp[i] +
                     "\n  this run would write: " + want_cmp[i];
        } else if (have_found) {
            detail = "\n  the file has a column this run does not: " + found_cmp[i];
        } else {
            detail = "\n  this run has a column the file does not: " + want_cmp[i];
        }
        break;
    }
 
    Abort("Station '" + station.name + "': " + filename + " describes a different set of "
          "columns than this restart would write -- a changed field, location, height or "
          "units list -- so appending to it would produce columns its header does not "
          "describe.  Change the configuration back, or move the file aside and start a "
          "new series." + detail);
}

checkRestartFiles()

void StationSampler::checkRestartFiles

(

)

const

{
    if (!m_is_restart) { return; }
    if (!ParallelDescriptor::IOProcessor()) { return; }
 
    for (const auto& station : m_stations) {
        const std::string filename = m_dir + "/" + station.name + ".dat";
        if (FileExists(filename)) { checkHeaderMatches(station, filename); }
    }
}

columnSignature()

std::string StationSampler::columnSignature ( const Station &  station ) const

private

{
    std::ostringstream os;
    os << std::setprecision(sigprecision);
    os << station_format_tag << '|' << station.name
       << "|ts=" << (m_write_timestamp ? 1 : 0);
 
    auto describe = [&](const StationVar& var, Real height, bool with_height)
    {
        os << '|' << var.name << ':' << var.units;
        if (with_height) {
            os << ':' << height << (station.heights_are_agl ? "agl" : "abs");
        } else {
            os << ":2d";
        }
        if (var.is_missing) { os << ":na"; }
    };
 
    for (const auto& loc : station.locs) {
        os << "|@";
        if (loc.use_latlon) {
            os << "ll=" << loc.req_lat << ',' << loc.req_lon;
        } else {
            os << "xy=" << loc.x << ',' << loc.y;
        }
        for (const auto& var : station.vars) {
            if (var.is_2d) { describe(var, Real(0.0), false); }
        }
        for (const auto height : station.heights) {
            for (const auto& var : station.vars) {
                if (!var.is_2d) { describe(var, height, true); }
            }
        }
    }
    return std::string(station_format_tag) + " " + signature_hash(os.str());
}

empty()

bool StationSampler::empty ( ) const

inline

{ return m_stations.empty(); }

flush()

void StationSampler::flush

(

)

{
    if (!ParallelDescriptor::IOProcessor()) { return; }
 
    const int nrow = static_cast<int>(m_times.size());
    if (nrow == 0) { return; }
 
    if (!m_opened) {
        if (!UtilCreateDirectory(m_dir, 0755)) { CreateDirectoryFailed(m_dir); }
    }
 
    for (const auto& station : m_stations)
    {
        const std::string filename = m_dir + "/" + station.name + ".dat";
 
        // A run that is restarting appends to the file it wrote before, so the
        // series is continuous; the restart is marked in a comment rather than
        // by starting a new file.  A run that is not restarting starts afresh.
        const bool append = m_opened || (m_is_restart && FileExists(filename));
 
        // The first thing a restart writes into an existing file: drop anything
        // the earlier run wrote past the checkpoint this run restarted from.
        // That the file is the same series was settled at setup, by
        // checkRestartFiles.
        if (append && !m_opened) {
            const int dropped = truncate_station_file(filename, m_times[0]);
            if (dropped > 0) {
                Print() << "Station output: dropped " << dropped << " row(s) at or after t = "
                        << m_times[0] << " from " << filename << ", written by the run before "
                        << "the restart past the checkpoint it restarted from" << std::endl;
            }
        }
 
        std::ofstream os(filename, append ? (std::ios::out | std::ios::app)
                                          : (std::ios::out | std::ios::trunc));
        if (!os.good()) { FileOpenFailed(filename); }
 
        if (!m_opened) {
            if (append) {
                os << "# restarted run resumes at t = "
                   << std::setprecision(timeprecision) << m_times[0] << "\n";
            } else {
                writeHeader(station, os);
            }
        }
 
        for (int irow = 0; irow < nrow; ++irow) {
            os << std::setw(datwidth) << std::setprecision(timeprecision) << m_times[irow];
            if (m_write_timestamp) {
                // Every column has to be one whitespace-delimited field for the
                // file to be readable by column, and the usual datetime formats
                // put a space between the date and the time.
                std::string stamp = getTimestamp(static_cast<double>(m_epoch_times[irow]),
                                                 m_datetime_format, false);
                std::replace(stamp.begin(), stamp.end(), ' ', 'T');
                os << std::setw(datwidth) << stamp;
            }
            for (const int icol : station.columns) {
                os << std::setw(datwidth) << std::setprecision(datprecision)
                   << m_rows[static_cast<std::size_t>(irow)*m_ncolumns + icol];
            }
            os << "\n";
        }
        os.close();
    }
 
    m_opened = true;
    m_times.clear();
    m_epoch_times.clear();
    m_rows.clear();
}

Here is the call graph for this function:

levelsReported()

bool StationSampler::levelsReported ( ) const

inline

{ return m_levels_reported; }

lowHeightWarned()

bool StationSampler::lowHeightWarned ( ) const

inline

{ return m_low_height_warned; }

numColumns()

int StationSampler::numColumns ( ) const

inline

{ return m_ncolumns; }

rememberStencilGrids()

void StationSampler::rememberStencilGrids	(	const amrex::Vector< amrex::BoxArray > &	ba,
		int	finest_level
	)

{
    m_stencil_grids.resize(finest_level+1);
    for (int lev = 0; lev <= finest_level; ++lev) { m_stencil_grids[lev] = ba[lev]; }
}

setLevelsReported()

void StationSampler::setLevelsReported ( )

inline

{ m_levels_reported = true; }

setLowHeightWarned()

void StationSampler::setLowHeightWarned ( )

inline

{ m_low_height_warned = true; }

setRestart()

void StationSampler::setRestart ( bool  is_restart )

inline

{ m_is_restart = is_restart; }

setWriteTimestamp()

void StationSampler::setWriteTimestamp ( bool  write_timestamp, const std::string &  datetime_format  )

inline

                                                                                    {
        m_write_timestamp = write_timestamp;
        m_datetime_format = datetime_format;
    }

stations() [1/2]

amrex::Vector<Station>& StationSampler::stations ( )

inline

{ return m_stations; }

stations() [2/2]

const amrex::Vector<Station>& StationSampler::stations ( ) const

inline

{ return m_stations; }

stencilsValidFor()

bool StationSampler::stencilsValidFor	(	const amrex::Vector< amrex::BoxArray > &	ba,
		int	finest_level
	)		const

{
    if (static_cast<int>(m_stencil_grids.size()) != finest_level+1) { return false; }
    for (int lev = 0; lev <= finest_level; ++lev) {
        if (m_stencil_grids[lev] != ba[lev]) { return false; }
    }
    return true;
}

writeHeader()

void StationSampler::writeHeader ( const Station &  station, std::ostream &  os  ) const

private

{
    // Coordinates are written with enough digits that a reader can tell two
    // neighbouring stations apart
    os << std::setprecision(datprecision);
    os << "# ERF station time series\n";
    // The one line a restart compares.  Everything below it is for the reader.
    os << "# format: " << columnSignature(station) << "\n";
    os << "# station: " << station.name << "\n";
 
    int col = 1;
    os << "# column " << col++ << ": time [s]\n";
    if (m_write_timestamp) {
        os << "# column " << col++ << ": UTC timestamp (" << m_datetime_format
           << ", with 'T' where the format has a space, so the column is one field)\n";
    }
 
    const bool station_heights_are_agl = station.heights_are_agl;
 
    auto describe = [&](const StationLoc& loc, const StationVar& var, Real height, bool with_height)
    {
        os << "# column " << col++ << ": " << var.name;
        if (!var.units.empty()) { os << " [" << var.units << "]"; }
        if (loc.use_latlon) {
            os << " at requested lat=" << loc.req_lat << " long=" << loc.req_lon
               << " (sampled lat=" << loc.got_lat << " long=" << loc.got_lon << ")";
        }
        os << " at x=" << loc.x << " y=" << loc.y;
        if (with_height) {
            os << ", height=" << height
               << (station_heights_are_agl ? " m above local terrain" : " m (absolute, model z)");
            if (height < loc.first_cc_height) {
                os << " -- below the first cell centre, so sampled there, at height="
                   << loc.first_cc_height << " m when this file was started";
            }
        }
        if (var.is_missing) { os << "  [NOT AVAILABLE in this run: constant " << var.missing_value << "]"; }
        if (var.earth_wind != StationEarthWind::None) {
            os << " (earth-relative; grid x axis at "
               << std::atan2(loc.sin_alpha, loc.cos_alpha) * Real(180.0) / PI
               << " deg counterclockwise from east";
            if (var.earth_wind == StationEarthWind::WindDir) {
                os << "; degrees clockwise from north, wind from; " << var.missing_value << " when calm";
            }
            os << ")";
        }
        os << "\n";
    };
 
    for (const auto& loc : station.locs) {
        for (const auto& var : station.vars) {
            if (var.is_2d) { describe(loc, var, Real(0.0), false); }
        }
        for (const auto height : station.heights) {
            for (const auto& var : station.vars) {
                if (!var.is_2d) { describe(loc, var, height, true); }
            }
        }
    }
    os << "#\n";
    os << "# Values are bilinearly interpolated in the horizontal and linearly\n";
    os << "# interpolated in the vertical, from the finest level that covers the\n";
    os << "# interpolation stencil.\n";
    os << "# A requested height below the first cell centre of that level is\n";
    os << "# reported as the value at the first cell centre; the columns where that\n";
    os << "# happens say so above.  Their sampled height is the one when this file\n";
    os << "# was started, and can change if regridding changes the level sampled from.\n";
}

Member Data Documentation

m_allow_unavailable

bool StationSampler::m_allow_unavailable = false

private

Referenced by allowUnavailable().

m_buffer_steps

int StationSampler::m_buffer_steps = 100

private

Referenced by bufferSteps().

m_datetime_format

std::string StationSampler::m_datetime_format

private

Referenced by setWriteTimestamp().

m_dir

std::string StationSampler::m_dir = "Output_Stations"

private

m_epoch_times

amrex::Vector<amrex::Real> StationSampler::m_epoch_times

private

m_is_restart

bool StationSampler::m_is_restart = false

private

Referenced by setRestart().

m_levels_reported

bool StationSampler::m_levels_reported = false

private

Referenced by levelsReported(), and setLevelsReported().

m_low_height_warned

bool StationSampler::m_low_height_warned = false

private

Referenced by lowHeightWarned(), and setLowHeightWarned().

m_ncolumns

int StationSampler::m_ncolumns = 0

private

Referenced by numColumns().

m_opened

bool StationSampler::m_opened = false

private

m_rows

amrex::Vector<amrex::Real> StationSampler::m_rows

private

m_stations

amrex::Vector<Station> StationSampler::m_stations

private

Referenced by empty(), and stations().

m_stencil_grids

amrex::Vector<amrex::BoxArray> StationSampler::m_stencil_grids

private

m_times

amrex::Vector<amrex::Real> StationSampler::m_times

private

m_write_timestamp

bool StationSampler::m_write_timestamp = false

private

Referenced by setWriteTimestamp().

---

The documentation for this class was generated from the following files:

• Source/IO/ERF_StationSampler.H
• Source/IO/ERF_StationSampler.cpp