/*************************************************************************** * 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 "cloudsatfiledata.h" const string CLOUDSATFileData::time_name="Profile_time"; const string CLOUDSATFileData::time_origin_name="TAI_start"; const float CLOUDSATFileData::colocation_tolerance=0.01; CLOUDSATFileData::CLOUDSATFileData(const string &name /*= ""*/, const string &mode/*= "r"*/) : FileData(name,mode),SatelliteFileData(name,mode),HDFFileData(name,mode) { init(); } void CLOUDSATFileData::init() { lat_data=NULL; lon_data=NULL; time_data=NULL; product=""; 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(); } CLOUDSATFileData::~CLOUDSATFileData() { free_geolocation_data(); free_hdf_metadata(); free_hdf_file(); } void CLOUDSATFileData::set_time_coverage( ) { //-- Read the end time in the hdf file attributes char schar_end_time[14+1]; // date has format "YYYYMMDDHHMMSS" -> 14 characters memset(schar_end_time,'\0',15); hdf_file->get_attr("end_time").get_values(schar_end_time); string s_endtime(schar_end_time); //-- Build a date with this end time Date end_time; end_time.set_date_str(s_endtime,"%Y%m%d%H%M%S"); //-- Substrack end and start time to compute the time coverage time_coverage=end_time.get_epoch_time()-date->get_epoch_time(); } void CLOUDSATFileData::free_geolocation_data() { delete[] lat_data; lat_data=NULL; delete[] lon_data; lon_data=NULL; delete[] time_data; time_data=NULL; } void CLOUDSATFileData::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) { if ( product != "MODIS-AUX" ) // stored in a VDATA lat_data=(float*)read_vdata(lat_data,latitude_name.c_str(),latitude_name.c_str()); else // stored in a SDS lat_data=(float*)read_data(lat_data,latitude_name.c_str()); } if (lon_data==NULL) { if ( product != "MODIS-AUX" ) // stored in a VDATA lon_data=(float*)read_vdata(lon_data,longitude_name.c_str(),longitude_name.c_str()); else // stored in a SDS lon_data=(float*)read_data(lon_data,longitude_name.c_str()); } if (time_data==NULL) { // read the timestamp of the first profile double time_start=0.; read_vdata(&time_start,time_origin_name.c_str(),time_origin_name.c_str()); // read the times passed after the first profile float32* tmp_time_data=NULL; try{ tmp_time_data=(float32*)read_vdata(time_data,time_name.c_str(),time_name.c_str()); // set the profiles times in TAI93 convention time_data=new double[lat_lon_index_max+1]; for (int i=0;i<=lat_lon_index_max;++i) time_data[i]=tmp_time_data[i]+time_start; delete[] tmp_time_data; } catch (exception &e) { delete[] tmp_time_data; delete[] time_data, time_data=NULL; throw; } } load_v_pixel(); if (!hdf_file_already_loaded) free_hdf_file(); } } void CLOUDSATFileData::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 < int > ipix(1); // index of a pixel for (int record_idx=0;record_idx v_dim(0); if ( product != "MODIS-AUX" ) // stored in a VDATA v_dim = get_vdata_dimension(latitude_name); else // stored in a SDS v_dim = get_sds_dimension(latitude_name); lat_lon_index_max=v_dim[0]-1; } void CLOUDSATFileData::parse_filename( const string& short_filename ) { check_filename(short_filename); // throw a parse_CLOUDSAT_filename_error if the filename doesn't matches the pattern int infilename_size=short_filename.size(); using namespace MyTools; // for to_num function // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time) date->set_date_str(short_filename.substr(0,13),"%Y%j%H%M%S"); // orbit number orbit_nb=to_num(short_filename.substr(14,5)); // product name size_t underscore_pos=short_filename.find_first_of("_",23); product=short_filename.substr(23,underscore_pos-23); // geolocation SDS names depend of the product if ( product == "MODIS-AUX" ) { latitude_name="MODIS_latitude"; longitude_name="MODIS_longitude"; } else { latitude_name="Latitude"; longitude_name="Longitude"; } // Version : Reprocessing number+Epoch number version=short_filename.substr(infilename_size-11,7); } bool CLOUDSATFileData::check_filename( const string & short_filename ) { int infilename_size=short_filename.size(); using namespace MyTools; // for has_type function // timestamp : the first 14 characters are digits and represents the orbit start date (in UTC time) string s_date=short_filename.substr(0,13); if (!has_type(s_date)) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } // orbit number string s_orbit_nb=short_filename.substr(14,5); if (!has_type(s_orbit_nb)) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } // Mission string s_mission=short_filename.substr(20,2); if (s_mission!=string("CS")) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } // product name size_t underscore_pos=short_filename.find_first_of("_",23); string s_product=short_filename.substr(23,underscore_pos-23); if (s_product!="1B-CPR" && s_product!="1B-CPR-FL" && s_product!="2B-GEOPROF" && s_product!="2B-GEOPROF-LIDAR" && s_product!="2B-LWC-RO" && s_product!="2B-IWC-RO" && s_product!="2B-CWC-RO" && s_product!="2B-LWC-RVOD" && s_product!="2B-IWC-RVOD" && s_product!="2B-CWC-RVOD" && s_product!="2B-TAU" && s_product!="2B-FLXHR" && s_product!="2B-CLDCLASS" && s_product!="2B-CLDCLASS-LIDAR" && s_product!="ECMWF-AUX" && s_product!="MODIS-AUX" ) { 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; } // Epoch number string s_epoch_nb=short_filename.substr(infilename_size-6,2); if (short_filename[infilename_size-7]!='E' || !has_type(s_epoch_nb)) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } // Reprocessing number string s_reprocessing_nb=short_filename.substr(infilename_size-10,2); if (short_filename[infilename_size-11]!='R' || !has_type(s_reprocessing_nb)) { invalid_filename e(__FILE__,__LINE__,short_filename,"CLOUDSAT"); throw e; } return true; } const bool CLOUDSATFileData::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 CLOUDSATFileData::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 CLOUDSATFileData::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 CLOUDSATFileData::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; }