/*************************************************************************** * Copyright (C) 2011 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 "dardarfiledata.h" const float DARDARFileData::colocation_tolerance=0.01; const int16 DARDARFileData::dardar_height [] = { 25080, 25020, 24960, 24900, 24840, 24780, 24720, 24660, 24600, 24540, 24480, 24420, 24360, 24300, 24240, 24180, 24120, 24060, 24000, 23940, 23880, 23820, 23760, 23700, 23640, 23580, 23520, 23460, 23400, 23340, 23280, 23220, 23160, 23100, 23040, 22980, 22920, 22860, 22800, 22740, 22680, 22620, 22560, 22500, 22440, 22380, 22320, 22260, 22200, 22140, 22080, 22020, 21960, 21900, 21840, 21780, 21720, 21660, 21600, 21540, 21480, 21420, 21360, 21300, 21240, 21180, 21120, 21060, 21000, 20940, 20880, 20820, 20760, 20700, 20640, 20580, 20520, 20460, 20400, 20340, 20280, 20220, 20160, 20100, 20040, 19980, 19920, 19860, 19800, 19740, 19680, 19620, 19560, 19500, 19440, 19380, 19320, 19260, 19200, 19140, 19080, 19020, 18960, 18900, 18840, 18780, 18720, 18660, 18600, 18540, 18480, 18420, 18360, 18300, 18240, 18180, 18120, 18060, 18000, 17940, 17880, 17820, 17760, 17700, 17640, 17580, 17520, 17460, 17400, 17340, 17280, 17220, 17160, 17100, 17040, 16980, 16920, 16860, 16800, 16740, 16680, 16620, 16560, 16500, 16440, 16379, 16320, 16260, 16200, 16139, 16080, 16020, 15960, 15900, 15840, 15780, 15720, 15660, 15600, 15540, 15480, 15420, 15360, 15300, 15240, 15180, 15120, 15060, 15000, 14940, 14880, 14820, 14760, 14700, 14640, 14580, 14520, 14460, 14400, 14340, 14280, 14220, 14160, 14100, 14040, 13980, 13920, 13860, 13800, 13740, 13680, 13620, 13560, 13500, 13440, 13380, 13320, 13260, 13200, 13140, 13080, 13020, 12960, 12900, 12840, 12780, 12720, 12660, 12600, 12540, 12480, 12420, 12360, 12300, 12240, 12180, 12120, 12060, 12000, 11940, 11880, 11820, 11760, 11700, 11640, 11580, 11520, 11460, 11400, 11340, 11280, 11220, 11160, 11100, 11040, 10980, 10920, 10860, 10800, 10740, 10680, 10620, 10560, 10500, 10440, 10380, 10320, 10260, 10200, 10140, 10080, 10020, 9960, 9900, 9840, 9780, 9720, 9660, 9600, 9540, 9480, 9420, 9360, 9300, 9240, 9180, 9120, 9060, 9000, 8940, 8880, 8820, 8760, 8700, 8640, 8580, 8520, 8460, 8400, 8340, 8280, 8220, 8160, 8100, 8040, 7980, 7920, 7860, 7800, 7740, 7680, 7620, 7560, 7500, 7440, 7380, 7320, 7260, 7200, 7140, 7080, 7020, 6960, 6900, 6840, 6780, 6720, 6660, 6600, 6540, 6480, 6420, 6360, 6300, 6240, 6180, 6120, 6060, 6000, 5940, 5880, 5820, 5760, 5700, 5640, 5580, 5520, 5460, 5400, 5340, 5280, 5220, 5160, 5100, 5040, 4980, 4920, 4860, 4800, 4740, 4680, 4620, 4560, 4500, 4440, 4380, 4320, 4260, 4200, 4140, 4080, 4020, 3960, 3900, 3840, 3780, 3720, 3660, 3600, 3540, 3480, 3420, 3360, 3300, 3240, 3180, 3120, 3060, 3000, 2940, 2880, 2820, 2760, 2700, 2640, 2580, 2520, 2460, 2400, 2340, 2280, 2220, 2160, 2100, 2040, 1980, 1920, 1860, 1800, 1740, 1680, 1620, 1560, 1500, 1440, 1380, 1320, 1260, 1200, 1140, 1080, 1020, 960, 900, 840, 780, 720, 660, 600, 540, 480, 420, 360, 300, 240, 180, 120, 60, -0, -60, -120, -180, -240, -300, -360, -420, -480, -540, -600, -660, -720, -780, -840, -900, -960, -1020 }; DARDARFileData::DARDARFileData(const string & name, const string & mode) : FileData(name,mode),SatelliteFileData(name,mode),HDFFileData(name,mode) { init(); } DARDARFileData::~DARDARFileData() { free_geolocation_data(); free_hdf_metadata(); free_hdf_file(); } void DARDARFileData::init() { lat_data=NULL; lon_data=NULL; time_data=NULL; product=""; product_id = UNDEFINED; time_coverage=-1.; orbit_nb=-1; // parse the file name to know the product, the acquisition start time... string short_filename = get_tail(get_name()); parse_filename(short_filename); // Read characteristics in the hdf file load_hdf_file(); // set the attibutes set_lat_lon_index_max(); set_time_coverage(); free_hdf_file(); } void DARDARFileData::parse_filename(const string & short_filename) { check_filename(short_filename); using namespace MyTools; // for to_num function // Product type size_t pos = short_filename.find_first_of("_"); product = short_filename.substr (0, pos); if ( product == "DARDAR-MASK" ) product_id = DARDAR_MASK; else if ( product == "DARDAR-CLOUD" ) product_id = DARDAR_CLOUD; else { string msg ( "Unsupported product " + product ) ; g_exception e( __FILE__ , __LINE__ , msg ); throw e; } // version ++pos; version = short_filename.substr( pos, 6); // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time) pos += 7; // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time) string s_date = short_filename.substr (pos, 13); date->set_date_str (s_date, "%Y%j%H%M%S"); // orbit number ++pos; orbit_nb = to_num ( short_filename.substr (pos, 5) ); // geolocation SDS names depend of the product switch ( product_id ) { case DARDAR_MASK : latitude_name = "CLOUDSAT_Latitude"; longitude_name = "CLOUDSAT_Longitude"; time_name = "CLOUDSAT_TAI_Time"; break; case DARDAR_CLOUD : latitude_name = "latitude"; longitude_name = "longitude"; time_name = "time"; break; default : string msg ( "Unsupported product " + product ) ; g_exception e( __FILE__ , __LINE__ , msg ); throw e; break; } } bool DARDARFileData::check_filename(const string & short_filename) { int infilename_size=short_filename.size(); using namespace MyTools; // for has_type function // Product type size_t pos = short_filename.find_first_of("_"); string s_prod_id = short_filename.substr (0, pos); if ( s_prod_id != "DARDAR-MASK" && s_prod_id != "DARDAR-CLOUD" ) { invalid_filename e(__FILE__,__LINE__,short_filename,"DARDAR"); throw e; } // version ++pos; string s_version = short_filename.substr( pos, 6); // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time) pos += 7; string s_date = short_filename.substr (pos, 13); if (!has_type(s_date)) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } // orbit number ++pos; string s_orbit_nb = short_filename.substr (pos, 5); if (!has_type(s_orbit_nb)) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } // extension string s_extension=short_filename.substr(infilename_size-4,4); if (s_extension!=string(".hdf")) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } return true; } void DARDARFileData::set_time_coverage( ) { float64 * d_time_data = NULL; float32 * f_time_data = NULL; switch ( product_id ) { case DARDAR_MASK : d_time_data = (float64*) read_data ( d_time_data, time_name.c_str() ); time_coverage = d_time_data [lat_lon_index_max] - d_time_data [0]; break; case DARDAR_CLOUD : f_time_data = (float32*) read_data ( f_time_data, time_name.c_str() ); time_coverage = f_time_data [lat_lon_index_max] - f_time_data [0]; break; default: string msg ( "Unsupported product " + product ) ; g_exception e( __FILE__ , __LINE__ , msg ); throw e; break; } delete[] d_time_data; delete[] f_time_data; } void DARDARFileData::free_geolocation_data() { delete[] lat_data; lat_data=NULL; delete[] lon_data; lon_data=NULL; delete[] time_data; time_data=NULL; } void DARDARFileData::load_geolocation_data() { if (!is_geolocation_data_loaded()) { // hdf_file is needed, but the method assume to conserve the caller's state bool hdf_file_already_loaded = is_hdf_file_loaded(); if (!hdf_file_already_loaded) load_hdf_file(); //-- Read lat,lon,time data if (lat_data==NULL) lat_data=(float*)read_data(lat_data,latitude_name.c_str()); if (lon_data==NULL) lon_data=(float*)read_data(lon_data,longitude_name.c_str()); if (time_data==NULL) { double nsec_tai93_day_start; // get TAI to UTC leap seconds double leap_sec; // TAI93 = epoch UTC - dtime double dtime; int sz = lat_lon_index_max + 1; Date d; float32* _time_data = NULL; switch ( product_id ) { case DARDAR_MASK : time_data=(double*)read_data(time_data, time_name.c_str()); break; case DARDAR_CLOUD : // --- UTC to TAI93 conversion _time_data = (float32*) read_data ( _time_data, time_name.c_str() ); // get number of seconds from start of day in TAI93 convention d.set_date ( date->get_year (), date->get_month (), date->get_day () ); nsec_tai93_day_start = d.get_TAI93_time(); // get TAI to UTC leap seconds leap_sec = Date::get_tai93_leap_sec ( d ); // TAI93 = epoch UTC - dtime dtime = nsec_tai93_day_start + leap_sec; time_data = new double [sz]; for ( int itime = 0 ; itime < sz ; ++itime ) time_data [itime] = _time_data [itime] + dtime; break; default : string msg ( "Unsupported product " + product ) ; g_exception e( __FILE__ , __LINE__ , msg ); throw e; break; } delete[] _time_data; } load_v_pixel(); if (!hdf_file_already_loaded) free_hdf_file(); } } const bool DARDARFileData::contain_location(const float &lat,const float &lon, const double &tolerance){ int nearest_idx=-1; return get_index(lat,lon,nearest_idx,colocation_tolerance); } const bool DARDARFileData::contain_data( const float & lat, const float & lon, const double & time, const double &colocation_tolerance ) { return contain_location(lat,lon,colocation_tolerance)&&contain_time(time); } const bool DARDARFileData::get_index( const float & lat, const float & lon, int &nearest_pix_idx, const float colocation_tolerance ) { assert(lat_data!=NULL && lon_data!=NULL && time_data!=NULL); bool already_loaded_data=is_geolocation_data_loaded(); if (!already_loaded_data) load_geolocation_data(); double delta_lat,delta_lon,delta,best_delta; delta_lat=delta_lon=delta=best_delta=9999.; // 9999 is bigger than the biggest delta that can be int _nearest_pix_idx=-1; /*** use a generic search on sorted datas ***/ // bounds of the colocation frame float lat_min=lat-colocation_tolerance; // lower latitude acceptable float lon_min=lon-colocation_tolerance; // lower longitude acceptable PixelType pix_min(&lat_min,&lon_min,vector(0)); float lat_max=lat+colocation_tolerance; // biggest latitude acceptable float lon_max=lon+colocation_tolerance; // biggest longitude acceptable PixelType pix_max(&lat_max,&lon_max,vector(0)); // cout << "----- (" << lat << ", " << lon << ") -----" << endl; // set the range of colocated points candidates I_Pixel i_pix_min=lower_bound(v_pixel.begin(),v_pixel.end(),pix_min); I_Pixel i_pix_max=upper_bound(v_pixel.begin(),v_pixel.end(),pix_max); // cout<<"Pix min ["<get_val()[0]<<"]"<<" ("<<(float)(i_pix_min->get_lat())<<","<<(float)(i_pix_min->get_lon())<<")"<get_val()[0]<<"]"<<" ("<get_lat()<<","<get_lon()<<")"<get_lat()-lat); delta_lon=abs(i_pix->get_lon()-lon); // cout<<"Possible point ("<get_lat()<<","<get_lon()<<") ["<get_val()[0]<<"]"<get_val()[0]; } } } if (!already_loaded_data) free_geolocation_data(); nearest_pix_idx=_nearest_pix_idx; return (nearest_pix_idx!=-1); } const float DARDARFileData::get_nearest_point_distance( const float & lat, const float & lon, const float coloc_tolerance ) { bool data_already_loaded=is_geolocation_data_loaded(); if (!data_already_loaded) load_geolocation_data(); int nearest_point_idx; double nearest_point_distance=-1.; if (get_index(lat,lon,nearest_point_idx,coloc_tolerance)) // a point in the coincidence frame has been found nearest_point_distance=sqrt(pow(lat_data[nearest_point_idx]-lat,2) + pow(lon_data[nearest_point_idx]-lon,2)); if (!data_already_loaded) free_geolocation_data(); return nearest_point_distance; } void DARDARFileData::set_lat_lon_index_max() { // read the dimension of the latitude VData. Assume that longitude and time ones have the same number of values vector v_dim(0); v_dim = get_sds_dimension(latitude_name); lat_lon_index_max=v_dim[0]-1; } void DARDARFileData::load_v_pixel() { //-- Build the sorted pixels list used for finding coincidences int record_max=lat_lon_index_max+1; // nb track records along the orbit vector ipix(1); // index of a pixel for (int record_idx=0;record_idx