/* grid.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 "common.h" #include "filetypes.h" #include "grid.h" #include "hdf_utils.h" #include "reproject.h" #include "VInputs.h" void reset_grid(grid_type *grid) { PDEBUG; assert(grid); int nrows = grid->nrows; int ncols = grid->ncols; for (int irow = 0 ; irow < nrows ; irow++) { for (int icol = 0 ; icol < ncols ; icol++) { grid->lat[irow*ncols + icol] = g_coord_fill_value; grid->lon[irow*ncols + icol] = g_coord_fill_value; grid->tim[irow*ncols + icol] = g_time_fill_value; grid->data[irow*ncols + icol] = g_data_fill_value; if (g_extended_output) { assert(grid->distance_from_ref); assert(grid->time_from_ref); grid->distance_from_ref[irow*ncols + icol] = g_distance_fill_value; grid->time_from_ref[irow*ncols + icol] = g_time_fill_value; } } } grid->first_lat = g_coord_fill_value; grid->first_lon = g_coord_fill_value; grid->last_lat = g_coord_fill_value; grid->last_lon = g_coord_fill_value; if (grid->is_target) { assert(grid->src_irows); assert(grid->src_icols); for (int irow = 0 ; irow < nrows ; irow++) { for (int icol = 0 ; icol < ncols ; icol++) { grid->src_irows[irow*ncols + icol] = irow; grid->src_icols[irow*ncols + icol] = icol; } } } } void create_grid(grid_type *grid, int nrows, int ncols, float64 slope, float64 offset) { PDEBUG; assert(grid); grid->nrows = nrows; grid->ncols = ncols; grid->first_lat = g_coord_fill_value; grid->first_lon = g_coord_fill_value; grid->last_lat = g_coord_fill_value; grid->last_lon = g_coord_fill_value; grid->lat = new coord_type[nrows*ncols]; assert(grid->lat); grid->lon = new coord_type[nrows*ncols]; assert(grid->lon); grid->tim = new time_type[nrows*ncols]; assert(grid->tim); grid->data = new data_type[nrows*ncols]; assert(grid->data); if (g_extended_output) { grid->distance_from_ref = new distance_type[nrows*ncols]; assert(grid->distance_from_ref); grid->time_from_ref = new time_type[nrows*ncols]; assert(grid->time_from_ref); } grid->slope = slope; grid->offset = offset; if (grid->is_target) { grid->src_irows = new int[nrows*ncols]; assert(grid->src_irows); grid->src_icols = new int[nrows*ncols]; assert(grid->src_icols); } assert(grid->data); reset_grid(grid); } void destroy_grid(grid_type *grid) { PDEBUG; assert(grid); if (grid->lat) { delete[] grid->lat; grid->lat = NULL; } if (grid->lon) { delete[] grid->lon; grid->lon = NULL; } if (grid->tim) { delete[] grid->tim; grid->tim = NULL; } if (grid->data) { delete[] grid->data; grid->data = NULL; } if (g_extended_output) { if (grid->distance_from_ref) { delete[] grid->distance_from_ref; grid->distance_from_ref = NULL; } if (grid->time_from_ref) { delete[] grid->time_from_ref; grid->time_from_ref = NULL; } } if (grid->src_irows) { delete[] grid->src_irows; grid->src_irows = NULL; } if (grid->src_icols) { delete[] grid->src_icols; grid->src_icols = NULL; } memset(&grid, 0, sizeof(grid)); } // looks for the first and last valid pixel coordinates in the grid static void find_grid_extremities(grid_type *grid) { PDEBUG; assert(grid); const int nrows = grid->nrows; const int ncols = grid->ncols; const int npixels = nrows*ncols; int first_ipixel = 0; int last_ipixel = npixels - 1; coord_type first_lat, first_lon; // first valid coordinates in the grid coord_type last_lat, last_lon; // last valid coordinates in the grid do { first_lat = grid->lat[first_ipixel]; first_lon = grid->lon[first_ipixel]; } while ((++first_ipixel < npixels) && ((MIN_LAT <= first_lat && first_lat <= MAX_LAT) == false || (MIN_LON <= first_lon && first_lon <= MAX_LON) == false)); do { last_lat = grid->lat[last_ipixel]; last_lon = grid->lon[last_ipixel]; } while ((last_ipixel-- > 0) && ((MIN_LAT <= last_lat && last_lat <= MAX_LAT) == false || (MIN_LON <= last_lon && last_lon <= MAX_LON) == false)); assert(first_ipixel <= last_ipixel); grid->first_lat = first_lat; grid->first_lon = first_lon; grid->last_lat = last_lat; grid->last_lon = last_lon; } static void load_grid_helper(grid_type *grid, VInput *vinput) { PDEBUG; assert(grid); assert(vinput); int nrows; int ncols; vector dimensions = vinput->dimensions(); grid->rank = 2; assert(dimensions.size() == static_cast(grid->rank)); grid->nrows = nrows = dimensions[0]; grid->ncols = ncols = dimensions[1]; if (dynamic_cast(vinput)) { VFile *vfile = static_cast(vinput); try { vfile->get_calibration(grid->slope, grid->offset); } catch (...) { ostringstream s; s << APPNAME << ": " << vfile->filename() << ": " << vfile->dataset() << ": error while reading calibration factors" << endl; throw s.str().c_str(); } } else { grid->slope = 1.; grid->offset = 0.; } assert(nrows*ncols >= 1); create_grid(grid, nrows, ncols, grid->slope, grid->offset); grid->content = vinput->content(); grid->wavelength = vinput->wavelength(); grid->data_type_code = vinput->data_type_code(); memcpy(grid->lat, vinput->lat(), nrows*ncols*sizeof(*(grid->lat))); memcpy(grid->lon, vinput->lon(), nrows*ncols*sizeof(*(grid->lon))); memcpy(grid->tim, vinput->time(), nrows*ncols*sizeof(*(grid->tim))); if (vinput->data()) { coord_type *grid_lat = grid->lat; coord_type *grid_lon = grid->lon; data_type *grid_data = grid->data; data_type *vinput_data = static_cast(vinput->data()); for (int irow = 0 ; irow < nrows ; ++irow) { for (int icol = 0 ; icol < ncols ; ++icol) { const int ipixel = irow*ncols + icol; // we set all invalid values to g_data_fill_value const data_type pixel_data = vinput_data[ipixel]; if ((MIN_LAT <= grid_lat[ipixel] && grid_lat[ipixel] <= MAX_LAT) == false || (MIN_LON <= grid_lon[ipixel] && grid_lon[ipixel] <= MAX_LON) == false || (MIN_VALID_DATA <= pixel_data && pixel_data <= MAX_VALID_DATA) == false) { grid_data[ipixel] = g_data_fill_value; } else { grid_data[ipixel] = pixel_data; } } } rescale(grid, g_output_slope, g_output_offset, g_data_fill_value, g_data_fill_value); } find_grid_extremities(grid); } int load_grid(grid_type *grid) { using namespace std; PDEBUG; assert(grid); char *file = grid->file; char *dataset = grid->input_dataset; int ichannel = grid->ichannel; filetype_type filetype; grid->slope = 1.0; grid->offset = 0.0; if (file[0] == '\0') { filetype = FILETYPE_EQUIRECT; } else { filetype = get_filetype(file); } try { switch (filetype) { case FILETYPE_EQUIRECT : { assert(MIN_LAT <= g_equirect_latitude_orig && g_equirect_latitude_orig <= MAX_LAT); assert(MIN_LON <= g_equirect_longitude_orig && g_equirect_longitude_orig <= MAX_LON); assert(g_equirect_resolution > 0.); assert(g_equirect_nrows > 0); assert(g_equirect_ncols > 0); VEquiRect vinput(g_equirect_latitude_orig, g_equirect_longitude_orig, g_equirect_resolution, g_equirect_nrows, g_equirect_ncols); load_grid_helper(grid, &vinput); break; } case FILETYPE_MOD021KM : case FILETYPE_MYD021KM : { VModis vinput(file, dataset, ichannel); load_grid_helper(grid, &vinput); break; } case FILETYPE_MOD03 : case FILETYPE_MYD03 : { VMOD03 vinput(file, dataset, ichannel); load_grid_helper(grid, &vinput); break; } case FILETYPE_CAL_IIR_L1 : { VIIR vinput(file, dataset, ichannel); load_grid_helper(grid, &vinput); break; } case FILETYPE_HDF_SEVIRI : { assert(ichannel == 0); VHdf_Seviri vinput(file, dataset); load_grid_helper(grid, &vinput); break; } case FILETYPE_XRIT_SEVIRI : { assert(ichannel == 0); VXRIT_SEVIRI vinput(file, dataset); load_grid_helper(grid, &vinput); break; } case FILETYPE_PARASOL : { assert(ichannel == 0); VParasol vinput(file, dataset); load_grid_helper(grid, &vinput); break; } case FILETYPE_UNKNOWN : { cerr << APPNAME ": " << file << ": unknown or not yet implemented format" << endl; return -1; break; } default : { cerr << APPNAME ": internal error: unexpected filetype value: " << filetype << endl; exit (EXIT_FAILURE); break; } } // switch (filetype) } catch (const char *err_msg) { cerr << APPNAME ": " << err_msg << endl; return -1; } catch (const string &err_msg) { cerr << APPNAME ": " << err_msg << endl; return -1; } catch (exception &e) { cerr << APPNAME ": " << e.what() << endl; return -1; } return 0; } int save_latlontime(grid_type *grid) { using namespace std; PDEBUG; assert(grid); assert(grid->lat); assert(grid->lon); assert(grid->tim); int32 rank = 2; int32 start[2] = { 0, 0 }; int32 edges[2]; int32 *stride = NULL; int err; edges[0] = grid->nrows; edges[1] = grid->ncols; err = hdf_add_empty_sds(grid->file, "Latitude" , HDF_COORD_TYPE, rank, edges, &g_coord_fill_value); err = hdf_add_empty_sds(grid->file, "Longitude", HDF_COORD_TYPE, rank, edges, &g_coord_fill_value); err = hdf_add_empty_sds(grid->file, "Time_TAI93", HDF_TIME_TYPE, rank, edges, &g_time_fill_value); try { Hdf_file hdf_file(grid->file, DFACC_WRITE); hdf_file.write_sds("Latitude", grid->lat, start, stride, edges); hdf_file.write_sds("Longitude", grid->lon, start, stride, edges); hdf_file.write_sds("Time_TAI93", grid->tim, start, stride, edges); } catch (string &err_msg) { cerr << APPNAME ": " << grid->file << ": " << err_msg << endl; return -1; } catch (sd_write_data_error &e) { cerr << APPNAME ": " << grid->file << ": " << e.what() << endl; return -1; } return 0; } int save_grid(grid_type *grid, bool save_delta) { using namespace std; PDEBUG; assert(grid); int32 rank = 2; int32 start[2] = { 0, 0 }; int32 edges[2]; int32 *stride = NULL; int err; char sds_distance_from_ref[STRING_MAXLEN + 1]; char sds_time_from_ref[STRING_MAXLEN + 1]; memset(sds_distance_from_ref, 0, sizeof(sds_distance_from_ref)); memset(sds_time_from_ref, 0, sizeof(sds_time_from_ref)); edges[0] = grid->nrows; edges[1] = grid->ncols; err = hdf_add_empty_sds(grid->file, grid->output_dataset, HDF_DATA_TYPE, rank, edges, &g_data_fill_value, grid->slope, grid->offset); if (save_delta) { snprintf(sds_distance_from_ref, STRING_MAXLEN, "%s{distance_from_ref}", grid->output_dataset); snprintf(sds_time_from_ref, STRING_MAXLEN, "%s{time_from_ref}", grid->output_dataset); err = hdf_add_empty_sds(grid->file, sds_distance_from_ref, HDF_DISTANCE_TYPE, rank, edges, &g_distance_fill_value); err = hdf_add_empty_sds(grid->file, sds_time_from_ref, HDF_TIME_TYPE, rank, edges, &g_time_fill_value); } try { assert(grid->data); Hdf_file hdf_file(grid->file, DFACC_WRITE); hdf_file.write_sds(grid->output_dataset, grid->data, start, stride, edges); if (save_delta) { assert(grid->distance_from_ref); assert(grid->time_from_ref); hdf_file.write_sds(sds_distance_from_ref, grid->distance_from_ref, start, stride, edges); hdf_file.write_sds(sds_time_from_ref, grid->time_from_ref, start, stride, edges); } } catch (string &err_msg) { cerr << APPNAME ": " << grid->file << ": " << err_msg << endl; return -1; } catch (sd_write_data_error &e) { cerr << APPNAME ": " << grid->file << ": " << e.what() << endl; return -1; } return 0; } bool have_same_grid_coordinates(grid_type *grid1, grid_type *grid2) { PDEBUG; assert(grid1); assert(grid2); /* bool ret = (grid1->nrows == grid2->nrows && grid1->ncols == grid2->ncols && grid1->first_lat == grid2->first_lat && grid1->first_lon == grid2->first_lon && grid1->last_lat == grid2->last_lat && grid1->last_lon == grid2->last_lon); */ const coord_type epsilon = 0.1; bool ret = (grid1->nrows == grid2->nrows && grid1->ncols == grid2->ncols && fabs(grid1->first_lat - grid2->first_lat) < epsilon && fabs(grid1->first_lon - grid2->first_lon) < epsilon && fabs(grid1->last_lat - grid2->last_lat) < epsilon && fabs(grid1->last_lon - grid2->last_lon) < epsilon); return ret; } bool have_same_grid_footprint(grid_type *grid1, grid_type *grid2) { PDEBUG; assert(grid1); assert(grid2); if (grid1->ichannel != grid2->ichannel) return false; if (grid1->rank != grid2->rank) return false; if (have_same_grid_coordinates(grid1, grid2) == false) { return false; } if (grid1->slope != grid2->slope) return false; if (grid1->offset != grid2->offset) return false; if (grid1->content != grid2->content) return false; if (grid1->wavelength != grid2->wavelength) return false; // string comparisons are the most expensive operations, so we do them at last if (strncmp(grid1->file, grid2->file, STRING_MAXLEN) != 0) return false; if (strncmp(grid1->input_dataset, grid2->input_dataset, STRING_MAXLEN) != 0) return false; if (strncmp(grid1->output_dataset, grid2->output_dataset, STRING_MAXLEN) != 0) return false; return true; } void copy_grid_footprint(grid_type *target_grid, grid_type *src_grid) { PDEBUG; assert(src_grid); assert(target_grid); strncpy(target_grid->file, src_grid->file, STRING_MAXLEN); strncpy(target_grid->input_dataset, src_grid->input_dataset, STRING_MAXLEN); strncpy(target_grid->output_dataset, src_grid->output_dataset, STRING_MAXLEN); target_grid->nrows = src_grid->nrows; target_grid->ncols = src_grid->ncols; target_grid->ichannel = src_grid->ichannel; target_grid->rank = src_grid->rank; target_grid->first_lat = src_grid->first_lat; target_grid->first_lon = src_grid->first_lon; target_grid->last_lat = src_grid->last_lat; target_grid->last_lon = src_grid->last_lon; target_grid->slope = src_grid->slope; target_grid->offset = src_grid->offset; target_grid->content = src_grid->content; target_grid->wavelength = src_grid->wavelength; assert(have_same_grid_footprint(target_grid, src_grid)); } void copy_grid(grid_type *target_grid, grid_type *src_grid) { PDEBUG; assert(src_grid); assert(target_grid); int nrows = src_grid->nrows; int ncols = src_grid->ncols; float64 slope = src_grid->slope; float64 offset = src_grid->offset; destroy_grid(target_grid); create_grid(target_grid, nrows, ncols, slope, offset); copy_grid_footprint(target_grid, src_grid); for (int irow = 0 ; irow < nrows ; irow++) { for (int icol = 0 ; icol < ncols ; icol++) { target_grid->lat[irow*ncols + icol] = src_grid->lat[irow*ncols + icol]; target_grid->lon[irow*ncols + icol] = src_grid->lon[irow*ncols + icol]; target_grid->tim[irow*ncols + icol] = src_grid->tim[irow*ncols + icol]; target_grid->data[irow*ncols + icol] = src_grid->data[irow*ncols + icol]; if (g_extended_output) { assert(src_grid->distance_from_ref); assert(src_grid->time_from_ref); assert(target_grid->distance_from_ref); assert(target_grid->time_from_ref); target_grid->distance_from_ref[irow*ncols + icol] = src_grid->distance_from_ref[irow*ncols + icol]; target_grid->time_from_ref[irow*ncols + icol] = src_grid->time_from_ref[irow*ncols + icol]; } } } if (src_grid->is_target) { assert(target_grid->is_target); assert(src_grid->src_irows); assert(src_grid->src_icols); assert(target_grid->src_irows); assert(target_grid->src_icols); for (int irow = 0 ; irow < nrows ; irow++) { for (int icol = 0 ; icol < ncols ; icol++) { target_grid->src_irows[irow*ncols + icol] = src_grid->src_irows[irow*ncols + icol]; target_grid->src_icols[irow*ncols + icol] = src_grid->src_icols[irow*ncols + icol]; } } } } void grid_print(const grid_type *grid) { assert(grid); int nrows = grid->nrows; int ncols = grid->ncols; coord_type *lat = grid->lat; coord_type *lon = grid->lon; data_type *data = grid->data; for (int irow = 0 ; irow < nrows ; ++irow) { cout << irow << ": " << "\t" << lat[irow*ncols] << "\t" << lon[irow*ncols] << "\t" << data[irow*ncols] << endl; } }