WSL/SLF GitLab Repository

Commit 04bcddbd authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The Date::getJulianDate() called has been renamed as Date::getJulian() for simplicity.

Some extra documentation has been added for the UPS coordinate system. The coordinates test has been expanded. A missing default value has been added to the Suntrajectory.
parent 457d5e17
......@@ -229,7 +229,7 @@ double BufferedIOHandler::getAvgSamplingRate()
sum += (unsigned long)vec_buffer_meteo[ii].size();
}
if (sum > 0){
double days = buffer_end.getJulianDate() - buffer_start.getJulianDate();
double days = buffer_end.getJulian() - buffer_start.getJulian();
return ((double)sum / days);
}
}
......@@ -378,8 +378,8 @@ std::ostream& operator<<(std::ostream& os, const BufferedIOHandler& data)
os << "Config& cfg = " << hex << &data.cfg << dec << "\n";
os << "IOHandler &iohandler = " << hex << &data.iohandler << dec << "\n";
os << "Buffering " << data.chunks << " chunk(s) of " <<data.chunk_size.getJulianDate() << " day(s) with "
<< data.buff_before.getJulianDate() << " day(s) pre-buffering\n";
os << "Buffering " << data.chunks << " chunk(s) of " <<data.chunk_size.getJulian() << " day(s) with "
<< data.buff_before.getJulian() << " day(s) pre-buffering\n";
os << "Current buffer content (" << data.vec_buffer_meteo.size() << " stations, "
<< data.mapBufferedGrids.size() << " grids):\n";
......
......@@ -1243,7 +1243,8 @@ void Coords::UTM_to_WGS84(double east_in, double north_in, double& lat_out, doub
/**
* @brief Coordinate conversion: from WGS84 Lat/Long to Universal Polar Stereographic grid
* see J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2
* see J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2.
* This is valid above latitudes 84N or above 80S.
* @param[in] lat_in Decimal Latitude
* @param[in] long_in Decimal Longitude
* @param[out] east_out easting coordinate (Swiss system)
......@@ -1275,7 +1276,8 @@ void Coords::WGS84_to_UPS(double lat_in, double long_in, double& east_out, doubl
/**
* @brief Coordinate conversion: from Universal Polar Stereographic grid to WGS84 Lat/Long
* see J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2
* see J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2.
* This is valid above latitudes 84N or above 80S.
* @param[in] east_in easting coordinate (UTM)
* @param[in] north_in northing coordinate (UTM)
* @param[out] lat_out Decimal Latitude
......
......@@ -147,7 +147,7 @@ void Date::setDate(const Date& in_date)
dst = false;
undef = true;
} else {
setDate(in_date.getJulianDate(), in_date.getTimeZone(), in_date.getDST());
setDate(in_date.getJulian(), in_date.getTimeZone(), in_date.getDST());
}
}
......@@ -288,7 +288,7 @@ bool Date::getDST() const {
* @param gmt convert returned value to GMT? (default: false)
* @return julian date in the current timezone / in GMT depending on the gmt parameter
*/
double Date::getJulianDate(const bool& gmt) const {
double Date::getJulian(const bool& gmt) const {
if(undef==true)
throw UnknownValueException("Date object is undefined!", AT);
......@@ -750,7 +750,7 @@ std::ostream& operator<<(std::ostream &os, const Date &date) {
else {
os << date.toString(Date::ISO) << "\n";
os << "TZ=GMT" << showpos << date.timezone << noshowpos << "\t\t" << "DST=" << date.dst << "\n";
os << "julian:\t\t\t" << setprecision(10) << date.getJulianDate() << "\t(GMT=" << date.getJulianDate(true) << ")\n";
os << "julian:\t\t\t" << setprecision(10) << date.getJulian() << "\t(GMT=" << date.getJulian(true) << ")\n";
os << "ModifiedJulian:\t\t" << date.getModifiedJulianDate() << "\n";
os << "TruncatedJulian:\t" << date.getTruncatedJulianDate() << "\n";
os << "MatlabJulian:\t\t" << date.getMatlabDate() << "\n";
......
......@@ -114,7 +114,7 @@ class Date {
bool isUndef() const;
double getTimeZone() const;
bool getDST() const;
double getJulianDate(const bool& gmt=false) const;
double getJulian(const bool& gmt=false) const;
double getModifiedJulianDate(const bool& gmt=false) const;
double getTruncatedJulianDate(const bool& gmt=false) const;
time_t getUnixDate() const;
......
......@@ -488,7 +488,7 @@ std::string IOManager::toString() const {
os << point_cache.begin()->first.toString(Date::ISO) << " - 1 timestep\n";
}
if(count>1) {
const double avg_sampling = ( (point_cache.rbegin()->first.getJulianDate()) - (point_cache.begin()->first.getJulianDate()) ) / (double)(count-1);
const double avg_sampling = ( (point_cache.rbegin()->first.getJulian()) - (point_cache.begin()->first.getJulian()) ) / (double)(count-1);
os << "Resampled cache content (";
if(max_stations==min_stations)
......
......@@ -639,9 +639,9 @@ size_t IOUtils::seek(const Date& soughtdate, const std::vector<MeteoData>& vecM,
}
size_t first = 0, last = vecM.size()-1;
const double start_val=vecM[first].date.getJulianDate(true);
const double end_val=vecM[last].date.getJulianDate(true);
const double curr_val = soughtdate.getJulianDate(true);
const double start_val=vecM[first].date.getJulian(true);
const double end_val=vecM[last].date.getJulian(true);
const double curr_val = soughtdate.getJulian(true);
//if we reach this point: at least one element in buffer
if (vecM[first].date > soughtdate) {
......@@ -659,8 +659,8 @@ size_t IOUtils::seek(const Date& soughtdate, const std::vector<MeteoData>& vecM,
const double raw_pos = (curr_val-start_val) / (end_val-start_val);
const size_t start = MAX( (size_t)(floor(raw_pos*last*.8)), first);
const size_t end = MIN( (size_t)ceil(raw_pos*last*1.2), last);
if(vecM[start].date.getJulianDate(true)<curr_val) first=start;
if(vecM[end].date.getJulianDate(true)>=curr_val) last=end;
if(vecM[start].date.getJulian(true)<curr_val) first=start;
if(vecM[end].date.getJulian(true)>=curr_val) last=end;
//if we reach this point: the date is spanned by the buffer and there are at least two elements
if (exactmatch){
......
......@@ -164,7 +164,7 @@ void ResamplingAlgorithms::NearestNeighbour(const size_t& index, const Resamplin
const double& val1 = vecM[indexP1](paramindex);
const double& val2 = vecM[indexP2](paramindex);
if (IOUtils::checkEpsilonEquality(diff1.getJulianDate(true), diff2.getJulianDate(true), 0.1/1440)){ //within 6 seconds
if (IOUtils::checkEpsilonEquality(diff1.getJulian(true), diff2.getJulian(true), 0.1/1440)){ //within 6 seconds
md(paramindex) = Interpol1D::weightedMean(val1, val2, 0.5);
} else if (diff1 < diff2){
md(paramindex) = val1;
......@@ -255,11 +255,11 @@ void ResamplingAlgorithms::LinearResampling(const size_t& index, const Resamplin
//At this point indexP1 and indexP2 point to values that are different from IOUtils::nodata
const double& val1 = vecM[indexP1](paramindex);
const double jul1 = vecM[indexP1].date.getJulianDate(true);
const double jul1 = vecM[indexP1].date.getJulian(true);
const double& val2 = vecM[indexP2](paramindex);
const double jul2 = vecM[indexP2].date.getJulianDate(true);
const double jul2 = vecM[indexP2].date.getJulian(true);
md(paramindex) = linearInterpolation(jul1, val1, jul2, val2, resampling_date.getJulianDate(true));
md(paramindex) = linearInterpolation(jul1, val1, jul2, val2, resampling_date.getJulian(true));
}
/**
......@@ -313,7 +313,7 @@ void ResamplingAlgorithms::Accumulate(const size_t& index, const ResamplingPosit
//find start of accumulation period and initialize the sum
double sum = IOUtils::nodata;
const Date dateStart(resampling_date.getJulianDate() - accumulate_period/(24.*3600.), resampling_date.getTimeZone());
const Date dateStart(resampling_date.getJulian() - accumulate_period/(24.*3600.), resampling_date.getTimeZone());
bool found_start=false;
size_t start_idx; //this is the index of the first point of the window that contains dateStart
for (start_idx=index+1; (start_idx--) > 0; ) {
......@@ -378,13 +378,13 @@ double ResamplingAlgorithms::funcval(size_t pos, const size_t& paramindex, const
const double valend = vecM[end](paramindex);
if (valend == IOUtils::nodata) return IOUtils::nodata;
const double jul1 = vecM[pos].date.getJulianDate(true);
const double jul2 = vecM[end].date.getJulianDate(true);
const double jul1 = vecM[pos].date.getJulian(true);
const double jul2 = vecM[end].date.getJulian(true);
if(start_pt)
return valend - linearInterpolation(jul1, 0., jul2, valend, date.getJulianDate(true));
return valend - linearInterpolation(jul1, 0., jul2, valend, date.getJulian(true));
else
return linearInterpolation(jul1, 0., jul2, valend, date.getJulianDate(true));
return linearInterpolation(jul1, 0., jul2, valend, date.getJulian(true));
}
/**
......
......@@ -52,8 +52,8 @@ void FilterRate::process(const unsigned int& param, const std::vector<MeteoData>
double& curr_value = ovec[ii](param);
const double& prev_value = ovec[last_good](param);
const double curr_time = ovec[ii].date.getJulianDate();
const double prev_time = ovec[last_good].date.getJulianDate();
const double curr_time = ovec[ii].date.getJulian();
const double prev_time = ovec[last_good].date.getJulian();
if (curr_value == IOUtils::nodata)
continue;
......
......@@ -203,8 +203,8 @@ const ProcessingProperties& ProcessingBlock::getProperties() const {
}
std::ostream& operator<<(std::ostream& os, const ProcessingProperties& data) {
const double h_before = data.time_before.getJulianDate()*24.;
const double h_after = data.time_after.getJulianDate()*24.;
const double h_before = data.time_before.getJulian()*24.;
const double h_after = data.time_after.getJulian()*24.;
const unsigned int p_before = data.points_before;
const unsigned int p_after = data.points_after;
......
......@@ -65,7 +65,7 @@ class SunTrajectory {
static double getAngleOfIncidence(const double& sun_azi, const double& sun_elev,
const double& slope_azi, const double& slope_elev);
static double projectHorizontalToSlope(const double& sun_azi, const double& sun_elev,
const double& slope_azi, const double& slope_elev, const double& H_radiation, const double& elev_threshold);
const double& slope_azi, const double& slope_elev, const double& H_radiation, const double& elev_threshold=5.);
static double projectSlopeToHorizontal(const double& sun_azi, const double& sun_elev,
const double& slope_azi, const double& slope_elev, const double& S_radiation);
static double projectHorizontalToBeam(const double& sun_elev, const double& H_radiation);
......
......@@ -430,8 +430,8 @@ void A3DIO::read2DMeteo(std::vector< std::vector<MeteoData> >& vecMeteo)
//1D and 2D data must correspond, that means that if there is 1D data
//for a certain date (e.g. 1.1.2006) then 2D data must exist (prec2006.txt etc),
//otherwise throw FileNotFoundException
Date startDate(vecMeteo[0][0].date.getJulianDate(), in_tz, false); //so that the correct filenames for the TZ will be constructed
Date endDate(vecMeteo[0][vecMeteo[0].size()-1].date.getJulianDate(), in_tz, false);
Date startDate(vecMeteo[0][0].date.getJulian(), in_tz, false); //so that the correct filenames for the TZ will be constructed
Date endDate(vecMeteo[0][vecMeteo[0].size()-1].date.getJulian(), in_tz, false);
constructMeteo2DFilenames(startDate, endDate, filenames);//get all files for all years
const size_t stations = getNrOfStations(filenames, hashStations);
......@@ -816,7 +816,7 @@ int A3DIO::create1DFile(const std::vector< std::vector<MeteoData> >& data)
file.flags ( std::ios::fixed );
for(size_t j=0; j<size; j++) {
int yyyy, mm, dd, hh;
Date tmp_date(data[ii][j].date.getJulianDate(true), out_tz);
Date tmp_date(data[ii][j].date.getJulian(true), out_tz);
tmp_date.getDate(yyyy, mm, dd, hh);
file.fill('0');
file << setw(4) << yyyy << " " << setw(2) << mm << " " << setw(2) << dd << " " << setw(2) << hh << " ";
......@@ -910,14 +910,14 @@ int A3DIO::write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data,
std::ofstream file;
int startyear, year, month, day, hour;
Date tmp_date(data[0][0].date.getJulianDate(true), out_tz);
Date tmp_date(data[0][0].date.getJulian(true), out_tz);
tmp_date.getDate(startyear, month, day, hour);
open2DFile(data, fileprefix, label, startyear, file);
file.flags ( ios::fixed );
for(size_t ii=0; ii<nb_timesteps; ii++) {
tmp_date.setDate(data[0][ii].date.getJulianDate(true), out_tz);
tmp_date.setDate(data[0][ii].date.getJulian(true), out_tz);
tmp_date.getDate(year, month, day, hour);
if(year!=startyear) {
//if the year has changed, we need to write to a new file
......
......@@ -415,6 +415,7 @@ bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelT
Date base_date;
double P1, P2;
getDate(h, base_date, P1, P2);
std::cerr << "Found date=" << base_date.toString(Date::ISO) << "\n";
//see WMO code table5 for definitions of timeRangeIndicator. http://dss.ucar.edu/docs/formats/grib/gribdoc/timer.html
long timeRange;
......@@ -495,7 +496,7 @@ void GRIBIO::indexFile(const std::string& filename)
throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
}
llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this sets llcorner, cellsize and bearing_offset
//llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this sets llcorner, cellsize and bearing_offset
if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
throw InvalidArgumentException("Cells can not be represented by square cells. This is not supported!", AT);
}
......@@ -549,6 +550,9 @@ void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, con
indexFile(filename);
}
listFields(filename);
std::cerr << "fields have been listed\n";
//Basic meteo parameters
if(parameter==MeteoGrids::P) read2DGrid_indexed(1.2, 1, 0, date, grid_out); //PS
if(parameter==MeteoGrids::TA) read2DGrid_indexed(11.2, 105, 2, date, grid_out); //T_2M
......@@ -677,7 +681,7 @@ void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, con
void GRIBIO::readDEM(DEMObject& dem_out)
{
const Date d; //ie: undef. This will be caught when reading the GRIB file
const Date d(2012,9,19,0,0,0); //ie: undef. This will be caught when reading the GRIB file
std::string filename;
cfg.getValue("DEMFILE", "Input", filename);
......
......@@ -464,9 +464,9 @@ void ImisIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
vecMeteoAnetz.insert(vecMeteoAnetz.begin(), vecAnetzStation.size(), vector<MeteoData>());
//date_anetz_start/end must be changed to be a multiple of 6h before the original dateStart, dateEnd
Date date_anetz_start = Date(floor(dateStart.getJulianDate(true) * 4.0) / 4.0, 0.);
Date date_anetz_start = Date(floor(dateStart.getJulian(true) * 4.0) / 4.0, 0.);
date_anetz_start.setTimeZone(in_tz);
Date date_anetz_end = Date(floor(dateEnd.getJulianDate(true) * 4.0) / 4.0, 0.);
Date date_anetz_end = Date(floor(dateEnd.getJulian(true) * 4.0) / 4.0, 0.);
date_anetz_end.setTimeZone(in_tz);
//read Anetz Data
......@@ -512,7 +512,7 @@ void ImisIO::assimilateAnetzData(const Date& dateStart, const AnetzData& ad,
for (size_t jj=0; jj<vecMeteo[stationindex].size(); jj++){
while (vecMeteo[stationindex][jj].date > (current_slice_date+0.2485)){
counter++;
double julian = floor((current_slice_date.getJulianDate(true) +0.25001) * 4.0) / 4.0;
double julian = floor((current_slice_date.getJulian(true) +0.25001) * 4.0) / 4.0;
current_slice_date = Date(julian, 0.);
current_slice_date.setTimeZone(in_tz);
}
......@@ -569,7 +569,7 @@ void ImisIO::calculatePsum(const Date& dateStart, const Date& dateEnd,
const std::vector< std::vector<MeteoData> >& vecMeteoAnetz,
std::vector< std::vector<double> >& vec_of_psums)
{
const unsigned int nr_of_slices = (unsigned int)((dateEnd.getJulianDate(true) - dateStart.getJulianDate(true) + 0.00001) * 4.0) + 1;
const unsigned int nr_of_slices = (unsigned int)((dateEnd.getJulian(true) - dateStart.getJulian(true) + 0.00001) * 4.0) + 1;
for (size_t ii=0; ii<vecMeteoAnetz.size(); ii++){
double tmp_psum = 0.0;
......@@ -769,7 +769,7 @@ void ImisIO::readSWE(const Date& dateStart, const Date& dateEnd, std::vector< st
IOUtils::convertString(curr_swe, rs->getString(2));
if(prev_swe!=IOUtils::nodata && curr_swe!=IOUtils::nodata) {
//valid values for Delta computation
if((d.getJulianDate()-prev_date.getJulianDate())<=max_interval) {
if((d.getJulian()-prev_date.getJulian())<=max_interval) {
//data not too far apart, so we accept it for Delta SWE
//looking for matching timestamp in the vecMeteo
while(ii_serie<serie_len && vecMeteo[stationindex][ii_serie].date<d) ii_serie++;
......
......@@ -458,7 +458,7 @@ void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
if (myreader.contains_timestamp()){
myreader.read(dateStart.toString(Date::ISO), dateEnd.toString(Date::ISO), mytimestamps, mydata);
} else {
myreader.read(dateStart.getJulianDate(), dateEnd.getJulianDate(), mydata);
myreader.read(dateStart.getJulian(), dateEnd.getJulian(), mydata);
}
copy_data(myreader, mytimestamps, mydata, vecMeteo[ii]);
......@@ -516,9 +516,9 @@ void SMETIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMete
if(out_dflt_TZ!=IOUtils::nodata) {
Date tmp_date(vecMeteo[ii][jj].date);
tmp_date.setTimeZone(out_dflt_TZ);
julian = tmp_date.getJulianDate();
julian = tmp_date.getJulian();
} else {
julian = vecMeteo[ii][jj].date.getJulianDate();
julian = vecMeteo[ii][jj].date.getJulian();
}
vec_data.push_back(julian);
}
......
......@@ -608,7 +608,7 @@ void SNIO::writeStationMeteo(const std::vector<MeteoData>& vecmd, const std::str
Date tmp_date(vecmd[jj].date);
tmp_date.setTimeZone(out_tz);
tmp_date.getDate(YYYY, MM, DD, HH, MI);
const double sn_julian = tmp_date.getJulianDate() - sn_julian_offset + 0.5;
const double sn_julian = tmp_date.getJulian() - sn_julian_offset + 0.5;
const double ta = vecmd[jj](MeteoData::TA);
const double rh = vecmd[jj](MeteoData::RH);
const double hnw = vecmd[jj](MeteoData::HNW);
......
......@@ -65,5 +65,33 @@ int main() {
exit(1);
}
//Now to UPS
cout << "Lat/long to UPS\n";
point1.setProj("UPS","N");
point1.setLatLon(84.2872338889, -132.247989167, 0.);
const double ups_x=point1.getEasting(), ups_y=point1.getNorthing();
if(!IOUtils::checkEpsilonEquality(ups_x, 1530125.78038, 1e-4) || !IOUtils::checkEpsilonEquality(ups_y, 2426773.59547, 1e-4)) {
cerr << setprecision(12) << "calculated Easting=" << ups_x << " calculated northing=" << ups_y << "\n";
exit(1);
}
//Now to UPS
cout << "UPS to Lat/long\n";
point1.setProj("UPS","S");
point1.setXY(2500000., 1500000., 0.);
const double ups_lat=point1.getLat(), ups_lon=point1.getLon();
if(!IOUtils::checkEpsilonEquality(ups_lat, -83.6373175, 1e-4) || !IOUtils::checkEpsilonEquality(ups_lon, 135., 1e-4)) {
cerr << setprecision(12) << "calculated lat=" << ups_lat << " calculated lon=" << ups_lon << "\n";
exit(1);
}
cout << "Pretty printing of negative latitudes\n";
stringstream ss2;
ss2 << point1.printLatLon();
if(ss2.str() != string("(-83°38'14.343220\" , 135°0'0.000000\")")) {
cerr << "Pretty printing failed: " << ss2.str() << "\n";
exit(1);
}
return 0;
}
......@@ -32,7 +32,7 @@ bool writeSun24h(ofstream& os, const mio::Date start_date, const double iswr_ref
os << date.toString(Date::ISO) << "\t";
os << std::setprecision(10);
Sun.setDate(date.getJulianDate(), date.getTimeZone()); //local julian date and timezone
Sun.setDate(date.getJulian(), date.getTimeZone()); //local julian date and timezone
Sun.calculateRadiation(TA, RH, mean_albedo);
//radiation in the beam'
......
Supports Markdown
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