WSL/SLF GitLab Repository

Commit 69354526 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A bug has been fixed when doing nearest neighbour 2D interpolation. The GRIB...

A bug has been fixed when doing nearest neighbour 2D interpolation. The GRIB plugin can now handle non-square cells (by resampling the domain). The cellsize and geolocalization in NetCDF are now handled the same way as in GRIB: the error is minimized in the center of the domain and the llcorner is back-calculated. Therefore very large domains are handled better. The first elements for supporting ECMWF grids have been implemented, the DEMs can now be read (but their remain as geopotential instead of height).
parent bdc7c9e1
......@@ -31,7 +31,6 @@
namespace mio {
/**
* @class Config
* @brief A class that reads a key/value file. These files (typically named *.ini) follow the INI file format standard (see http://en.wikipedia.org/wiki/INI_file) and have the following structure:
......
......@@ -75,6 +75,7 @@ void legend::simpleLegend(const unsigned int &height, const double &minimum, con
grid.resize(total_width, height, IOUtils::nodata);
const double level_inc = (maximum-minimum)/(double)(nb_labels-1); //the infamous interval thing...
//vertically center legend in free space
if(height>=total_height) {
const unsigned int free_space = height-total_height;
const unsigned int start_legend = free_space/2; //we will start from the bottom
......@@ -124,6 +125,7 @@ void legend::smartLegend(const unsigned int &height, const double &minimum, cons
nb_labels_norm=1;
}
//vertically center legend in free space
const unsigned int smart_height = nb_labels_norm*legend::label_height+legend::interline;
if(height>=smart_height) {
const unsigned int free_space = height-smart_height;
......@@ -382,7 +384,7 @@ void Gradient::getColor(const double& val, unsigned char& r, unsigned char& g, u
a=false;
double r_d,g_d,b_d;
double val_norm;
if(autoscale && val<min) val_norm=0.;
if(autoscale && val<min) val_norm=0.;
else if(autoscale && val>max) val_norm=1.;
else val_norm = (val-min)/delta;
......@@ -636,14 +638,16 @@ gr_terrain::gr_terrain(const double& /*i_min*/, const double& i_max, const bool&
X.push_back(1.); v_h.push_back(22.); v_s.push_back(.2); v_v.push_back(.5); //light maroon
} else {
const double snow_line = i_max - 500.;
X.push_back(0.); v_h.push_back(198.); v_s.push_back(.50); v_v.push_back(.74); //sea, light blue
X.push_back(0.); v_h.push_back(144.); v_s.push_back(.50); v_v.push_back(.39); //sea level, dark green
X.push_back(.4); v_h.push_back(46.); v_s.push_back(.54); v_v.push_back(.86); //yellow, 1200m
X.push_back(.73333); v_h.push_back(4.); v_s.push_back(.71); v_v.push_back(.53); //dark red, 2200m
X.push_back(.9); v_h.push_back(22.); v_s.push_back(.88); v_v.push_back(.41); //maroon, 2700m
X.push_back(.98333); v_h.push_back(22.); v_s.push_back(.36); v_v.push_back(.79); //light maroon, 2950m
X.push_back(1.); v_h.push_back(0.); v_s.push_back(0.); v_v.push_back(.7); //light gray == permanent snow line, 3000m
for(size_t i=0; i<X.size(); i++) X[i] = X[i]*snow_line/i_max; //snow line is the reference
for(size_t ii=2; ii<X.size(); ii++) X[ii] = X[ii]*snow_line/i_max; //snow line is the reference
X.push_back(1.); v_h.push_back(0.); v_s.push_back(0.); v_v.push_back(.95); //almost white == fully glaciated line
}
......
......@@ -25,10 +25,7 @@ using namespace std;
namespace mio {
/**
* @brief Bilinear spatial data resampling
*/
const Grid2DObject ResamplingAlgorithms2D::BilinearResampling(const Grid2DObject &i_grid, const double &factor)
const Grid2DObject ResamplingAlgorithms2D::NearestNeighbour(const Grid2DObject &i_grid, const double &factor)
{
if (factor<=0) {
ostringstream ss;
......@@ -40,11 +37,14 @@ const Grid2DObject ResamplingAlgorithms2D::BilinearResampling(const Grid2DObject
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);
Bilinear(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
NearestNeighbour(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
return o_grid;
}
const Grid2DObject ResamplingAlgorithms2D::cubicBSplineResampling(const Grid2DObject &i_grid, const double &factor)
/**
* @brief Bilinear spatial data resampling
*/
const Grid2DObject ResamplingAlgorithms2D::BilinearResampling(const Grid2DObject &i_grid, const double &factor)
{
if (factor<=0) {
ostringstream ss;
......@@ -56,11 +56,11 @@ const Grid2DObject ResamplingAlgorithms2D::cubicBSplineResampling(const Grid2DOb
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);
cubicBSpline(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
Bilinear(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
return o_grid;
}
const Grid2DObject ResamplingAlgorithms2D::NearestNeighbour(const Grid2DObject &i_grid, const double &factor)
const Grid2DObject ResamplingAlgorithms2D::cubicBSplineResampling(const Grid2DObject &i_grid, const double &factor)
{
if (factor<=0) {
ostringstream ss;
......@@ -72,7 +72,7 @@ const Grid2DObject ResamplingAlgorithms2D::NearestNeighbour(const Grid2DObject &
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);
NearestNeighbour(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
cubicBSpline(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
return o_grid;
}
......@@ -128,12 +128,10 @@ void ResamplingAlgorithms2D::NearestNeighbour(Array2D<double> &o_grid, const Arr
const double scale_y = (double)dest_ny / (double)org_ny;
for (size_t jj=0; jj<dest_ny; jj++) {
size_t org_jj = (size_t) Optim::round( (double)jj/scale_y );
if(org_jj>=org_ny) org_jj=org_ny-1;
const size_t org_jj = std::min( (size_t) Optim::floor( (double)jj/scale_y ) , org_ny-1 );
for (size_t ii=0; ii<dest_nx; ii++) {
size_t org_ii = (size_t) Optim::round( (double)ii/scale_x );
if(org_ii>=org_nx) org_ii=org_nx-1;
const size_t org_ii = std::min( (size_t) Optim::floor( (double)ii/scale_x ) , org_nx-1 );
o_grid(ii,jj) = i_grid(org_ii, org_jj);
}
}
......
......@@ -342,16 +342,6 @@ Coords GRIBIO::getGeolocalization(grib_handle* h, double &cell_x, double &cell_y
const double cellsize=.5*(cell_x+cell_y);
cntr.moveByXY(-.5*(double)Ni*cellsize, -.5*(double)Nj*cellsize);
//checking that cellsize does not vary too much across the grid
/*const double cellsize_x_ll = Coords::lon_degree_lenght(ll_latitude)*d_j;
const double cellsize_x_ur = Coords::lon_degree_lenght(ur_latitude)*d_j;
if( fabs(cellsize_x_ll-cellsize_x_ur)/cellsize_x > 1./100.) {
ostringstream ss;
ss << "Cell size varying too much in the x direction between lower left and upper right corner: ";
ss << cellsize_x_ll << "m to " << cellsize_x_ur << "m";
throw IOException(ss.str(), AT);
}*/
return cntr; //this is now the ll corner
}
......@@ -373,13 +363,7 @@ 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);
/*if( fabs(cellsize_x-cellsize_y)/cellsize_x > 1./100.) {
free(values);
throw InvalidArgumentException("Unsupported geometry: cells can not be represented by square cells!", AT);
}*/
}
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);
size_t i=0;
for(size_t jj=0; jj<(unsigned)Nj; jj++) {
......
......@@ -228,6 +228,12 @@ void NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& fi
if (dimid.size()!=3 || dimlen[0]<1 || dimlen[lat_index]<2 || dimlen[lon_index]<2)
throw IOException("Variable '" + varname + "' may only have three dimensions, all have to at least have length 1", AT);
} else if (dimid.size()==3 && dimlen[0]==1) { //in case the variable is associated with a 1 element time dimension
lat_index = 1;
lon_index = 2;
if (dimlen[lat_index]<2 || dimlen[lon_index]<2)
throw IOException("All dimensions for variable '" + varname + "' have to at least have length 1", AT);
} else if (dimid.size()!=2 || dimlen[lat_index]<2 || dimlen[lon_index]<2) {
throw IOException("Variable '" + varname + "' may only have two dimensions and both have to have length >1", AT);
}
......@@ -249,27 +255,66 @@ void NetCDFIO::read2DGrid_internal(Grid2DObject& grid_out, const std::string& fi
read_data(ncid, varname, varid, grid);
}
copy_grid(dimlen[lat_index], dimlen[lon_index], lat, lon, grid, grid_out);
double missing_value=plugin_nodata;
if (ncpp::check_attribute(ncid, varid, "missing_value"))
ncpp::get_attribute(ncid, varname, varid, "missing_value", missing_value);
close_file(filename, ncid);
copy_grid(dimlen[lat_index], dimlen[lon_index], lat, lon, grid, missing_value, grid_out);
//handle data packing if necessary
if (ncpp::check_attribute(ncid, varid, "scale_factor")) {
double scale_factor=1.;
ncpp::get_attribute(ncid, varname, varid, "scale_factor", scale_factor);
grid_out.grid2D *= scale_factor;
}
if (ncpp::check_attribute(ncid, varid, "add_offset")) {
double add_offset=0.;
ncpp::get_attribute(ncid, varname, varid, "add_offset", add_offset);
grid_out.grid2D += add_offset;
}
close_file(filename, ncid);
delete[] lat; delete[] lon; delete[] grid;
}
void NetCDFIO::copy_grid(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
const double * const grid, Grid2DObject& grid_out)
const double * const grid, const double& nodata, Grid2DObject& grid_out)
{
Coords location(coordin, coordinparam);
location.setLatLon(lat[0], lon[0], grid[0]);
double resampling_factor_x = IOUtils::nodata, resampling_factor_y=IOUtils::nodata;
const double cellsize = calculate_cellsize(latlen, lonlen, lat, lon, resampling_factor_x, resampling_factor_y);
grid_out.set(lonlen, latlen, cellsize, location);
const double cntr_lat = .5*(lat[0]+lat[latlen-1]);
const double cntr_lon = .5*(lon[0]+lon[lonlen-1]);
for (size_t kk=0; kk < latlen; kk++) {
for (size_t ll=0; ll < lonlen; ll++) {
grid_out(ll, kk) = IOUtils::standardizeNodata(grid[kk*lonlen + ll], plugin_nodata);
//computing lower left corner by using the center point as reference
Coords cntr(coordin, coordinparam);
cntr.setLatLon(cntr_lat, cntr_lon, IOUtils::nodata);
cntr.moveByXY(-.5*(double)(lonlen-1)*cellsize, -.5*(double)(latlen-1)*cellsize);
grid_out.set(lonlen, latlen, cellsize, cntr);
//Handle the case of llcorner/urcorner swapped
if (lat[0]<=lat[latlen-1]) {
for (size_t kk=0; kk < latlen; kk++) {
const size_t row = kk*lonlen;
if (lon[0]<=lon[lonlen-1]) {
for (size_t ll=0; ll < lonlen; ll++)
grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + ll], nodata);
} else {
for (size_t ll=0; ll < lonlen; ll++)
grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + (lonlen -1) - ll], nodata);
}
}
} else {
for (size_t kk=0; kk < latlen; kk++) {
const size_t row = ((latlen-1) - kk)*lonlen;
if (lon[0]<=lon[lonlen-1]) {
for (size_t ll=0; ll < lonlen; ll++)
grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + ll], nodata);
} else {
for (size_t ll=0; ll < lonlen; ll++)
grid_out(ll, kk) = IOUtils::standardizeNodata(grid[row + (lonlen -1) - ll], nodata);
}
}
}
......@@ -290,24 +335,12 @@ void NetCDFIO::copy_grid(const size_t& latlen, const size_t& lonlen, const doubl
double NetCDFIO::calculate_cellsize(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
double& factor_x, double& factor_y)
{
Coords llcorner(coordin, coordinparam);
llcorner.setLatLon(lat[0], lon[0], IOUtils::nodata);
Coords urcorner(coordin, coordinparam);
urcorner.setLatLon(lat[latlen-1], lon[lonlen-1], IOUtils::nodata);
const double ll_easting=llcorner.getEasting(), ll_northing=llcorner.getNorthing();
const double ur_easting=urcorner.getEasting(), ur_northing=urcorner.getNorthing();
const double distanceX = ur_easting - ll_easting;
const double distanceY = ur_northing - ll_northing;
if(distanceX<0 || distanceY<0) {
ostringstream ss;
ss << "Can not compute cellsize: this is most probably due to an inappropriate input coordinate system (COORDSYS).";
ss << "Please configure one that can accomodate (" << llcorner.getLat() << "," << llcorner.getLon() << ") - ";
ss << "(" << urcorner.getLat() << "," << urcorner.getLon() << ")";
throw InvalidArgumentException(ss.str(), AT);
}
const double cntr_lat = .5*(lat[0]+lat[latlen-1]);
const double cntr_lon = .5*(lon[0]+lon[lonlen-1]);
double alpha;
const double distanceX = Coords::VincentyDistance(cntr_lat, lon[0], cntr_lat, lon[lonlen-1], alpha);
const double distanceY = Coords::VincentyDistance(lat[0], cntr_lon, lat[latlen-1], cntr_lon, alpha);
// lonlen, latlen are decremented by 1; n linearly connected points have (n-1) connections
const double cellsize_x = distanceX / (lonlen-1);
......
......@@ -88,7 +88,7 @@ class NetCDFIO : public IOInterface {
void write2DGrid_internal(const Grid2DObject& grid_in, const std::string& filename, const std::string& varname, const Date& date=Date());
void fill_data(const Grid2DObject& grid, double*& data);
void copy_grid(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
const double * const grid, Grid2DObject& grid_out);
const double * const grid, const double& nodata, Grid2DObject& grid_out);
double calculate_cellsize(const size_t& latlen, const size_t& lonlen, const double * const lat, const double * const lon,
double& factor_x, double& factor_y);
void calculate_dimensions(const Grid2DObject& grid, double*& lat_array, double*& lon_array);
......
......@@ -69,13 +69,13 @@ void get_attribute(const int& ncid, const std::string& varname, const int& varid
int status = nc_inq_attlen (ncid, varid, attr_name.c_str(), &attr_len);
if (status != NC_NOERR)
throw IOException("Could not retrieve attribute '" + attr_name + "'for var '" + varname + "': " + nc_strerror(status), AT);
throw IOException("Could not retrieve attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
char* value = new char[attr_len + 1]; // +1 for trailing null
status = nc_get_att_text(ncid, varid, attr_name.c_str(), value);
if (status != NC_NOERR)
throw IOException("Could not read attribute '" + attr_name + "'for var '" + varname + "': " + nc_strerror(status), AT);
throw IOException("Could not read attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
value[attr_len] = '\0';
attr_value = string(value);
......@@ -83,6 +83,19 @@ void get_attribute(const int& ncid, const std::string& varname, const int& varid
delete[] value;
}
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, double& attr_value)
{
size_t attr_len;
int status = nc_inq_attlen (ncid, varid, attr_name.c_str(), &attr_len);
if (status != NC_NOERR)
throw IOException("Could not retrieve attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
status = nc_get_att_double(ncid, varid, attr_name.c_str(), &attr_value);
if (status != NC_NOERR)
throw IOException("Could not read attribute '" + attr_name + "' for var '" + varname + "': " + nc_strerror(status), AT);
}
bool check_attribute(const int& ncid, const int& varid, const std::string& attr_name)
{
size_t attr_len;
......
......@@ -45,6 +45,7 @@ namespace ncpp {
void add_attribute(const int& ncid, const int& varid, const std::string& attr_name, const std::string& attr_value);
void add_attribute(const int& ncid, const int& varid, const std::string& attr_name, const double& attr_value);
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, std::string& attr_value);
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, double& attr_value);
bool check_attribute(const int& ncid, const int& varid, const std::string& attr_name);
//Adding dimensions
......
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