WSL/SLF GitLab Repository

Skip to content
GitLab
Commit d7fb729c authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Several warnings have been fixed. The (slightly) wrong geolocalization... 

Several warnings have been fixed. The (slightly) wrong geolocalization calculation such as found in fixed in NetCDF has also been fixed in GRIB. A bug that was preventing GRIB DEM from being read has been fixed.
parent 92454269
Hide whitespace changes
Inline Side-by-side
meteoio/plugins/GRIBIO.cc
View file @ d7fb729c
......@@ -100,8 +100,8 @@ GRIBIO::GRIBIO(const std::string& configfile)
meteo_ext(default_ext), grid2d_ext(default_ext), grid2d_prefix(), idx_filename(), coordin(), coordinparam(),
VW(), DW(), wind_date(), llcorner(), fp(NULL), idx(NULL),
latitudeOfNorthernPole(IOUtils::nodata), longitudeOfNorthernPole(IOUtils::nodata), bearing_offset(IOUtils::nodata),
cellsize_x(IOUtils::nodata), cellsize_y(IOUtils::nodata), factor_x(IOUtils::nodata), factor_y(IOUtils::nodata),
indexed(false), meteo_initialized(false), update_dem(false)
cellsize(IOUtils::nodata), factor_x(IOUtils::nodata), factor_y(IOUtils::nodata),
indexed(false), meteo_initialized(false), llcorner_initialized(false), update_dem(false)
{
setOptions();
}
......@@ -111,8 +111,8 @@ GRIBIO::GRIBIO(const Config& cfgreader)
meteo_ext(default_ext), grid2d_ext(default_ext), grid2d_prefix(), idx_filename(), coordin(), coordinparam(),
VW(), DW(), wind_date(), llcorner(), fp(NULL), idx(NULL),
latitudeOfNorthernPole(IOUtils::nodata), longitudeOfNorthernPole(IOUtils::nodata), bearing_offset(IOUtils::nodata),
cellsize_x(IOUtils::nodata), cellsize_y(IOUtils::nodata), factor_x(IOUtils::nodata), factor_y(IOUtils::nodata),
indexed(false), meteo_initialized(false), update_dem(false)
cellsize(IOUtils::nodata), factor_x(IOUtils::nodata), factor_y(IOUtils::nodata),
indexed(false), meteo_initialized(false), llcorner_initialized(false), update_dem(false)
{
setOptions();
}
......@@ -138,12 +138,12 @@ GRIBIO& GRIBIO::operator=(const GRIBIO& source) {
latitudeOfNorthernPole = source.latitudeOfNorthernPole;
longitudeOfNorthernPole = source.longitudeOfNorthernPole;
bearing_offset = source.bearing_offset;
cellsize_x = source.cellsize_x;
cellsize_y = source.cellsize_y;
cellsize = source.cellsize;
factor_x = source.factor_x;
factor_y = source.factor_y;
indexed = source.indexed;
meteo_initialized = source.meteo_initialized;
llcorner_initialized = source.llcorner_initialized;
update_dem = source.update_dem;
}
return *this;
......@@ -292,8 +292,8 @@ void GRIBIO::getDate(grib_handle* h, Date &base, double &d1, double &d2) {
throw InvalidFormatException(ss.str(), AT);
}
d1 = static_cast<double>(startStep*step_units);
d2 = static_cast<double>(endStep*step_units);
d1 = static_cast<double>(startStep)*step_units;
d2 = static_cast<double>(endStep)*step_units;
}
Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y)
......@@ -328,7 +328,7 @@ Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y
GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);
double bearing;
cell_x = Coords::VincentyDistance(cntr_lat, ll_lon, cntr_lat, ur_lon, bearing) / (double)Ni;
cell_y = Coords::VincentyDistance(cntr_lon, ll_lat, cntr_lon, ur_lat, bearing) / (double)Nj;
cell_y = Coords::VincentyDistance(ll_lat, cntr_lon, ur_lat, cntr_lon, bearing) / (double)Nj;
//determining bearing offset
double delta_lat, delta_lon; //geographic coordinates
......@@ -336,16 +336,14 @@ Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y
Coords::VincentyDistance(cntr_lat, cntr_lon, delta_lat, delta_lon, bearing_offset);
bearing_offset = fmod( bearing_offset + 180., 360.) - 180.; // turn into [-180;180)
//computing lower left corner by using the center point as reference
//returning the center point as reference
Coords cntr(coordin, coordinparam);
cntr.setLatLon(cntr_lat, cntr_lon, IOUtils::nodata);
const double cellsize=.5*(cell_x+cell_y);
cntr.moveByXY(-.5*(double)Ni*cellsize, -.5*(double)Nj*cellsize);
return cntr; //this is now the ll corner
return cntr;
}
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& read_geolocalization)
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out)
{
long Ni, Nj;
GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
......@@ -363,18 +361,36 @@ void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& rea
GRIB_CHECK(grib_get_double_array(h,"values",values,&values_len),0);
if(read_geolocalization) llcorner = getGeolocalization(h, cellsize_x, cellsize_y);
grid_out.set(static_cast<size_t>(Ni), static_cast<size_t>(Nj), .5*(cellsize_x+cellsize_y), llcorner);
if (!llcorner_initialized) {
//most important: get cellsize. llcorner will be finalized AFTER aspect ration correction
double cellsize_x, cellsize_y;
llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this is the center cell
cellsize = (double)Optim::round( std::min(cellsize_x, cellsize_y) * 100. ) / 100.; // round to 1cm precision for numerical stability
if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
factor_x = cellsize_x / cellsize;
factor_y = cellsize_y / cellsize;
}
}
grid_out.set(static_cast<size_t>(Ni), static_cast<size_t>(Nj), cellsize, llcorner);
size_t i=0;
for(size_t jj=0; jj<(unsigned)Nj; jj++) {
for(size_t ii=0; ii<(unsigned)Ni; ii++)
grid_out(ii,jj) = values[i++];
}
free(values);
//cells were not square, we have to resample
if (factor_x!=IOUtils::nodata && factor_y!=IOUtils::nodata) {
grid_out.grid2D = ResamplingAlgorithms2D::BilinearResampling(grid_out.grid2D, factor_x, factor_y);
}
free(values);
if (!llcorner_initialized) { //take into account aspect ration conversion for computing true llcorner
llcorner.moveByXY(-.5*(double)grid_out.getNx()*cellsize, -.5*(double)grid_out.getNy()*cellsize );
llcorner_initialized = true;
grid_out.llcorner = llcorner;
}
}
bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out)
......@@ -401,13 +417,12 @@ bool GRIBIO::read2DGrid_indexed(const double& in_marsParam, const long& i_levelT
long level=0;
if(i_level!=0) GRIB_CHECK(grib_get_long(h,"level", &level),0);
if(level==i_level) {
//geolocalization has been initialized when indexing the file, so we don't need to redo it
if( (i_date.isUndef()) ||
(timeRange==0 && i_date==base_date+P1) ||
(timeRange==1 && i_date==base_date) ||
((timeRange==2 || timeRange==3) && i_date>=base_date+P1 && i_date<=base_date+P2) ||
((timeRange==4 || timeRange==5) && i_date==base_date+P2) ) {
read2Dlevel(h, grid_out, false);
read2Dlevel(h, grid_out);
return true;
}
}
......@@ -434,7 +449,7 @@ void GRIBIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
}
read2Dlevel(h, grid_out, true);
read2Dlevel(h, grid_out);
grib_handle_delete(h);
} else {
cleanup();
......@@ -472,14 +487,6 @@ void GRIBIO::indexFile(const std::string& filename)
cleanup();
throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
}
llcorner = getGeolocalization(h, cellsize_x, cellsize_y); //this sets llcorner, cellsize and bearing_offset
if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
// round to 1cm precision for numerical stability
const double cellsize = Optim::round( std::min(cellsize_x, cellsize_y) * 100. ) / 100.;
factor_x = cellsize_x / cellsize;
factor_y = cellsize_y / cellsize;
}
grib_handle_delete(h);
} else {
cleanup();
......@@ -657,7 +664,7 @@ void GRIBIO::read2DGrid(const std::string& filename, Grid2DObject& grid_out, con
void GRIBIO::readDEM(DEMObject& dem_out)
{
const Date d(2012,9,19,0,0,0); //ie: undef. This will be caught when reading the GRIB file
const Date d; //ie: undef. This will be caught when reading the GRIB file
std::string filename;
cfg.getValue("DEMFILE", "Input", filename);
......
meteoio/plugins/GRIBIO.h
View file @ d7fb729c
......@@ -66,7 +66,7 @@ class GRIBIO : public IOInterface {
void listFields(const std::string& filename);
void getDate(grib_handle* h, Date &base, double &d1, double &d2);
Coords getGeolocalization(grib_handle* h, double &cellsize_x, double &cellsize_y);
void read2Dlevel(grib_handle* h, Grid2DObject& grid_out, const bool& read_geolocalization);
void read2Dlevel(grib_handle* h, Grid2DObject& grid_out);
bool read2DGrid_indexed(const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out);
void read2DGrid(const std::string& filename, Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date);
void readWind(const std::string& filename, const Date& date);
......@@ -74,8 +74,6 @@ class GRIBIO : public IOInterface {
void readStations(std::vector<Coords> &vecPoints);
void listKeys(grib_handle** h, const std::string& filename);
void scanMeteoPath();
/*void rotatedToTrueLatLon(const double& lat_rot, const double& lon_rot, double &lat_true, double &lon_true) const;
void trueLatLonToRotated(const double& lat_true, const double& lon_true, double &lat_rot, double &lon_rot) const;*/
void cleanup() throw();
bool removeDuplicatePoints(std::vector<Coords>& vecPoints, double *lats, double *lons);
......@@ -102,13 +100,14 @@ class GRIBIO : public IOInterface {
grib_index *idx; //because it needs to be kept between calls
double latitudeOfNorthernPole, longitudeOfNorthernPole; //for rotated coordinates
double bearing_offset; //to correct vectors coming from rotated lat/lon, we will add an offset to the bearing
double cellsize_x, cellsize_y, factor_x, factor_y;
double cellsize, factor_x, factor_y;
static const std::string default_ext;
static const double plugin_nodata; //plugin specific nodata value, e.g. -999
static const double tz_in; //GRIB time zone
bool indexed; //flag to know if the file has already been indexed
bool meteo_initialized; //set to true after we scanned METEOPATH, filed the cache, read the virtual stations from io.ini
bool llcorner_initialized; //set to true after we properly computed llcorner
bool update_dem;
};
......
meteoio/plugins/NetCDFIO.cc
View file @ d7fb729c
......@@ -161,13 +161,12 @@ void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const std::string& arguments)
}
}
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
void NetCDFIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& /*parameter*/, const Date& date)
{
const string filename = cfg.get("GRID2DFILE", "Input");
string varname="";
read2DGrid_internal(grid_out, filename, varname, date);
}
......
meteoio/plugins/libncpp.cc
View file @ d7fb729c
......@@ -520,7 +520,7 @@ void calculate_dimensions(const mio::Grid2DObject& grid, double*& lat_array, dou
// The idea is to use the difference in coordinates of the upper right and the lower left
// corner to calculate the lat/lon intervals between cells
Coords urcorner(grid.llcorner);
urcorner.setGridIndex(grid.getNx() - 1, grid.getNy() - 1, IOUtils::nodata, true);
urcorner.setGridIndex(static_cast<int>(grid.getNx() - 1), static_cast<int>(grid.getNy() - 1), IOUtils::nodata, true);
grid.gridify(urcorner);
const double lat_interval = (urcorner.getLat() - lat_array[0]) / static_cast<double>(grid.getNy()-1);
......
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