/***************************************************************************
* Copyright (C) 2005 by Nicolas PASCAL *
* nicolas.pascal@icare.univ-lille1.fr *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "parasolfiledata.h"
const double PARASOLFileData::lat2time_coef=3840./233.;
const int PARASOLFileData::nb_parasol_direction = 16;
const double PARASOLFileData::satellite_altitude = 705.e3;
const double PARASOLFileData::pix_sz_nadir = 6e3; // 6 km
const double PARASOLFileData::pix_obs_nadir_angle = atan2(PARASOLFileData::pix_sz_nadir / 2., PARASOLFileData::satellite_altitude);
const double PARASOLFileData::max_sat_nadir_angle = radians (75.);
const double PARASOLFileData::pix_sz_max = PARASOLFileData::satellite_altitude * (tan (PARASOLFileData::max_sat_nadir_angle) - tan (PARASOLFileData::max_sat_nadir_angle - PARASOLFileData::pix_obs_nadir_angle));
const int PARASOLFileData::nb_acq_seq_max = 130;
const int PARASOLFileData::nb_acq_img_max = 9;
PARASOLFileData::PARASOLFileData(const string& name, const string& mode):
FileData(name,mode), SatelliteFileData(name,mode) {
init();
}
void PARASOLFileData::init() {
instrument=-1;
level=-1;
processing_line=UNDEFINED;
product=' ';
orbit_cycle_nb=-1;
orbit_nb=-1;
reprocessing_nb=' ';
grid_factor=-1;
leader=NULL;
line_data=NULL;
col_data=NULL;
software_version_number = "";
// extract the product infos from the file name
string short_filename = get_tail(get_name());
parse_filename(short_filename);
// init grid paramtres
set_grid_factor();
sz_cell = (double)(grid_factor) / 18.;
sz_grid [0] = 3240 / grid_factor;
sz_grid [1] = 6480 / grid_factor;
// init the leader file reader
string filename_skeleton=get_directory(get_name())+get_radix();
// remove ending "L" or "D" character
if ((*filename_skeleton.rbegin())=='L' || (*filename_skeleton.rbegin())=='D')
filename_skeleton=filename_skeleton.substr(0,filename_skeleton.size()-1);
leader = new PARASOLLeader(filename_skeleton+"L");
data = new PARASOLData(filename_skeleton+"D");
if (leader!=NULL) {
leader->open_file();
//-- read the software version number
leader->read_record(LEADER_FILE_DESCRIPTOR);
software_version_number = leader->leader_file_descriptor->soft_version;
//-- read the spatio-temporal parameters
leader->read_record(SPATIO_TEMPORAL_CHARACTERISTICS);
// init the time coverage
date->set_date_str(leader->spatio_temp_char->first_acq_date,"%Y%m%d%H%M%S");
Date end_date;
end_date.set_date_str(leader->spatio_temp_char->last_acq_date,"%Y%m%d%H%M%S");
time_coverage=end_date.get_TAI93_time()-date->get_TAI93_time();
// init the spatial coverage
// unsigned short line_minmax[2];
// float lat_minmax[2];
// int grid_idx[2];
// grid_idx[1]=0; // dummy col number
// set the min latitude using the min line index
int irow = MyTools::string2int(leader->spatio_temp_char->line_nb_southern_pix);
lat_min = get_grid_lat(irow);
// line_minmax[0]=
// grid_idx[0]=line_minmax[0]; // min line
// grid_to_geolocation(grid_idx, lat_minmax);
// lat_min=lat_minmax[0];
// set the max latitude using the max line index
irow = MyTools::string2int(leader->spatio_temp_char->line_nb_northern_pix);
lat_max = get_grid_lat(irow);
// grid_idx[0]=line_minmax[1]; // max line
// grid_to_geolocation(grid_idx,lat_minmax);
// lat_max=lat_minmax[1];
// init the scaling factors
leader->read_record(LEADER_FILE_DESCRIPTOR);
// leader->print_leader_file_descriptor();
leader->read_record(SCALING_FACTORS);
// leader->print_scaling_factors();
// leader->read_record(SPATIO_TEMPORAL_CHARACTERISTICS);
// leader->print_spatio_temp_char();
// leader->read_record(TECHNOLOGICAL_PARAMETERS);
// leader->print_techno_param();
leader->close_file();
}
// geolocation buffers
line_data = NULL;
col_data = NULL;
// viewing directions buffers
v_alt_view = NULL;
v_saa = NULL;
v_vza = NULL;
v_raa = NULL;
v_ni = NULL;
v_sn = NULL;
v_x_sat = NULL;
v_y_sat = NULL;
v_z_sat = NULL;
}
PARASOLFileData::~PARASOLFileData() {
delete leader, leader = NULL;
delete data, data = NULL;
free_geolocation_data();
}
bool PARASOLFileData::check_filename( const string& short_filename ) const {
string parasol_filename(short_filename);
bool check_instrument=( parasol_filename[0]=='P' && (parasol_filename[1]=='1' ||
parasol_filename[1]=='2' ||
parasol_filename[1]=='3') );
bool check_level=check_instrument && ( parasol_filename[2]=='L' && (parasol_filename[3]=='1' ||
parasol_filename[3]=='2') );
if (parasol_filename[3]=='3')
cerr<<"In "<<__FILE__<<" at "<<__LINE__<<" : Sorry but level 3 PARASOL files can't be read yet. The PARASOLFileData class must be completed"<<endl;
bool check_processing_line=check_level && ( parasol_filename[4]=='T' &&
(parasol_filename[5]=='B' ||
parasol_filename[5]=='L' ||
parasol_filename[5]=='O' ||
parasol_filename[5]=='R') );
bool check_product=check_processing_line && ( parasol_filename[6]=='G' &&
(parasol_filename[7]=='A' ||
parasol_filename[7]=='B' ||
parasol_filename[7]=='C' ||
parasol_filename[7]=='1') ) ;
bool check_orbit=check_product;
if (parasol_filename[3]!='3') { // attributes specific to level 1 and 2 files
string s_orbit_cycle_nb = parasol_filename.substr(8,3);
string s_orbit_nb = parasol_filename.substr(11,3);
string s_reprocessing_nb = parasol_filename.substr(14,1);
char s_reprocessing_nb_char = s_reprocessing_nb[0];
check_orbit=check_orbit && MyTools::is_int(s_orbit_cycle_nb) && MyTools::is_int(s_orbit_nb);
check_orbit=check_orbit && (s_reprocessing_nb_char>='A' && s_reprocessing_nb_char<='Z');
} else {
// TODO level 3 attributes must be implemented here
}
bool all_checked = check_orbit && ( parasol_filename[15]=='L' || parasol_filename[15]=='D' );
if ( !all_checked ) {
invalid_filename e(__FILE__,__LINE__,parasol_filename,"PARASOL");
throw e;
}
return all_checked;
}
void PARASOLFileData::parse_filename( const string& short_filename ) {
check_filename(short_filename); // throw a parse_PARASOL_filename_error if the filename doesn't matches the pattern
try {
string parasol_filename(short_filename);
instrument = MyTools::char2int(parasol_filename[1]);
level=MyTools::char2int(parasol_filename[3]);
char processing_line_code=parasol_filename[5];
if (processing_line_code == 'B')
processing_line=BASIC;
else if (processing_line_code == 'L')
processing_line=LAND_SURFACE;
else if (processing_line_code == 'O')
processing_line=OCEAN_COLOR;
else if (processing_line_code == 'R')
processing_line=RADIATION_BUDGET;
else
processing_line=UNDEFINED;
product=parasol_filename[7];
if (level!=3) { // attributes specific to level 1 and 2 files
string s_orbit_cycle_nb = parasol_filename.substr(8,3);
orbit_cycle_nb=MyTools::string2int(s_orbit_cycle_nb);
string s_orbit_nb = parasol_filename.substr(11,3);
orbit_nb=MyTools::string2int(s_orbit_nb);
reprocessing_nb=parasol_filename[14];
} else {
// TODO level 3 attributes must be implemented here
}
} catch (invalid_filename e) {
cerr<<e.what()<<endl;
throw;
}
}
const string PARASOLFileData::get_processing_line_string( ) const
{
if (processing_line == BASIC)
return "basic";
else if (processing_line == LAND_SURFACE)
return "land surface";
else if (processing_line == OCEAN_COLOR)
return "ocean color";
else if (processing_line == RADIATION_BUDGET)
return "radiation budget";
else
return "undefined";
}
const char PARASOLFileData::get_processing_line_char( ) const
{
if (processing_line == BASIC)
return 'B';
else if (processing_line == LAND_SURFACE)
return 'L';
else if (processing_line == OCEAN_COLOR)
return 'O';
else if (processing_line == RADIATION_BUDGET)
return 'R';
else
return ' ';
}
const string PARASOLFileData::get_radix()
{
string ret="P"+MyTools::to_string(instrument)+
"L"+MyTools::to_string(level)+
"T"+get_processing_line_char()+
"G"+product+
MyTools::int2string(orbit_cycle_nb,3)+MyTools::int2string(orbit_nb,3)+
reprocessing_nb;
return ret;
}
const double PARASOLFileData::get_grid_lat(const int irow) const {
if (irow <= 0 || irow > sz_grid [0]) {
g_exception e(__FILE__,__LINE__, "Invalid row indice " + MyTools::to_string(irow));
throw e;
}
return 90. - ((double)(irow) - .5) * sz_cell;
}
const int PARASOLFileData::get_nrow() const {
return sz_grid [0];
}
const int PARASOLFileData::get_ncol(const int irow) const {
if (irow == -1)
return sz_grid [0];
int ni = (int)(round (sz_grid [0] * sin (radians((irow - 0.5) * sz_cell))));
return 2 * ni;
}
const int PARASOLFileData::get_ncol(const double lat) const {
int ni = (int)(round(sz_grid [0] * cos(radians(lat))));
return 2 * ni;
}
void PARASOLFileData::geolocation_to_grid(const float lat, const float lon, int & irow, int & icol) {
check_geolocation(lat, lon);
irow = (int)(round((1. / sz_cell) * (90. - lat) + 0.5));
int ni = get_ncol (irow) / 2;
icol = (int)(round((sz_grid[0] + 0.5) + (ni / 180.) * lon));
}
void PARASOLFileData::geolocation_to_grid( const float * lat_lon, int * grid_idx ){
geolocation_to_grid (lat_lon[0], lat_lon[1], grid_idx[0], grid_idx[1]);
}
void PARASOLFileData::grid_to_geolocation( const int igrid[], float lat_minmax[], float lon_minmax[]) {
// grid middle pixel lat,lon
float lat, lon;
grid_to_geolocation(igrid, lat, lon);
double dlat = 0.5 * sz_cell; // 2*delta = size of a grid cell side along latitudes
lat_minmax[0] = lat - dlat;
lat_minmax[1] = lat + dlat;
int ni = get_ncol (lat) / 2;
double dlon = 180. / ni;
lon_minmax[0] = lon - dlon;
lon_minmax[1] = lon + dlon;
}
void PARASOLFileData::grid_to_geolocation(const int irow, const int icol, float & lat, float & lon) {
check_grid_coord(irow, icol);
lat = get_grid_lat(irow);
int ni = get_ncol (lat) / 2;
double dlon = 180. / ni;
lon = dlon * (icol - (sz_grid [0] + .5));
}
void PARASOLFileData::grid_to_geolocation(const int igrid[], float & lat, float & lon) {
grid_to_geolocation (igrid[0], igrid[1], lat, lon);
}
void PARASOLFileData::grid_to_geolocation(const int igrid[], float pos[]) {
grid_to_geolocation (igrid[0], igrid[1], pos[0], pos[1]);
}
void PARASOLFileData::set_grid_factor()
{
// algo build on the PARASOL level 1 Standard product manual Appendix A
if (level==1) {
grid_factor=1;
} else if (level == 2) {
if (processing_line==RADIATION_BUDGET) {
grid_factor=3;
} else if (processing_line==OCEAN_COLOR) {
if (product=='A' || product=='B')
grid_factor=1;
else if (product=='C')
grid_factor=3;
} else if (processing_line==LAND_SURFACE) {
if (product=='A')
grid_factor=1;
else if (product=='C')
grid_factor=3;
}
} else if (level == 3) {
if (processing_line==RADIATION_BUDGET) {
grid_factor=3;
} else if (processing_line==OCEAN_COLOR) {
if (product=='B')
grid_factor=1;
else if (product=='C')
grid_factor=3;
} else if (processing_line==LAND_SURFACE) {
if (product=='A' || product=='B')
grid_factor=1;
else if (product=='C')
grid_factor=3;
}
}
}
bool PARASOLFileData::init_read_data_range( const int entry_idx, int * &start, int * &edges, int &rank )
{
assert((start==NULL && edges==NULL) || (start!=NULL && edges!=NULL));
bool allocation_done=false;
/* initialize the data selector range */
// EntryFormat *entry_format=const_cast<EntryFormat *>(get_entry_format(entry_idx));
vector<int> data_dimension=get_data_dimension(entry_idx);
int my_rank=data_dimension.size();
if (rank!=-1 && rank!=my_rank) {
bad_rank e(rank,my_rank);
throw e;
}
// method sets the rank itself
if (rank==-1)
rank=my_rank;
// Do allocations if needed
if (start==NULL) {
// read the whole data : start from 0
allocation_done=true;
start=new int [rank];
// init to metacharacter -1
for (int i=0;i<rank;++i)
start[i]=-1;
}
if (edges==NULL) {
// read the whole data : edges from 0
allocation_done=true;
edges=new int [2];
for (int i=0;i<rank;++i)
edges[i]=-1;
}
// replace -1 metacharacter by values : 0 for starts, max_nb_val for edges
for (int i=0;i<rank;++i) {
if (start[i]==-1) start[i]=0; // use default
if (edges[i]==-1) edges[i]=data_dimension[i]-start[i]; // use default
}
/* check data selector range validity */
for (int i=0;i<rank;++i) {
if (start[i]<0 || start[i]>(data_dimension[i]-1)) {
bad_index e(start[i],data_dimension[i]-1);
throw e;
}
if (edges[i]<0 || edges[i]>(data_dimension[i]-start[i])) {
bad_index e(edges[i],(data_dimension[i]-start[i]));
throw e;
}
}
return allocation_done;
}
const EntryFormat* PARASOLFileData::get_entry_format( const string & entry_name )
{
// parameter name to index
int entry_idx=this->data->get_data_file_format()->get_entry_index(entry_name);
return get_entry_format(entry_idx);
}
const EntryFormat* PARASOLFileData::get_entry_format( const int & entry_idx )
{
// parameter name to index
if (entry_idx!=-1) {
// retrieve the entry corresponding to var_name
RecordFormat *rec_format=const_cast<RecordFormat *>(data->get_data_file_format()->get_data_format());
return &(*rec_format)[entry_idx];
} else
// no entry with this name
return NULL;
}
void * PARASOLFileData::read_data(void * data_buf, const char * var_name, int * start, int * stride, int * edges, int rank ) {
// scaling shouldn't be applied to bitfieds types and some int types
string s_var_name(var_name); // for == comparaisons
// TODO add other bitfield and a clever way to discriminate when the "read_count_data" is preferable to "read_scaled_data". This can be done using the numeric type, the field name ???
if (s_var_name==(string)("dqx") || s_var_name==(string)("utc_hour") || s_var_name==(string)("utc_minute"))
return read_count_data((unsigned char*)data_buf, var_name, start, edges, rank);
else
return read_scaled_data((double*)data_buf, var_name, start, edges, rank);
}
double * PARASOLFileData::read_scaled_data( double * scaled_data, const char * var_name, int * start, int * edges, int rank ) {
// Sets the data selector range and do the needed allocation
int var_idx=data->get_data_entry_index(var_name);
bool data_range_allocated=init_read_data_range(var_idx,start,edges,rank);
// process depends of the data type
EntryFormat *var_format=const_cast<EntryFormat *>(get_entry_format(var_name));
// allocate scaled data if needed
int buffer_size=1;
for (int i=0;i<rank;++i)
buffer_size*=edges[i];
// allocate the buffer for reading all the data
if (scaled_data==NULL) {
scaled_data=new double[buffer_size];
}
// read the scaling factors
double slope=1.;
double offset=0.;
unsigned short nb_bytes=0;
get_scaling(var_name,slope,offset,nb_bytes);
if (var_format!=NULL) {
PARASOLDataType type=var_format->type;
// the number of read values
// count buffer type depends of the entry type
if (type== UINT8) {
unsigned char* count_data=NULL;
count_data=static_cast<unsigned char*>(read_count_data(count_data,var_name,start,edges));
data->apply_scaling(scaled_data,count_data,slope,offset,buffer_size);
delete[] count_data;
} else if (type== UINT16) {
unsigned short* count_data=NULL;
count_data=static_cast<unsigned short*>(read_count_data(count_data,var_name,start,edges));
data->apply_scaling(scaled_data,count_data,slope,offset,buffer_size);
delete[] count_data;
} else if (type== INT16) {
short* count_data=NULL;
count_data=static_cast<short*>(read_count_data(count_data,var_name,start,edges));
data->apply_scaling(scaled_data,count_data,slope,offset,buffer_size);
delete[] count_data;
} else if (type== UINT32) {
unsigned int * count_data=NULL;
count_data=static_cast<unsigned int *>(read_count_data(count_data,var_name,start,edges));
data->apply_scaling(scaled_data,count_data,slope,offset,buffer_size);
delete[] count_data;
} else if (type== INT32) {
int * count_data=NULL;
count_data=static_cast<int*>(read_count_data(count_data,var_name,start,edges));
data->apply_scaling(scaled_data,count_data,slope,offset,buffer_size);
delete[] count_data;
// } else if (type== BIT2) {
// bit<2>* count_data=NULL;
// count_data=read_count_data(count_data,var_name,start,edges,dir);
// } else if (type== BIT4) {
// apply_scaling< bit<4> >(scaled_data,var_name,start,edges);
// } else if (type== BIT8) {
// apply_scaling< bit<8> >(scaled_data,var_name,start,edges);
// } else if (type== BIT16) {
// apply_scaling< bit<16> >(scaled_data,var_name,start,edges);
// } else if (type== BIT32) {
// apply_scaling< bit<32> >(scaled_data,var_name,start,edges);
// } else {
// cerr<<"In "<<__FILE__<<" at "<<__LINE__<<" : Unimplemented variable type"<<endl;
} else if (type==B2) {
short* count_data=NULL;
count_data=static_cast<short*>(read_count_data(count_data,var_name,start,edges));
data->apply_scaling(scaled_data,count_data,slope,offset,buffer_size);
delete[] count_data;
}
// for other types, it is supposed the scaling is 1
}
// free selector range allocations
if (data_range_allocated && start!=NULL && edges!=NULL) {
delete[] start;
start=NULL;
delete[] edges;
edges=NULL;
}
return scaled_data;
}
void * PARASOLFileData::read_count_data( void * buf, const char * var_name, int * start, int * edges, int rank ) {
int entry_idx=this->data->get_data_file_format()->get_entry_index(var_name);
if (entry_idx!=-1)
buf=read_count_data(buf,entry_idx,start,edges);
return buf;
}
void * PARASOLFileData::read_count_data( void * count_data, const int entry_idx, int * start, int * edges, int rank ) {
// Sets the data selector range and do the needed allocation
bool data_range_allocated=init_read_data_range(entry_idx,start,edges,rank);
// process depends of the data type
EntryFormat *var_format=const_cast<EntryFormat *>(get_entry_format(entry_idx));
if (var_format!=NULL) {
// if not specified, set the authorized data range to be read
//init_read_data_range(start,edges);
// an entry with this name exists : find its type
PARASOLDataType type=var_format->type;
if (count_data==NULL) {
// allocate the buffer for reading the data
int buffer_size=1;
for (int i=0;i<rank;++i)
buffer_size*=edges[i];
// add the length of the entry -> for char and bit types
buffer_size*=var_format->len;
// allocate the buffer
count_data=PARASOLFileReader::custom_new(count_data,buffer_size,type);
}
// read the file
data->read_data(count_data,entry_idx,start,edges,rank);
}
// free selector range allocations
if (data_range_allocated && start!=NULL && edges!=NULL) {
delete[] start, start = NULL;
delete[] edges, edges = NULL;
}
return count_data;
}
const int PARASOLFileData::get_nb_data() {
return data->get_nb_data_record();
}
const bool PARASOLFileData::is_geolocation_data_loaded( ) const {
return (line_data!=NULL && col_data!=NULL);
}
void PARASOLFileData::free_geolocation_data() {
delete[] line_data, line_data = NULL;
delete[] col_data, col_data = NULL;
}
void PARASOLFileData::load_geolocation_data() {
if (!is_geolocation_data_loaded()) {
// cout<<get_filename()<<"->load_geolocation_data"<<endl;
bool already_opened_data_file = data->is_file_loaded();
if (!already_opened_data_file)
data->open_file();
if (data->get_nb_data_record()>0) {
line_data = static_cast<unsigned short *>(read_count_data(NULL,"line"));
col_data = static_cast<unsigned short *>(read_count_data(NULL,"col"));
}
// load grid to record indices table
load_pix2data_map();
load_v_pixel();
if (!already_opened_data_file)
data->close_file();
}
}
void PARASOLFileData::load_v_pixel() {
// unused at this time
;
}
const bool PARASOLFileData::contain_data( const float & lat, const float & lon, const double & time )
{
return ( (time==-1. || contain_time(time)) && contain_location(lat,lon) );
}
const bool PARASOLFileData::contain_location( const float & lat, const float & lon )
{
// convert (lat,lon) to (lin,col)
float lat_lon[]={lat,lon};
int grid_idx[2];
geolocation_to_grid(lat_lon,grid_idx);
if (grid_idx[0]!=-1 && grid_idx[1]!=-1) {
return contain_pixel(grid_idx);
} else return false;
}
const bool PARASOLFileData::contain_pixel( const int* pix_idx )
{
int irec;
return get_index(pix_idx,irec);
}
void PARASOLFileData::get_pixel( const int & irec, int pix[] ) {
bool is_geolocation_already_loaded=is_geolocation_data_loaded();
if (!is_geolocation_already_loaded)
load_geolocation_data();
pix[0]=line_data[irec];
pix[1]=col_data[irec];
if (!is_geolocation_already_loaded)
free_geolocation_data();
}
void PARASOLFileData::igrid2idata(const int igrid[], int & irec) {
get_index (igrid, irec);
}
void PARASOLFileData::igrid2idata(const vector< int > & igrid, int & irec) {
get_index (igrid, irec);
}
const bool PARASOLFileData::get_index(const vector <int> & v_ipix, int & irec) {
std::map <PixelIndice, DataIndice>::iterator it;
if (v_ipix.size() > 2) {
// directional indice present in pixel
vector <int> _v_ipix (2);
copy (v_ipix.begin(), v_ipix.begin() + 2, _v_ipix.begin());
it = ipix2idata.find (_v_ipix);
} else
it = ipix2idata.find (v_ipix);
bool ret;
if (it == ipix2idata.end()) { // unexisting pixel
irec = -1;
ret = false;
} else {
irec = (it->second)[0];
ret = true;
}
// cout << "irec " << irec << endl;
return ret;
}
const bool PARASOLFileData::get_index( const int* pix_idx, int & irec ) {
irec = -1;
vector <int> v_ipix (2);
v_ipix [0] = pix_idx[0];
v_ipix [1] = pix_idx[1];
return get_index (v_ipix, irec);
// if (pix_idx!=NULL) {
// unsigned short lin_idx=pix_idx[0];
// unsigned short col_idx=pix_idx[1];
// bool is_geolocation_data_already_loaded=is_geolocation_data_loaded();
// if (!is_geolocation_data_already_loaded)
// load_geolocation_data();
//
// // line indexes are sorted -> do a dichotomic search to find the range for which line index =
// size_t nb_data=get_nb_data();
// pair<unsigned short*, unsigned short*> result;
// unsigned short *found_pix=NULL;
// // the ordering of the line indexes depends of the level : 1->descending, other->ascending
// if (level==1) {
// result = equal_range(line_data, line_data + nb_data, lin_idx, PARASOLFileData::greater_comparator<unsigned short>);
// } else {
// result = equal_range(line_data, line_data + nb_data, lin_idx);
// }
// int line_range[2];
// line_range[0]= distance(line_data,result.first);
// line_range[1]= distance(line_data,result.second);
// if (line_range[0]<line_range[1]) {
// // the data contain this line index -> search for the column in equal lines range
// unsigned short* start_col_search = col_data + line_range[0];
// unsigned short* end_col_search = col_data + line_range[1];
// found_pix = find(start_col_search, end_col_search, col_idx);
// if (found_pix!=NULL && found_pix!=end_col_search) {
// irec=distance(col_data,found_pix);
// found=true;
// }
// }
// if (!is_geolocation_data_already_loaded)
// free_geolocation_data();
// } else
// cerr<<"In "<<__FILE__<<" at "<<__LINE__<<" NULL pixel index not allowed"<<endl;
}
const bool PARASOLFileData::get_index( const float & lat, const float & lon, int & irec ) {
// convert (lat,lon) to (lin,col)
float lat_lon[]={lat,lon};
int pix_idx[2];
geolocation_to_grid(lat_lon,pix_idx);
if (pix_idx[0]!=-1 && pix_idx[1]!=-1) {
return get_index(pix_idx,irec);
} else
return false;
}
void PARASOLFileData::get_scaling( const string & var_name, double & slope, double & offset, unsigned short & nb_bytes ) {
int var_idx=data->get_data_entry_index(var_name);
get_scaling(var_idx, slope, offset, nb_bytes);
}
void PARASOLFileData::get_scaling( const int var_idx, double & slope, double & offset, unsigned short & nb_bytes ) {
if (!leader->is_record_loaded(SCALING_FACTORS))
leader->read_record(SCALING_FACTORS);
int NB_PARAM_WITHOUT_SCALING=6; // the number of parametres in the data record that don't have a scaling set in the leader file. It seems to be the same value for all products (either l1, l2 and l3 ones)
assert(var_idx>=0 && var_idx<(MyTools::string2int((leader->scaling_factors->nb_pix_param))+NB_PARAM_WITHOUT_SCALING));
if (var_idx>=0 && var_idx<NB_PARAM_WITHOUT_SCALING) {
EntryFormat *var_format=const_cast<EntryFormat *>(get_entry_format(var_idx)); // record entry format description
nb_bytes=PARASOLFileFormat::get_size_of(var_format->type);
slope=1.;
offset=0.;
} else if (var_idx>=NB_PARAM_WITHOUT_SCALING && var_idx<(MyTools::string2int((leader->scaling_factors->nb_pix_param))+NB_PARAM_WITHOUT_SCALING) ) {
// the 6th first parameters in the data record don't have any scaling factors set in the leader file.
int param_nb=var_idx-NB_PARAM_WITHOUT_SCALING; // what is called the "param #" in the DPC
nb_bytes=MyTools::string2int(leader->scaling_factors->nb_bytes[param_nb]);
slope=MyTools::string2float(leader->scaling_factors->slope[param_nb]);
offset=MyTools::string2float(leader->scaling_factors->offset[param_nb]);
} else { // should never arrived here
nb_bytes=0;
slope=-DBL_MAX;
offset=-DBL_MAX;
}
}
const bool PARASOLFileData::is_file_loaded( ) {
if (leader!=NULL && leader->is_file_loaded())
return true;
else if (data!=NULL && data->is_file_loaded())
return true;
else
return false;
}
const int PARASOLFileData::get_max_direction( ) {
return MAX_PARASOL_DIRECTION;
}
void PARASOLFileData::close_file( ) {
if (leader!=NULL && leader->is_file_loaded())
leader->close_file();
if (data!=NULL && data->is_file_loaded())
data->close_file();
}
vector<int> PARASOLFileData::get_data_dimension( const string & var_name ) {
int entry_idx = data->get_data_entry_index(var_name); // entry index in the record
return get_data_dimension(entry_idx);
}
vector<int> PARASOLFileData::get_data_dimension( const int entry_idx ) {
vector<int> v;
v.reserve(2); // at most 2 dimensions
vector <EntryBlock> v_entry_block = data->get_entry_block(DATA, entry_idx);
for (vector <EntryBlock>::iterator it = v_entry_block.begin() ; it != v_entry_block.end() ; ++it) {
v.push_back(it->ntime);
}
return v;
}
void PARASOLFileData::idata2igrid(const int irec, int igrid[2]) {
// validate record_index
if (irec < 0 || irec >= data->get_nb_data_record()) {
g_exception e(__FILE__,__LINE__, "Invalid record indice " + MyTools::to_string(irec));
throw e;
}
if (! is_geolocation_data_loaded()) {
g_exception e(__FILE__,__LINE__, "Geolocation data not loaded. Consider to call load_geolocation_data()");
throw e;
}
igrid [0] = line_data [irec];
igrid [1] = col_data [irec];
}
void PARASOLFileData::get_record_coord( const int irec, float &lat, float &lon ) {
// read the col,line indexes of the record
int igrid[2];
idata2igrid(irec, igrid);
// col,line indexes to lat,lon
float lat_range[] = {-1.,-1.}, lon_range[] = {-1.,-1.};
grid_to_geolocation(igrid, lat_range, lon_range);
// compute the coordinates of the centre of the pixel
lat = 0.5 * (lat_range[0] + lat_range[1]);
lon = 0.5 * (lon_range[0] + lon_range[1]);
}
const bool PARASOLFileData::contain_location(const float &lat,const float &lon, const double &tolerance){
// convert (lat,lon) to (lin,col)
float lat_lon[]={lat,lon};
float lat_minmax[]={lat,lon};
float lon_minmax[]={lat,lon};
int grid_idx[2];
geolocation_to_grid(lat_lon,grid_idx);
if (grid_idx[0]!=-1 && grid_idx[1]!=-1) {
grid_to_geolocation(grid_idx,lat_minmax,lon_minmax);
float lat1 = (lat_minmax[0]+lat_minmax[1])/2.0;
float lon1 = (lon_minmax[0]+lon_minmax[1])/2.0;
if(lat1<lat-tolerance || lat+tolerance<lat1 ) return false;
if(lon1<lon-tolerance || lon+tolerance<lon1 ) return false;
return contain_pixel(grid_idx);
} else return false;
}
const bool PARASOLFileData::contain_data( const float &lat, const float &lon, const double &time, const double &tolerance )
{
return ( (time==-1. || contain_time(time)) && contain_location(lat,lon,tolerance) );
}
const bool PARASOLFileData::get_index( const float & lat, const float & lon, int & irec, const float tolerance)
{
// convert (lat,lon) to (lin,col)
float lat_lon[]={lat,lon};
float lat_minmax[]={lat,lon};
float lon_minmax[]={lat,lon};
int pix_idx[2];
irec=-1;
geolocation_to_grid(lat_lon,pix_idx);
if (pix_idx[0]!=-1 && pix_idx[1]!=-1) {
grid_to_geolocation(pix_idx,lat_minmax,lon_minmax);
float lat1 = (lat_minmax[0]+lat_minmax[1])/2.0;
float lon1 = (lon_minmax[0]+lon_minmax[1])/2.0;
if(lat1<lat-tolerance || lat+tolerance<lat1 ) return false;
if(lon1<lon-tolerance || lon+tolerance<lon1 ) return false;
return get_index(pix_idx,irec);
} else
return false;
}
const float PARASOLFileData::get_nearest_point_distance( const float & lat, const float & lon, const float coloc_tolerance ) {
// convert (lat,lon) to (lin,col)
float lat_lon[]={lat,lon};
float lat_minmax[]={lat,lon};
float lon_minmax[]={lat,lon};
int pix_idx[2];
double nearest_point_distance=-1.;
geolocation_to_grid(lat_lon,pix_idx);
if (pix_idx[0]!=-1 && pix_idx[1]!=-1) {
grid_to_geolocation(pix_idx,lat_minmax,lon_minmax);
float lat1 = (lat_minmax[0]+lat_minmax[1])/2.0;
float lon1 = (lon_minmax[0]+lon_minmax[1])/2.0;
if(lat1<lat-coloc_tolerance || lat+coloc_tolerance<lat1 ) return -1.;
if(lon1<lon-coloc_tolerance || lon+coloc_tolerance<lon1 ) return -1.;
nearest_point_distance=sqrt(pow(lat1-lat,2) + pow(lon1-lon,2));
return nearest_point_distance;
} else
return -1.;
}
void PARASOLFileData::load_pix2data_map() {
free_pix2data_map();
// // number of available directions
// char * v_ni = (char *)(read_count_data(NULL, "ni"));
// int nrec = get_nb_data();
// int idir, ndir;
// vector <int> ipix(3); // pixel indice in PARASOL grid [irow, icolumn, idirection]
// vector <int> idata(2); // data indice [irecord, idirection]
// // buffer linear indice
// int ibuf;
// for (int irec = 0 ; irec < nrec ; ++irec) {
// // set pixel indices [irow, icolumn, idirection]
// ndir = (int)(v_ni[irec]);
// ipix[0] = line_data [irec];
// ipix[1] = col_data [irec];
// idata [0] = irec;
// for (idir = 0 ; idir < ndir ; ++idir) {
// ibuf = irec * nb_parasol_direction + idir;
// ipix[2] = idir;
// idata [1] = idir;
// ipix2idata [ipix] = idata;
// }
// }
int nrec = get_nb_data();
vector <int> ipix(2); // pixel indice in PARASOL grid [irow, icolumn]
vector <int> idata(1); // data indice [irecord]
for (int irec = 0 ; irec < nrec ; ++irec) {
// set pixel indices [irow, icolumn, idirection]
ipix[0] = line_data [irec];
ipix[1] = col_data [irec];
idata [0] = irec;
ipix2idata [ipix] = idata;
}
}
bool PARASOLFileData::is_earth_grid(const int irow, const int icol) const {
// check row range
if (irow < 1 || irow > sz_grid[0])
return false;
// check col range. Validity dependqs of the row indice
int ncol = sz_grid[1];
int ni = get_ncol (irow)/2;
int icol_min = ncol/2 - ni + 1;
int icol_max = ncol/2 + ni;
if (icol < icol_min || icol > icol_max)
return false;
return true;
}
void PARASOLFileData::get_v_icol_neighbours(const int irow, const int icol, int dicol_min, int dicol_max, vector< int > & v_icol){
check_grid_coord (irow, icol);
int ncol = sz_grid[1];
int ni = get_ncol (irow)/2;
int icol_min = ncol/2 - ni + 1;
int icol_max = ncol/2 + ni;
// cout << "ni " << ni << " icol_min " << icol_min << " icol_max " << icol_max << endl;
v_icol.clear();
if ((dicol_max - dicol_min + 1) >= (2 * ni) ) {
// all columns requested
MyTools::arange (v_icol, icol_min, icol_max + 1);
} else {
vector <int> v_tmp(0);
// check icol min range
int imin = (icol + dicol_min);
if (imin >= icol_min) { // change date meridian not crossed
MyTools::arange (v_icol, imin, icol + 1);
} else { // change date meridian crossed
MyTools::arange (v_icol, icol_min, icol + 1);
// report remaining columns on east side
int d = (icol_min - imin);
MyTools::arange (v_tmp, icol_max - d + 1, icol_max + 1);
MyTools::push_back (v_tmp, v_icol);
}
// check icol max range
int imax = (icol + dicol_max);
if (imax <= icol_max) { // change date meridian not crossed
MyTools::arange (v_tmp, icol + 1, imax + 1);
MyTools::push_back (v_tmp, v_icol);
} else { // change date meridian crossed
MyTools::arange (v_tmp, icol + 1, icol_max + 1);
MyTools::push_back (v_tmp, v_icol);
// report remaining columns on west side
int d = (imax - icol_max);
MyTools::arange (v_tmp, icol_min, icol_min + d);
MyTools::push_back (v_tmp, v_icol);
}
}
sort(v_icol.begin(), v_icol.end());
// cout << "irow " << irow << " icol " << icol << " dicol [" << dicol_min << ", " << dicol_max << "] " << MyTools::vec2str(v_icol) << endl;
}
void PARASOLFileData::get_neighbours_coord(const int irow, const int icol, int dirow_min, int dirow_max, int dicol_min, int dicol_max, vector< int > & v_irow, vector< int > & v_icol) {
check_grid_coord (irow, icol);
// int _d_icol = min (d_icol, icol_max / 2);
// cout << "-----------------------------------------" << endl;
// cout << "ipix [" << irow << ", " << icol << "] d_irow " << d_irow << " d_icol " << d_icol << endl;
// cout << "irow range [" << irow_min << "," << irow_max << "] " << "icol range [" << icol_min << "," << icol_max << "]" << endl;
// rows
v_irow.clear();
int i;
int irow_min = max (1, irow + dirow_min);
int irow_max = min (sz_grid[0], irow + dirow_max);
for (i = irow_min ; i <= irow_max ; ++i)
v_irow.push_back(i);
// search the row that has the biggest number of columns (ie the nearest one from equator, at indice sz_grid[0]/2 + 0.5)
float irow_equator = sz_grid[0]/2 + 0.5;
float dirow_near = 999999.;
float dirow_equator;
int irow_near = -9999;
for (i = irow_min ; i <= irow_max ; ++i) {
dirow_equator = fabs (i - irow_equator);
if (dirow_equator < dirow_near) {
dirow_near = dirow_equator;
irow_near = i;
}
}
if (irow_near == -9999) {
g_exception e(__FILE__,__LINE__, "Invalid row indice range [" + MyTools::to_string(irow_min) + ", " + MyTools::to_string(irow_max) + "]");
throw e;
}
// columns
v_icol.clear();
get_v_icol_neighbours (irow_near, icol, dicol_min, dicol_max, v_icol);
// cout << "v_irow " << MyTools::vec2str(v_irow) << endl;
// cout << "v_icol " << MyTools::vec2str(v_icol) << endl;
}
void PARASOLFileData::check_grid_coord(const int irow, const int icol) const {
// int icol_max = get_ncol(irow);
// cout << irow << " " << icol << " " << icol_max << endl;
g_exception e;
if (irow < 1 || irow > sz_grid[0]) {
e.set (__FILE__, __LINE__, "Invalid row indice " + MyTools::to_string (irow) + ". Must be in [1, " + MyTools::to_string (sz_grid[0]) + "]");
throw e;
}
if (icol < 1 || icol > sz_grid[1]) {
e.set (__FILE__, __LINE__, "Invalid column indice " + MyTools::to_string (icol) + ". Must be in [1, " + MyTools::to_string (sz_grid[1]) + "]");
throw e;
}
}
void PARASOLFileData::free_viewing_directions_data() {
delete [] v_alt_view, v_alt_view = NULL;
delete [] v_saa, v_saa = NULL;
delete [] v_vza, v_vza = NULL;
delete [] v_raa, v_raa = NULL;
delete [] v_ni, v_ni = NULL;
delete [] v_sn, v_sn = NULL;
delete [] v_x_sat, v_x_sat = NULL;
delete [] v_y_sat, v_y_sat = NULL;
delete [] v_z_sat, v_z_sat = NULL;
}
void PARASOLFileData::load_viewing_directions_data() {
g_exception e;
if (is_viewing_directions_data_loaded()) {
e.set (__FILE__, __LINE__, "Viewing directions data have already been loaded");
throw e;
}
// requested parametres depend of the level
int sz;
switch (level) {
case 1:
// position of satellite set in the technological parameters
sz = NB_MAX_PARASOL_SEQUENCE * NB_MAX_PARASOL_IMAGE_PER_SEQUENCE;
v_x_sat = new double [sz];
v_y_sat = new double [sz];
v_z_sat = new double [sz];
leader->read_record(TECHNOLOGICAL_PARAMETERS);
leader->techno_param->get_xyz_sat(v_x_sat, v_y_sat, v_z_sat);
leader->free_record(TECHNOLOGICAL_PARAMETERS);
v_alt_view = (short *)(read_count_data(NULL, "alt"));
v_ni = (unsigned char *)(read_count_data(NULL, "ni"));
v_sn = (unsigned char *)(read_count_data(NULL, "seq_nb"));
break;
case 2:
v_alt_view = (short *)(read_count_data(NULL, "alt"));
v_saa = (double *)(read_data(NULL, "phis"));
v_vza = (double *)(read_data(NULL, "thetav"));
v_raa = (double *)(read_data(NULL, "phi"));
v_ni = (unsigned char *)(read_count_data(NULL, "ni"));
break;
case 3:
e.set( __FILE__ , __LINE__ , "Method not implemented for L3 products" );
throw e;
break;
default:
e.set( __FILE__ , __LINE__ , "Invalid level " + MyTools::to_string(level) );
throw e;
break;
};
}
bool PARASOLFileData::is_viewing_directions_data_loaded() {
return (v_alt_view != NULL);
}
void PARASOLFileData::get_viewing_directions(const vector< int > & ipix, vector< Observation > & v_obs) {
g_exception e;
if (! is_viewing_directions_data_loaded()) {
e.set (__FILE__, __LINE__, "Viewing directions data must have been previously loaded. Consider calling load_viewing_directions_data() before");
throw e;
}
v_obs.clear();
v_obs.reserve (nb_parasol_direction);
// retrieve record indice
int irec;
bool is_valid_pix = get_index(ipix, irec);
if (is_valid_pix) {
// cout << "ipix " << MyTools::vec2str(ipix) << " -> irec " << irec << endl;
// available parametres depend of the level
switch (level) {
case 1:
_get_viewing_directions_l1(ipix, irec, v_obs);
break;
case 2:
_get_viewing_directions_l2(ipix, irec, v_obs);
break;
case 3:
e.set( __FILE__ , __LINE__ , "Method not implemented for L3 products" );
throw e;
break;
default:
e.set( __FILE__ , __LINE__ , "Invalid level " + MyTools::to_string(level) );
throw e;
break;
};
}
}
void PARASOLFileData::_get_viewing_directions_l1(const vector< int > & ipix, const int irec, vector< Observation > & v_obs) {
// --- for each pixel identified by its indices [irow, icolumn, idirection], read
int idir = 0, ndir = nb_parasol_direction;
// latitude, longitude and altitude at center of viewed pixel
float lat_view, lon_view, alt_view;
// carthesian coordinates in ECR system of satellite position.
double x_sat, y_sat, z_sat; // as double
// sequence number
unsigned short sn;
// points representing the location of the viewed pixel and the satellite
Geodetic::Point3D point;
// Geocentric::Point3D point;
Carthesian::Point3D point_sat, point_view;
vector <int> igrid(3); // grid indice : [irow, icol, idir]
vector <int> idata(1); // data buffer indice [irec]
idata [0] = irec;
// buffer linear indice
int ibuf;
int i_img = 4; // indice of the middle of sequence image
// set pixel indices [irow, icolumn, idirection]
igrid [0] = line_data [irec];
igrid [1] = col_data [irec];
// set the latitude, longitude and altitude of the viewed pixel (center)
get_record_coord(irec, lat_view, lon_view);
alt_view = v_alt_view [irec];
point.set (lat_view, lon_view, alt_view);
point_view = point.to_carthesian();
// point_view = point.to_carthesian();
// cout << point << endl;
// cout << point.to_carthesian(Earth::WGS84_MODEL) << endl;
// cout << point.to_carthesian() << endl;
// cout << point.to_carthesian(Earth::WGS84_MODEL).to_geocentric() << endl;
// cout << point.to_carthesian().to_geocentric() << endl;
// cout << point_view << endl;
// cout << point_view.to_geocentric() << endl;
// read satellite position in each direction
ndir = (int)(v_ni[irec]);
Observation obs;
Carthesian::Segment3D viewing;
// set directions to read range
int idir_start = 0, idir_end = ndir;
if (ipix.size() > 2) {
idir_start = ipix[2], idir_end = idir_start + 1;
if (idir_start >= ndir) // direction does not exists : skip next loop
idir_start = idir_end;
}
// for (sn = 0 ; sn < 130 ; ++sn) {
// for (i_img = 0 ; i_img < 9 ; ++i_img) {
// x_sat = MyTools::to_num<double>((leader->techno_param->x_sat) [sn][i_img]) * 1e3; // km -> m
// y_sat = MyTools::to_num<double>((leader->techno_param->y_sat) [sn][i_img]) * 1e3; // km -> m
// z_sat = MyTools::to_num<double>((leader->techno_param->z_sat) [sn][i_img]) * 1e3; // km -> m
// printf ("* seq = %d im = %d sat_pos = [%f, %f, %f]\n", sn + 1, i_img + 1, x_sat, y_sat, z_sat);
// }
// }
for (idir = idir_start ; idir < idir_end ; ++idir) {
ibuf = irec * nb_parasol_direction + idir;
// --- read x, y, z coords of satellite. Central image is considered for all filters (image number = 4)
sn = (unsigned short)(v_sn [ibuf]) - 1; // - 1 because starts at 1 in file, and here is used as indice from 0
// x_sat = MyTools::to_num<double>((leader->techno_param->x_sat) [sn][i_img]) * 1e3; // km -> m
// y_sat = MyTools::to_num<double>((leader->techno_param->y_sat) [sn][i_img]) * 1e3; // km -> m
// z_sat = MyTools::to_num<double>((leader->techno_param->z_sat) [sn][i_img]) * 1e3; // km -> m
ibuf = sn * NB_MAX_PARASOL_IMAGE_PER_SEQUENCE + i_img;
x_sat = v_x_sat [ibuf];
y_sat = v_y_sat [ibuf];
z_sat = v_z_sat [ibuf];
point_sat.set (x_sat, y_sat, z_sat);
// --- construct viewing segment and store it
viewing.set (point_sat, point_view);
igrid[2] = idir;
obs.set (lat_view, lon_view, alt_view, -1, idata, igrid, viewing);
v_obs.push_back (obs);
// printf ("--- irec = %d ; pixel [%d, %d, %d] ---\n", irec, igrid [0], igrid [1], idir);
// printf (" * ni = %d sn = %d alt_view = %f\n", ndir, sn + 1, alt_view);
// point = point_view.to_geodetic();
// printf (" * view : (%f, %f, %f) ; (%f, %f, %f)\n", point.lat, point.lon, point.alt, point_view.x, point_view.y, point_view.z);
// point = point_sat.to_geodetic();
// printf (" * sat : (%f, %f, %f) ; (%f, %f, %f)\n", point.lat, point.lon, point.alt, point_sat.x, point_sat.y, point_sat.z);
// printf (" * segment : distance %f\n", viewing.length());
}
}
// void PARASOLFileData::load_viewing_directions() {
// g_exception e;
// if (! is_viewing_directions_data_loaded()) {
// e.set (__FILE__, __LINE__, "Viewing directions data must have been previously loaded. Considr calling load_viewing_directions_data() before");
// throw e;
// }
// // assert(is_geolocation_data_loaded());
// // init container if needed
// viewing_directions.clear();
//
// // available parametres depend of the level
// switch (level) {
// case 1:
// _load_viewing_directions_l1();
// break;
// case 2:
// _load_viewing_directions_l2();
// break;
// case 3:
// e.set( __FILE__ , __LINE__ , "Method not implemented for L3 products" );
// throw e;
// break;
// default:
// e.set( __FILE__ , __LINE__ , "Invalid level " + MyTools::to_string(level) );
// throw e;
// break;
// };
// }
//
// void PARASOLFileData::_load_viewing_directions_l1() {
// // --- for each pixel identified by its indices [irow, icolumn, idirection], read
// int irec = 0, nrec = get_nb_data();
// int irow, icol;
// int idir = 0, ndir = nb_parasol_direction;
// // latitude, longitude and altitude at center of viewed pixel
// float lat_view, lon_view, alt_view;
// // carthesian coordinates in ECR system of satellite position.
// double x_sat, y_sat, z_sat; // as double
// // sequence number
// short sn;
//
// // points representing the location of the viewed pixel and the satellite
// Geodetic::Point3D point;
// Carthesian::Point3D point_sat, point_view;
// vector <int> ipix(3);
// // buffer linear indice
// int ibuf;
// int i_img = 4; // indice of the middle of sequence image
// for (irec = 0 ; irec < nrec ; ++irec) {
// // set pixel indices [irow, icolumn, idirection]
// irow = line_data [irec];
// icol = col_data [irec];
// // set the latitude, longitude and altitude of the viewed pixel (center)
// get_record_coord(irec, lat_view, lon_view);
// alt_view = v_alt_view [irec];
//
// point.set (lat_view, lon_view, alt_view);
// point_view = point.to_carthesian();
//
// // read satellite position in each direction
// ndir = (int)(v_ni[irec]);
// for (idir = 0 ; idir < ndir ; ++idir) {
// ibuf = irec * nb_parasol_direction + idir;
// // --- read x, y, z coords of satellite. Central image is considered for all filters (image number = 4)
// sn = (short)(v_sn [ibuf]) - 1; // - 1 because starts at 1 in file, and here is used as indice from 0
// x_sat = MyTools::to_num<double>((leader->techno_param->x_sat) [sn][i_img]) * 1e3; // km -> m
// y_sat = MyTools::to_num<double>((leader->techno_param->y_sat) [sn][i_img]) * 1e3; // km -> m
// z_sat = MyTools::to_num<double>((leader->techno_param->z_sat) [sn][i_img]) * 1e3; // km -> m
// point_sat.set (x_sat, y_sat, z_sat);
//
// // --- construct viewing segment and store it
// ipix[0] = irow, ipix[1] = icol, ipix[2] = idir;
// viewing_directions[ipix] = Carthesian::Segment3D (point_sat, point_view);
//
// // if (irec == 97552) {
// // printf ("--- irec = %d ; pixel [%d, %d, %d] ---\n", irec, irow, icol, idir);
// // printf (" * ni = %d sn = %d\n", ndir, sn + 1);
// // point = point_view.to_geocentric();
// // printf (" * view : (%f, %f, %f) ; (%f, %f, %f)\n", point.lat, point.lon, point.alt, point_view.x, point_view.y, point_view.z);
// // point = point_sat.to_geocentric();
// // printf (" * sat : (%f, %f, %f) ; (%f, %f, %f)\n", point.lat, point.lon, point.alt, point_sat.x, point_sat.y, point_sat.z);
// // printf (" * segment : distance %f\n", viewing_directions[ipix].length());
// // }
//
// }
// }
// // delete[] v_alt_view;
// // delete[] v_ni;
// // delete[] v_sn;
// }
void PARASOLFileData::_get_viewing_directions_l2(const vector< int > & ipix, const int irec, vector< Observation > & v_obs) {
// --- for each pixel identified by its indices [irow, icolumn, idirection], read
int idir = 0, ndir = nb_parasol_direction;
// latitude, longitude and altitude at center of viewed pixel
float lat_view, lon_view, alt_view;
// fixed satellite altitude
float alt_sat = satellite_altitude;
// relative azimuth, viewing azimuth, viewing zenith and solar azimuth angles
double raa, vaa, vza, saa;
// points representing the location of the viewed pixel and the satellite
Geodetic::Point3D point_view;
Carthesian::Point3D point_sat, point_view_xyz;
vector <int> igrid(3); // grid indice : [irow, icol, idir]
vector <int> idata(1); // data buffer indice [irec]
idata [0] = irec;
// buffer linear indice
int ibuf;
// set pixel indices [irow, icolumn, idirection]
igrid [0] = line_data [irec];
igrid [1] = col_data [irec];
// set the latitude, longitude and altitude of the viewed pixel (center)
get_record_coord(irec, lat_view, lon_view);
alt_view = v_alt_view [irec];
point_view.set (lat_view, lon_view, alt_view);
point_view_xyz = point_view.to_carthesian();
// constructs viewin points
saa = v_saa [irec];
// compute directionnal paramtres
ndir = (int)(v_ni[irec]);
Observation obs;
Carthesian::Segment3D viewing;
for (idir = 0 ; idir < ndir ; ++idir) {
ibuf = irec * nb_parasol_direction + idir;
// --- compute lat, lon, alt of the satellite from viewing directions
raa = v_raa [ibuf];
vaa = (saa - raa);
vza = v_vza [ibuf];
// --- constructs satellite position from viewing angles
Viewing::earth_to_sat_from_alt(point_view, vza, vaa, alt_sat, point_sat);
// --- construct viewing segment and store it
viewing.set (point_sat, point_view_xyz);
igrid[2] = idir;
obs.set (lat_view, lon_view, alt_view, -1, idata, igrid, viewing);
v_obs.push_back (obs);
// Geocentric::Point3D point_geo = point_sat.to_geocentric();
// printf ("--- irec = %d ; pixel [%d, %d, %d] ---\n", irec, irow, icol, idir);
// printf (" * ni = %d ; vra = %f ; saa = %f ; vaa = %f ; vza = %f\n", ndir, v_raa [ibuf], v_saa [irec], v_saa [irec] - v_raa [ibuf], v_vza [ibuf]);
// // printf (" dlat = %f ; dlon = %f ; alpha = %f\n", dlat, dlon, degrees (alpha));
// printf (" * view : (%f, %f, %f) ; (%f, %f, %f)\n", lat_view, lon_view, alt_view, point_view_xyz.x, point_view_xyz.y, point_view_xyz.z);
// printf (" * sat : (%f, %f, %f) ; (%f, %f, %f)\n", point_geo.lat, point_geo.lon, point_geo.alt, point_sat.x, point_sat.y, point_sat.z);
// printf (" * segment : distance %f\n", viewing_directions[ipix].length());
}
}
// void PARASOLFileData::_load_viewing_directions_l2() {
// // --- preload needed variables
// // // viewed pixel altitude
// // short * v_alt_view = (short *)(read_count_data(NULL, "alt"));
// // // solar azimuth angles. 1 per record
// // double * v_saa = (double *)(read_data(NULL, "phis"));
// // // viewing zenith angles. N_directions per record
// // double * v_vza = (double *)(read_data(NULL, "thetav"));
// // // relative azimuth angles. N_directions per record
// // double * v_raa = (double *)(read_data(NULL, "phi"));
// // // number of available directions
// // char * v_ni = (char *)(read_count_data(NULL, "ni"));
//
// // --- for each pixel identified by its indices [irow, icolumn, idirection], read
// int irec = 0, nrec = get_nb_data();
// int irow, icol;
// int idir = 0, ndir = nb_parasol_direction;
// // latitude, longitude and altitude at center of viewed pixel
// float lat_view, lon_view, alt_view;
// // fixed satellite altitude
// float alt_sat = satellite_altitude;
// // relative azimuth, viewing azimuth, viewing zenith and solar azimuth angles
// double raa, vaa, vza, saa;
// // points representing the location of the viewed pixel and the satellite
// Geocentric::Point3D point_view;
// Carthesian::Point3D point_sat;
// vector <int> ipix(3);
// // buffer linear indice
// int ibuf;
// for (irec = 0 ; irec < nrec ; ++irec) {
// // set pixel indices [irow, icolumn, idirection]
// irow = line_data [irec];
// icol = col_data [irec];
// // set the latitude, longitude and altitude of the viewed pixel (center)
// get_record_coord(irec, lat_view, lon_view);
// alt_view = v_alt_view [irec];
// point_view.set (lat_view, lon_view, alt_view);
// // constructs viewin points
// saa = v_saa [irec];
// // compute directionnal paramtres
// ndir = (int)(v_ni[irec]);
// for (idir = 0 ; idir < ndir ; ++idir) {
// ibuf = irec * nb_parasol_direction + idir;
// // --- compute lat, lon, alt of the satellite from viewing directions
// raa = v_raa [ibuf];
// vaa = (saa - raa);
// vza = v_vza [ibuf];
// // --- constructs satellite position from viewing angles
// Viewing::earth_to_sat_from_alt(point_view, vza, vaa, alt_sat, point_sat);
// // --- construct viewing segment and store it
// ipix[0] = irow, ipix[1] = icol, ipix[2] = idir;
// viewing_directions[ipix] = Carthesian::Segment3D (point_sat, point_view.to_carthesian());
//
// if (irec == 499) {
// Geocentric::Point3D point_geo = point_sat.to_geocentric();
// Carthesian::Point3D point_car = point_view.to_carthesian();
// printf ("--- irec = %d ; pixel [%d, %d, %d] ---\n", irec, irow, icol, idir);
// printf (" * ni = %d ; vra = %f ; saa = %f ; vaa = %f ; vza = %f\n", ndir, v_raa [ibuf], v_saa [irec], v_saa [irec] - v_raa [ibuf], v_vza [ibuf]);
// // printf (" dlat = %f ; dlon = %f ; alpha = %f\n", dlat, dlon, degrees (alpha));
// printf (" * view : (%f, %f, %f) ; (%f, %f, %f)\n", lat_view, lon_view, alt_view, point_car.x, point_car.y, point_car.z);
// printf (" * sat : (%f, %f, %f) ; (%f, %f, %f)\n", point_geo.lat, point_geo.lon, point_geo.alt, point_sat.x, point_sat.y, point_sat.z);
// printf (" * segment : distance %f\n", viewing_directions[ipix].length());
// }
//
// }
// }
// // delete[] v_alt_view;
// // delete[] v_saa;
// // delete[] v_vza;
// // delete[] v_raa;
// }
unsigned char* PARASOLFileData::get_ndir_data() const {
return v_ni;
}