Commit 8bdcf826 authored by Lefebvre, Jordan's avatar Lefebvre, Jordan

Adding draft fmsscore. Prefixing HYSPLIT named utilities with haem to avoid confusion/conflict.

parent f60d2a48
Pipeline #111521 failed with stages
in 5 minutes and 32 seconds
TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
SUBPACKAGES_DIRS_CLASSIFICATIONS_OPTREQS
exe haemexe SS REQUIRED
model haemmodel SS REQUIRED
exe haemexe SS REQUIRED
)
......@@ -11,17 +11,24 @@ SET(SOURCE
#SET(HEADERS
#)
#TRIBITS_ADD_EXECUTABLE(radixsnd2arl
# NOEXEPREFIX
# NOEXESUFFIX
# SOURCES radixsnd2arl.cc
# INSTALLABLE
#)
TRIBITS_ADD_EXECUTABLE(haemsnd2arl
NOEXEPREFIX
NOEXESUFFIX
SOURCES haemsnd2arl.cc
INSTALLABLE
)
TRIBITS_ADD_EXECUTABLE(haemconprob
NOEXEPREFIX
NOEXESUFFIX
SOURCES haemconprob.cc
INSTALLABLE
)
TRIBITS_ADD_EXECUTABLE(conprob
TRIBITS_ADD_EXECUTABLE(fmsscore
NOEXEPREFIX
NOEXESUFFIX
SOURCES conprob.cc
SOURCES fmsscoreex.cc
INSTALLABLE
)
......
TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
#LIB_REQUIRED_PACKAGES aptoollib
#LIB_REQUIRED_TPLS Qt5Widgets Qt5Positioning Qt5Network Qt5Gui Qt5Core
#TEST_REQUIRED_PACKAGES testframework
LIB_REQUIRED_PACKAGES haemmodel
)
/*
* Example utility to calculate figure of merit in space (FMS) score between 2
* cdump files
*/
#include <algorithm>
#include <cmath>
#include <ctime>
#include <iostream>
#include <string>
#include <vector>
#include "radixio/fmsscore.hh"
using namespace radix;
int main(int argc, char **argv)
{
// std::cout << "************************" << std::endl;
// std::cout << "******* fmsscore *******" << std::endl;
// std::cout << "************************" << std::endl;
if (argc != 8)
{
std::cerr << "Usage: " << argv[0]
<< " <cdump_file_observed> <pollutant_index> "
"<cdump_file_calculated> <pollutant_index> <level_index> "
"<time_index> <scale_factor>"
<< std::endl;
return 1;
}
const char *fi[2] = {argv[1], argv[3]};
int polIndex0 = atoi(argv[2]);
int polIndex1 = atoi(argv[4]);
int levelIndex = atoi(argv[5]);
int timeIndex = atoi(argv[6]);
float scaleFactor = static_cast<float>(atof(argv[7]));
HysplitCDump::SP cd[2];
for (int i = 0; i < 2; i++)
{
cd[i] = std::make_shared<HysplitCDump>();
HysplitCDumpStream<HysplitCDump::SP> stream1(cd[i]);
if (!stream1.read_from(fi[i]))
{
std::cerr << "ERROR: Could not open file: " << fi[i] << std::endl;
return -1.f;
}
}
float score =
fmsscore(cd, polIndex0, polIndex1, levelIndex, timeIndex, scaleFactor);
if (score < 0.f)
{
std::cerr << "ERROR occured in fmsscore()" << std::endl;
return 1;
}
// std::cout << "FMS score = " << score << std::endl;
std::cout << score << std::endl;
return 0;
}
......@@ -13,6 +13,7 @@ concentrationcontrol.hh
concentrationworkflow.hh
control.hh
emission.hh
fmsscore.hh
meteorology.hh
outputgridrequest.hh
persistencemanager.hh
......@@ -33,6 +34,7 @@ concentrationcontrol.cc
concentrationworkflow.cc
control.cc
emission.cc
fmsscore.cc
meteorology.cc
outputgridrequest.cc
persistencemanager.cc
......
#include "haemmodel/fmsscore.hh"
namespace haem
{
// see: "User-oriented two-dimensional measure of effectiveness for the
// evaluation of transport and dispersion models." Warner, et al 2003.
//
// fi = filenames of 2 cdump files to score (observed first, calculated second)
// polIndex0, polIndex1 = pollutant indexes
// timeIndex = time index into cdump files
// levelIndex = time index into cdump files
// scale = scale factor by which to multiply concentrations from calculated
// cdump file (fi[1])
float fmsscore(HysplitCDump::SP *cd, int polIndex0, int polIndex1,
int timeIndex, int levelIndex, float scale)
{
// int numLocs = cd[0]->numLocations();
// if (cd[1]->numLocations() != numLocs)
// {
// std::cerr << "ERROR: cdump files have different num locations: " <<
// numLocs << ", " << cd[1]->numLocations() << std::endl;
// return -1.f;
// }
int numLat = cd[0]->numLat();
int numLon = cd[0]->numLon();
int numStartTimes = cd[0]->numStartTimes();
int numLevels = cd[0]->numLevels();
bool ok = true;
if (timeIndex < 0 || timeIndex >= numStartTimes)
{
std::cerr << "ERROR: invalid time index: " << timeIndex << ". "
<< "valid range is 0 to " << numStartTimes - 1 << std::endl;
ok = false;
}
if (levelIndex < 0 || levelIndex >= numLevels)
{
std::cerr << "ERROR: invalid level index: " << levelIndex << ". "
<< "valid range is 0 to " << numLevels - 1 << std::endl;
ok = false;
}
if (cd[1]->numLat() != numLat || cd[1]->numLon() != numLon)
{
std::cerr << "ERROR: cdump files have different lat x lon dimensions: "
<< numLat << " x " << numLon << ", " << cd[1]->numLat() << " x "
<< cd[1]->numLon() << std::endl;
ok = false;
}
if (polIndex0 < 0 || polIndex0 >= cd[0]->numPollutants())
{
std::cerr << "ERROR: invalid pollutant index for 1st cdump file"
<< std::endl;
ok = false;
}
if (polIndex1 < 0 || polIndex1 >= cd[1]->numPollutants())
{
std::cerr << "ERROR: invalid pollutant index for 1st cdump file"
<< std::endl;
ok = false;
}
if (!ok) return 1;
char *gridflags = new char[numLat * numLon];
float *gridf0 = new float[numLat * numLon];
float *gridf1 = new float[numLat * numLon];
for (int i = 0; i < numLat * numLon; i++)
{
gridflags[i] = 0;
gridf0[i] = 0.f;
gridf1[i] = 0.f;
}
for (int i = 0; i < 2; i++)
{
// std::cout << "num pollutants = " << cd[i]->numPollutants() << std::endl;
// NOTE: assuming pollutant index for for both files = 0
int pi = (0 == i) ? polIndex0 : polIndex1;
// std::cout << "deltaLat = " << cd[i]->deltaLat() << std::endl;
// std::cout << "deltaLon = " << cd[i]->deltaLon() << std::endl;
// std::cout << "cornerLat = " << cd[i]->cornerLat() << std::endl;
// std::cout << "cornerLon = " << cd[i]->cornerLon() << std::endl;
int nnzp = cd[i]->numNonZeroPoints(timeIndex, pi, levelIndex);
// std::cout << "nnzp = " << nnzp << std::endl;
for (int ni = 0; ni < nnzp; ++ni)
{
short xi, yi;
float conc;
cd[i]->point(timeIndex, pi, levelIndex, ni, xi, yi, conc);
if (i == 1) conc *= scale;
// printf("%d %d %f\n", xi, yi, (double)conc);
// set bit for this cdump
int index = yi * numLon + xi;
gridflags[index] |= 1 << i;
if (i == 0)
gridf0[index] = conc;
else if (i == 1)
gridf1[index] = conc;
}
}
float gridIntersection = 0, gridUnion = 0;
for (int i = 0; i < numLat * numLon; i++)
{
// calculate intersection when both bits set
// A_OV from Warner
if (gridflags[i] == 3)
{
if (gridf0[i] < gridf1[i])
gridIntersection += gridf0[i];
else
gridIntersection += gridf1[i];
}
// calculate union if at least one bit set
// A_OV + A_FN + A_FP from Warner
if (gridflags[i] > 0)
{
if (gridf0[i] > gridf1[i])
gridUnion += gridf0[i];
else
gridUnion += gridf1[i];
}
}
// std::cout << "intersection = " << gridIntersection << std::endl;
// std::cout << "union = " << gridUnion << std::endl;
delete[] gridflags;
delete[] gridf0;
delete[] gridf1;
return 100.f * gridIntersection / gridUnion;
}
} // namespace haem
#ifndef HAEM_HAEMMODEL_FMSSCORE_HH_
#define HAEM_HAEMMODEL_FMSSCORE_HH_
#include "radixcore/visibility.hh"
#include "radixio/hysplitcdump.hh"
namespace haem
{
RADIX_PUBLIC float fmsscore(HysplitCDump::SP *cd, int polIndex0, int polIndex1,
int timeIndex, int levelIndex, float scale);
}
#endif /** HAEM_HAEMMODEL_FMSSCORE_HH_ */
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment