/* reproject.cpp */ /* remap Copyright (C) 2006 Fabrice Ducos, fabrice.ducos@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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ #include #include #include "reproject.h" #include "common.h" #include "Pixel/Pixel.h" #include "hdffiledata.h" typedef pair rowcol_type; typedef struct { } empty_type; void reproject_coordinates(grid_type *src_grid, grid_type *target_grid) { assert(src_grid); assert(target_grid); using namespace std; int src_nrows = src_grid->nrows; int src_ncols = src_grid->ncols; int target_nrows = target_grid->nrows; int target_ncols = target_grid->ncols; coord_type *src_lat = src_grid->lat; coord_type *src_lon = src_grid->lon; coord_type *target_lat = target_grid->lat; coord_type *target_lon = target_grid->lon; int *src_irows = target_grid->src_irows; int *src_icols = target_grid->src_icols; typedef Pixel_base Pixel; multiset src_pixels; // indexes each source pixel to prepare reprojection (indexation is performed in the constructor of Pixel) Pixel::set_resolution(g_distance); clock_t t0 = clock(); for (int src_irow = 0 ; src_irow < src_nrows ; src_irow++) { for (int src_icol = 0 ; src_icol < src_ncols ; src_icol++) { const int isrc = src_irow * src_ncols + src_icol; const coord_type lat = src_lat[isrc]; if ((MIN_LAT <= lat && lat <= MAX_LAT) == false) continue; const coord_type lon = src_lon[isrc]; if ((MIN_LON <= lon && lon <= MAX_LON) == false) continue; const Pixel src_pixel(lat, lon, make_pair(src_irow, src_icol), Pixel::DEGREES); src_pixels.insert(src_pixel); } // src_icol } // src_irow clock_t t1 = clock(); double t01 = (t1 - t0)/static_cast(CLOCKS_PER_SEC); if (g_verbose == true) { std::clog << APPNAME << ": indexation time: " << t01 << " s" << std::endl; } // for each pixel in the target grid, finds the nearest source pixel and stores its measure in the target pixel const int npixels = target_nrows*target_ncols; const int nalert = npixels / 50; // will display a dot every 2% of processing (verbose mode only) int ipixel = 0; if (g_verbose == true) { std::clog << APPNAME << ": "; } clock_t t2 = clock(); for (int target_irow = 0 ; target_irow < target_nrows; target_irow++) { for (int target_icol = 0 ; target_icol < target_ncols ; target_icol++) { vector neighbours; const int itarget = target_irow * target_ncols + target_icol; src_irows[itarget] = -1; src_icols[itarget] = -1; if (g_verbose == true) { ipixel++; if (ipixel % nalert == 0) clog << "."; // alert (display a dot) every 'nalert' pixels } const coord_type lat = target_lat[itarget]; if ((MIN_LAT <= lat && lat <= MAX_LAT) == false) continue; const coord_type lon = target_lon[itarget]; if ((MIN_LON <= lon && lon <= MAX_LON) == false) continue; const Pixel target_pixel(lat, lon, make_pair(target_irow, target_icol), Pixel::DEGREES); target_pixel.get_neighbours(neighbours, src_pixels, g_distance); if (neighbours.size() > 0) { Pixel nearest = *min_element(neighbours.begin(), neighbours.end(), Nearer_from(target_pixel)); const rowcol_type src_rowcol = nearest.val(); const rowcol_type target_rowcol = nearest.val(); const int src_irow = src_rowcol.first; const int src_icol = src_rowcol.second; src_irows[itarget] = src_irow; src_icols[itarget] = src_icol; } } // target_icol } // target_irow if (g_verbose == true) { clock_t t3 = clock(); double t23 = (t3 - t2)/static_cast(CLOCKS_PER_SEC); std::clog << endl; std::clog << APPNAME << ": reprojection time: " << t23 << " s" << std::endl; } } /* checks if one tries to reproject into an already existing dataset in the output file * in that case, reloads the existing data from the output file * in the other case, instructs the code to initialize the dataset to a fill value. */ void initialize_target_grid(grid_type *src_grid, grid_type *target_grid) { int target_nrows = target_grid->nrows; int target_ncols = target_grid->ncols; data_type *target_data = target_grid->data; time_type *target_time_from_ref = target_grid->time_from_ref; distance_type *target_distance_from_ref = target_grid->distance_from_ref; bool measures_were_initialized = false; strncpy(target_grid->output_dataset, src_grid->output_dataset, STRING_MAXLEN); target_grid->ichannel = src_grid->ichannel; try { float64 slope, offset; float64 slope_err, offset_err; // unused int32 sds_type; // unused HDFFileData f_data(target_grid->file, "r"); vector data_dimensions = f_data.get_sds_dimension(target_grid->output_dataset); assert(data_dimensions[0] == target_grid->nrows); assert(data_dimensions[1] == target_grid->ncols); // target_grid->output_dataset already exists in target_grid->file, load its data in target_grid->data target_grid->data = f_data.read_data(target_grid->data, target_grid->output_dataset); f_data.get_hdf_file()->get_calibration(target_grid->output_dataset, slope, slope_err, offset, offset_err, sds_type); target_grid->slope = slope; target_grid->offset = offset; measures_were_initialized = true; if (g_extended_output) { char sds_distance_from_ref[STRING_MAXLEN + 1]; char sds_time_from_ref[STRING_MAXLEN + 1]; snprintf(sds_distance_from_ref, STRING_MAXLEN, "%s{distance_from_ref}", target_grid->output_dataset); snprintf(sds_time_from_ref, STRING_MAXLEN, "%s{time_from_ref}", target_grid->output_dataset); HDFFileData f_data(target_grid->file, "r"); target_grid->distance_from_ref = static_cast(f_data.read_data(static_cast(target_grid->distance_from_ref), sds_distance_from_ref)); target_grid->time_from_ref = static_cast(f_data.read_data(static_cast(target_grid->time_from_ref), sds_time_from_ref)); } } catch (bad_sds_name &e) { //assert(measures_were_initialized == false); /* target_grid->output_dataset doesn't exist yet, do nothing */ } catch (bad_file &e) { Debug(cerr << APPNAME << ": " << __PRETTY_FUNCTION__ << ": " << e.what() << endl;); throw; } catch (...) { cerr << APPNAME << ": " __FILE__ << ": " << __LINE__ << ": unexpected exception" << endl; throw; } if (measures_were_initialized == false) { target_grid->slope = src_grid->slope; target_grid->offset = src_grid->offset; for (int target_irow = 0 ; target_irow < target_nrows; target_irow++) { for (int target_icol = 0 ; target_icol < target_ncols ; target_icol++) { const int itarget = target_irow*target_ncols + target_icol; target_data[itarget] = g_data_fill_value; if (g_extended_output) { target_time_from_ref[itarget] = g_time_fill_value; target_distance_from_ref[itarget] = g_distance_fill_value; } } } } } /* reprojects measures without recomputing nearest neighbours coordinates (uses previous results from reproject_coordinates) */ void reproject_measures(grid_type *src_grid, grid_type *target_grid) { assert(src_grid); assert(target_grid); using namespace std; typedef Pixel_base Pixel; int src_nrows = src_grid->nrows; int src_ncols = src_grid->ncols; int target_nrows = target_grid->nrows; int target_ncols = target_grid->ncols; data_type *src_data = src_grid->data; time_type *src_tim = src_grid->tim; data_type *target_data = target_grid->data; time_type *target_tim = target_grid->tim; time_type *target_time_from_ref = target_grid->time_from_ref; distance_type *target_distance_from_ref = target_grid->distance_from_ref; const coord_type *src_lat = src_grid->lat; const coord_type *src_lon = src_grid->lon; const coord_type *target_lat = target_grid->lat; const coord_type *target_lon = target_grid->lon; int *src_irows = target_grid->src_irows; int *src_icols = target_grid->src_icols; // for each pixel in the target grid, finds the nearest source pixel and stores its measure in the target pixel initialize_target_grid(src_grid, target_grid); for (int target_irow = 0 ; target_irow < target_nrows; target_irow++) { for (int target_icol = 0 ; target_icol < target_ncols ; target_icol++) { const int itarget = target_irow*target_ncols + target_icol; if ((MIN_LAT <= target_lat[itarget] && target_lat[itarget] <= MAX_LAT) == false) continue; if ((MIN_LON <= target_lon[itarget] && target_lon[itarget] <= MAX_LON) == false) continue; const int src_irow = src_irows[itarget]; const int src_icol = src_icols[itarget]; const int isrc = src_irow * src_ncols + src_icol; if ((MIN_LAT <= src_lat[isrc] && src_lat[isrc] <= MAX_LAT) == false) continue; if ((MIN_LON <= src_lon[isrc] && src_lon[isrc] <= MAX_LON) == false) continue; if ((0 <= src_irow && src_irow < src_nrows) == false || (0 <= src_icol && src_icol < src_ncols) == false) continue; const time_type src_time = src_tim[isrc]; const time_type target_time = target_tim[itarget]; const time_type dtime = fabs(target_time - src_time); if ((g_max_dtime < INFINITY) && (dtime <= g_max_dtime) == false) continue; target_data[itarget] = src_data[isrc]; if (g_extended_output) { empty_type empty; const Pixel src_pixel(src_lat[isrc], src_lon[isrc], empty, Pixel::DEGREES); const Pixel target_pixel(target_lat[itarget], target_lon[itarget], empty, Pixel::DEGREES); target_time_from_ref[itarget] = dtime; target_distance_from_ref[itarget] = target_pixel.distance(src_pixel); } } // target_icol } // target_irow } void rescale(grid_type *grid, float64 new_slope, float64 new_offset, data_type old_undef_count, data_type new_undef_count) { PDEBUG; int nrows = grid->nrows; int ncols = grid->ncols; float64 old_slope = grid->slope; float64 old_offset = grid->offset; data_type *count = grid->data; if (new_slope == old_slope && new_offset == old_offset) { return; // rescaling is useless } for (int irow = 0 ; irow < nrows ; ++irow) { for (int icol = 0 ; icol < ncols ; ++icol) { data_type old_count = count[irow*ncols + icol]; data_type new_count; if (old_count != old_undef_count) { double calibrated_value = old_slope*(old_count - old_offset); new_count = lround(calibrated_value/new_slope + new_offset); } else { new_count = new_undef_count; } count[irow*ncols + icol] = new_count; } } grid->slope = new_slope; grid->offset = new_offset; }