WSL/SLF GitLab Repository

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

Two new methods have been implemented as static of Coords, to compute the...

Two new methods have been implemented as static of Coords, to compute the ground length of one degree of latitude or longitude at any given latitude.

A check for 0x0 sized grids has been added to PNGIO. The ability to specify the number of colors (ie: discretization) in the io.ini has been added.

A breakthrough occured with GRIBIO, that starts to be able to do some things. This is still not usable, but definitely moving forward!
parent 9e0d9c67
......@@ -842,6 +842,40 @@ std::string Coords::decimal_to_dms(const double& decimal) {
return dms.str();
}
/**
* @brief Lenght of one degree of latitude
* This returns the lenght in meters of one degree of latitude around the given latitude
* (ie: latitude-.5, latitude+.5). See https://en.wikipedia.org/wiki/Latitude#The_length_of_a_degree_of_latitude
* @param[in] latitude latitude where to perform the computation
* @return lenght of one degree of latitude
*/
double Coords::lat_degree_lenght(const double& latitude) {
const double to_rad = Cst::PI / 180.0;
const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis
const double b = ellipsoids[E_WGS84].b; //minor ellipsoid semi-axis
const double e2 = (a*a-b*b) / (a*a); //ellispoid eccentricity, squared
const double degree_length = (Cst::PI*a*(1.-e2)) / ( 180.*pow(1.-e2*IOUtils::pow2(sin(latitude*to_rad)), 1.5) );
return fabs( degree_length );
}
/**
* @brief Lenght of one degree of longitude
* This returns the lenght in meters of one degree of longitude around the given latitude
* (ie: latitude-.5, latitude+.5). See https://en.wikipedia.org/wiki/Latitude#The_length_of_a_degree_of_latitude
* @param[in] latitude latitude where to perform the computation
* @return lenght of one degree of longitude
*/
double Coords::lon_degree_lenght(const double& latitude) {
const double to_rad = Cst::PI / 180.0;
const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis
const double b = ellipsoids[E_WGS84].b; //minor ellipsoid semi-axis
const double e2 = (a*a-b*b) / (a*a); //ellispoid eccentricity, squared
const double degree_length = (Cst::PI*a*cos(latitude*to_rad)) / ( 180.*pow(1.-e2*IOUtils::pow2(sin(latitude*to_rad)), .5) );
return fabs( degree_length );
}
/**
* @brief Coordinate conversion: from WGS84 Lat/Long to Swiss grid
* See http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf for more.
......
......@@ -111,6 +111,9 @@ class Coords {
static std::string decimal_to_dms(const double& decimal);
static void parseLatLon(const std::string& coordinates, double& lat, double& lon);
static double lat_degree_lenght(const double& latitude);
static double lon_degree_lenght(const double& latitude);
private:
//Coordinates conversions
void convert_to_WGS84(double easting, double northing, double& latitude, double& longitude) const;
......
......@@ -17,11 +17,11 @@
*/
#include "GRIBIO.h"
#include <meteoio/ResamplingAlgorithms2D.h>
#include <errno.h>
#include <grib_api.h>
//#include "GRIBIO_names.cc" //HACK
using namespace std;
namespace mio {
......@@ -97,7 +97,7 @@ void GRIBIO::listKeys(grib_handle** h, const std::string& filename)
grib_keys_iterator_delete(kiter);
}
void GRIBIO::listContent(const std::string& filename)
void GRIBIO::listContent_old(const std::string& filename)
{
grib_handle* h=NULL;
int err=0, grib_count=0;
......@@ -129,26 +129,13 @@ void GRIBIO::listContent(const std::string& filename)
if(table_id==201 && param_id==133) {
std::cout << table_id << ":" << param_id << " -> Date=" << date.toString(Date::ISO) << " ";
/*int Ni_missing=1;
grib_is_missing(h, "Ni", &Ni_missing);
if(Ni_missing==1) std::cout << "Key Ni missing\n";
else {*/
//grib_get_native_type(h,options->print_keys[i].name,&(options->print_keys[i].type));
/*char Ni[100]="";
size_t len=100;
grib_get_string(h, "Ni", Ni, &len);*/
/*long Ni=-1;
grib_get_long(h,"numberOfPointsAlongAParallel",&Ni);
std::cout << Ni << " Ni ";*/
//}
/*long Ni, Nj;
long Ni, Nj;
GRIB_CHECK(grib_get_long(h,"Nx",&Ni),0);
GRIB_CHECK(grib_get_long(h,"numberOfPointsAlongAMeridian",&Nj),0);
std::cout << Ni << "x" << Nj << "\n";*/
/*long nb_pts;
std::cout << Ni << "x" << Nj << "\n";
long nb_pts;
GRIB_CHECK(grib_get_long(h,"numberOfDataPoints",&nb_pts),0);
std::cout << "nb_pts=" << nb_pts << " ";*/
std::cout << "nb_pts=" << nb_pts << " ";
std::cout << "\n";
/*double *values;
......@@ -161,51 +148,172 @@ void GRIBIO::listContent(const std::string& filename)
printf("%d %g\n",i+1,values[i]);
free(values);*/
listKeys(&h, filename);
//listKeys(&h, filename);
} else {
//std::cout << table_id << ":" << param_id << " -> Date=" << date.toString(Date::ISO) << "\n";
}
}
}
//printf("table2Version=%ld parameter=%ld\n",table_id, param_id);
/*std::cout << table_id << ":" << param_id << " -> ";
if(table_id==2)
std::cout << parm_table_dwd_002[param_id].name << " (" << parm_table_dwd_002[param_id].comment << ")";
else if(table_id==201)
std::cout << parm_table_dwd_201[param_id].name << " (" << parm_table_dwd_201[param_id].comment << ")";
else if(table_id==202)
std::cout << parm_table_dwd_202[param_id].name << " (" << parm_table_dwd_202[param_id].comment << ")";
else if(table_id==203)
std::cout << parm_table_dwd_203[param_id].name << " (" << parm_table_dwd_203[param_id].comment << ")";
else if(table_id==204)
std::cout << parm_table_dwd_204[param_id].name << " (" << parm_table_dwd_204[param_id].comment << ")";
else if(table_id==205)
std::cout << parm_table_dwd_205[param_id].name << " (" << parm_table_dwd_205[param_id].comment << ")";
else
std::cout << "unknown: " << param_id << "\n";
std::cout << "\n";*/
//GRIB_CHECK(grib_get_size(h,"values",&values_len),0); //gridded data in "values" key
Date GRIBIO::getDate(grib_handle* h) {
long dataDate, dataTime;
GRIB_CHECK(grib_get_long(h,"dataDate",&dataDate),0);
GRIB_CHECK(grib_get_long(h,"dataTime",&dataTime),0);
const int year=dataDate/10000, month=dataDate/100-year*100, day=dataDate-month*100-year*10000;
const int hour=dataTime/100, minutes=dataTime-hour*100;
return Date(year, month, day, hour, minutes, tz_in);
}
Coords GRIBIO::getGeolocalization(grib_handle* h, double &cellsize_x, double &cellsize_y)
{
double angleOfRotationInDegrees;
GRIB_CHECK(grib_get_double(h,"angleOfRotationInDegrees",&angleOfRotationInDegrees),0);
if(angleOfRotationInDegrees!=0.) {
throw InvalidArgumentException("Rotated grids not supported!", AT);
}
double ll_latitude;
double ll_longitude;
GRIB_CHECK(grib_get_double(h,"latitudeOfFirstGridPointInDegrees",&ll_latitude),0);
GRIB_CHECK(grib_get_double(h,"longitudeOfFirstGridPointInDegrees",&ll_longitude),0);
double latitudeOfSouthernPole;
double longitudeOfSouthernPole;
GRIB_CHECK(grib_get_double(h,"latitudeOfSouthernPoleInDegrees",&latitudeOfSouthernPole),0);
GRIB_CHECK(grib_get_double(h,"longitudeOfSouthernPoleInDegrees",&longitudeOfSouthernPole),0);
Coords llcorner(coordin, coordinparam);
llcorner.setLatLon( ll_latitude+(90.+latitudeOfSouthernPole), ll_longitude+longitudeOfSouthernPole, 0.);
double ur_latitude;
double ur_longitude;
GRIB_CHECK(grib_get_double(h,"latitudeOfLastGridPointInDegrees",&ur_latitude),0);
GRIB_CHECK(grib_get_double(h,"longitudeOfLastGridPointInDegrees",&ur_longitude),0);
const double cntr_latitude = .5*(ll_latitude+ur_latitude)+(90.+latitudeOfSouthernPole);
double d_i, d_j;
GRIB_CHECK(grib_get_double(h,"jDirectionIncrementInDegrees",&d_j),0);
GRIB_CHECK(grib_get_double(h,"iDirectionIncrementInDegrees",&d_i),0);
cellsize_x = Coords::lon_degree_lenght(cntr_latitude)*d_i;
cellsize_y = Coords::lat_degree_lenght(cntr_latitude)*d_j;
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.) {
stringstream 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 llcorner;
}
void GRIBIO::read2DGrid(Grid2DObject& /*grid_out*/, const std::string& i_name)
void GRIBIO::read2Dlevel(grib_handle* h, Grid2DObject& grid_out)
{
const std::string filename = grid2dpath_in+"/"+i_name;
fp = fopen(filename.c_str(),"r");
if(fp==NULL) {
long Ni, Nj;
GRIB_CHECK(grib_get_long(h,"Ni",&Ni),0);
GRIB_CHECK(grib_get_long(h,"Nj",&Nj),0);
double *values;
size_t values_len= 0;
GRIB_CHECK(grib_get_size(h,"values",&values_len),0);
if(values_len!=(unsigned)(Ni*Nj)) {
stringstream ss;
ss << "Error openning file \"" << filename << "\", possible reason: " << strerror(errno);
throw FileAccessException(ss.str(), AT);
ss << "Declaring grid of size " << Ni << "x" << Nj << "=" << Ni*Nj << " ";
ss << "but containing " << values_len << " values. This is inconsistent!";
throw InvalidArgumentException(ss.str(), AT);
}
values = (double*)malloc(values_len*sizeof(double));
GRIB_CHECK(grib_get_double_array(h,"values",values,&values_len),0);
double cellsize_x, cellsize_y;
const Coords llcorner = getGeolocalization(h, cellsize_x, cellsize_y);
//Grid2DObject tmp_grid(static_cast<unsigned int>(Ni), static_cast<unsigned int>(Nj), cellsize_x, llcorner);
grid_out.set(Ni, Nj, cellsize_x, llcorner); //HACK
int i=0;
for(unsigned int jj=0; jj<(unsigned)Nj; jj++) {
for(unsigned int ii=0; ii<(unsigned)Ni; ii++)
//tmp_grid(ii,jj) = values[i++];
grid_out(ii,jj) = values[i++];
}
free(values);
//grid_out = ResamplingAlgorithms2D::BilinearResampling(tmp_grid, 1., cellsize_x/cellsize_y);
}
listContent(filename);
cleanup();
void GRIBIO::read2DGrid_intern(const std::string& filename, const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out)
{
grib_handle* h=NULL;
int err=0;
//loop over the messages
while((h = grib_handle_new_from_file(0,fp,&err)) != NULL) {
if(!h) {
cleanup();
throw IOException("Unable to create grib handle for \""+filename+"\"", AT);
}
double marsParam;
GRIB_CHECK(grib_get_double(h,"marsParam",&marsParam),0);
if(marsParam==in_marsParam) {
const Date date = getDate(h);
long levelType;
GRIB_CHECK(grib_get_long(h,"indicatorOfTypeOfLevel", &levelType),0); //sfc (surface), pl (pressure level), ml (model level)
long level=0;
if(levelType==105) GRIB_CHECK(grib_get_long(h,"level", &level),0);
if(levelType==i_levelType && level==i_level && date==i_date) {
read2Dlevel(h, grid_out);
return;
}
}
}
}
void GRIBIO::read2DGrid(Grid2DObject& /*grid_out*/, const std::string& /*i_name*/)
{
//Nothing so far
throw IOException("Nothing implemented here", AT);
}
void GRIBIO::read2DGrid(Grid2DObject& /*grid_out*/, const MeteoGrids::Parameters& /*parameter*/, const Date& /*date*/)
void GRIBIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
{
//Nothing so far
throw IOException("Nothing implemented here", AT);
//in 115.260 and 116.260, swiss coordinates for each cell
const std::string prefix="laf";
const std::string ext=".grb";
const std::string filename = grid2dpath_in+"/"+prefix+date.toString(Date::NUM).substr(0,10)+"f"+date.toString(Date::NUM).substr(10,2)+ext;
fp = fopen(filename.c_str(),"r");
if(fp==NULL) {
stringstream ss;
ss << "Error openning file \"" << filename << "\", possible reason: " << strerror(errno);
throw FileAccessException(ss.str(), AT);
}
if(parameter==MeteoGrids::DW) read2DGrid_intern(filename, 31.2, 105, 10, date, grid_out); //10m wind direction, level 10, type 105
if(parameter==MeteoGrids::VW) read2DGrid_intern(filename, 32.2, 105, 10, date, grid_out); //10m wind speed, level 10, type 105
if(parameter==MeteoGrids::TA) read2DGrid_intern(filename, 17.2, 105, 2, date, grid_out); //2m TA, type 105, level 2
if(parameter==MeteoGrids::RH) read2DGrid_intern(filename, 52.2, 105, 2, date, grid_out); //type 105, level 2
if(parameter==MeteoGrids::TSS) read2DGrid_intern(filename, 11.2, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::HNW) read2DGrid_intern(filename, 61.2, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::ILWR) read2DGrid_intern(filename, 25.201, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::ISWR) read2DGrid_intern(filename, 111.250, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::P) read2DGrid_intern(filename, 1.2, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::HS) read2DGrid_intern(filename, 66.2, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::DEM) read2DGrid_intern(filename, 8.2, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::SLOPE) read2DGrid_intern(filename, 98.202, 1, 0, date, grid_out); //type 1
if(parameter==MeteoGrids::AZI) read2DGrid_intern(filename, 99.202, 1, 0, date, grid_out); //type 1
cleanup();
}
void GRIBIO::readDEM(DEMObject& /*dem_out*/)
......
......@@ -63,7 +63,11 @@ class GRIBIO : public IOInterface {
private:
void setOptions();
void listContent(const std::string& filename);
void listContent_old(const std::string& filename); //HACK: to delete when done
Date getDate(grib_handle* h);
Coords getGeolocalization(grib_handle* h, double &cellsize_x, double &cellsize_y);
void read2Dlevel(grib_handle* h, Grid2DObject& grid_out);
void read2DGrid_intern(const std::string& filename, const double& in_marsParam, const long& i_levelType, const long& i_level, const Date i_date, Grid2DObject& grid_out);
void listKeys(grib_handle** h, const std::string& filename);
void cleanup() throw();
......
......@@ -33,6 +33,7 @@ namespace mio {
* No data read has been implemented, because reading an existing file would require the exact knowlege of the color gradient that has been used
* to create it. When writing grids, various color gradients will be used depending on the parameter that the data represents. Nodata values
* are represented by transparent pixels (transparency is acheived through a transparent color instead of a true alpha channel for size and performance).
* If a grid containing no data (ie: size 0x0) is sent to the plugin, then no file will be written.
* Finally, the naming scheme for meteo grids should be: YYYY-MM-DDTHH.mm_{MeteoGrids::Parameters}.png
*
* @section template_units Units
......@@ -49,7 +50,10 @@ namespace mio {
* - PNG_SCALING: scaling algorithm, either nearest or bilinear (default=bilinear)
* - PNG_AUTOSCALE: autoscale for the color gradient? (default=true)
* - PNG_WORLD_FILE: create world file with each file? (default=false)
*
* Advanced parameters (ie: don't mess up with them if you don't know what you're doing):
* - PNG_INDEXED: create an indexed PNG? (default=true)
* - PNG_NR_LEVELS: number of colors to use (less=smaller files, but it must be at least 5 and less than 255. default=30)
* - PNG_SPEED_OPTIMIZE: optimize file creation for speed? (default=true, otherwise optimize for file size)
*
* The size are specified as width followed by height, with the separator being either a space, 'x' or '*'. If a minimum and a maximum size are given, the average of the smallest and largest permissible sizes will be used.
......@@ -92,7 +96,6 @@ const double PNGIO::plugin_nodata = -999.; //plugin specific nodata value. It ca
const unsigned char PNGIO::channel_depth = 8;
const unsigned char PNGIO::channel_max_color = 255;
const unsigned char PNGIO::transparent_grey = channel_max_color;
const unsigned char PNGIO::nr_levels = 30;
PNGIO::PNGIO(void (*delObj)(void*), const Config& i_cfg) : IOInterface(delObj), cfg(i_cfg)
{
......@@ -142,6 +145,13 @@ void PNGIO::setOptions()
cfg.getValue("PNG_INDEXED", "Output", indexed_png, Config::nothrow);
optimize_for_speed = true;
cfg.getValue("PNG_SPEED_OPTIMIZE", "Output", optimize_for_speed, Config::nothrow);
nr_levels = 30;
unsigned int tmp=0;
cfg.getValue("PNG_NR_LEVELS", "Output", tmp, Config::nothrow);
if(tmp>255 || tmp<5) {
throw InvalidFormatException("PNG_NR_LEVELS must be between 5 and 255!", AT);
}
if(tmp>0) nr_levels=static_cast<unsigned char>(tmp);
}
void PNGIO::parse_size(const std::string& size_spec, unsigned int& width, unsigned int& height)
......@@ -455,6 +465,8 @@ void PNGIO::write2DGrid(const Grid2DObject& grid_in, const std::string& filename
//scale input image
const Grid2DObject grid = scaleGrid(grid_in);
const unsigned int ncols = grid.ncols, nrows = grid.nrows;
if(ncols==0 || nrows==0) return;
const double min = grid.grid2D.getMin();
const double max = grid.grid2D.getMax();
......@@ -498,6 +510,8 @@ void PNGIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameter
//scale input image
Grid2DObject grid = scaleGrid(grid_in);
const unsigned int ncols = grid.ncols, nrows = grid.nrows;
if(ncols==0 || nrows==0) return;
double min = grid.grid2D.getMin();
double max = grid.grid2D.getMax();
......
......@@ -84,6 +84,7 @@ class PNGIO : public IOInterface {
bool has_world_file; ///< create world file with each file?
bool optimize_for_speed; ///< optimize for speed instead of compression?
bool indexed_png; ///< write an indexed png?
unsigned char nr_levels; ///< number of levels to represent? (less-> smaller file size and faster)
std::string coordout, coordoutparam; //projection parameters
std::string grid2dpath;
......@@ -96,7 +97,6 @@ class PNGIO : public IOInterface {
static const unsigned char channel_depth;
static const unsigned char channel_max_color;
static const unsigned char transparent_grey;
static const unsigned char nr_levels; ///< number of levels to represent? (less-> smaller file size and faster)
};
} //namespace
......
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