Commit 1f8ae515 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapReadDataScrip: refactor.

parent b7424920
......@@ -24,7 +24,9 @@
#endif
#include "cdo_int.h"
#include "grid.h"
#include "remap_vars.h"
#include "remap.h" // SubmapType
#ifdef HAVE_LIBNETCDF
static void
......@@ -34,14 +36,396 @@ nce(int istat)
if (istat != NC_NOERR) cdoAbort(nc_strerror(istat));
}
void
nc_read_var_double(int nc_file_id, const char *name, double *array)
{
int ncvarid;
nce(nc_inq_varid(nc_file_id, name, &ncvarid));
nce(nc_get_var_double(nc_file_id, ncvarid, array));
}
void
nc_read_coordinate_radian(int nc_file_id, const char *name, size_t size, double *array)
{
int ncvarid;
nce(nc_inq_varid(nc_file_id, name, &ncvarid));
nce(nc_get_var_double(nc_file_id, ncvarid, array));
size_t attlen;
char grid_units[64];
nce(nc_get_att_text(nc_file_id, ncvarid, "units", grid_units));
nce(nc_inq_attlen(nc_file_id, ncvarid, "units", &attlen));
grid_units[attlen] = 0;
grid_to_radian(grid_units, size, array, name);
}
struct RemapGridW
{
int gridID;
int rank; // rank of the grid
size_t size; // total points on the grid
size_t num_cell_corners; // number of corners for each grid cell
bool lneed_cell_corners;
bool luse_cell_corners; // use corners for bounding boxes
bool lextrapolate;
bool non_global;
bool is_cyclic;
size_t dims[2]; // size of grid dimension
int nvgp; // size of vgpm
std::vector<int> vgpm; // flag which cells are valid
std::vector<int> mask; // flag which cells participate
std::vector<double> cell_center_lon; // lon/lat coordinates for
std::vector<double> cell_center_lat; // each grid center in radians
std::vector<double> cell_corner_lon; // lon/lat coordinates for
std::vector<double> cell_corner_lat; // each grid corner in radians
std::vector<double> cell_area; // total area of each grid cell
std::vector<double> cell_frac; // fractional area of grid cells participating in remapping
};
void
remapGridWInit(RemapGridW &grid)
{
grid.num_cell_corners = 0;
grid.luse_cell_corners = false;
grid.lneed_cell_corners = false;
grid.nvgp = 0;
}
void
remapGridWAlloc(RemapMethod mapType, RemapGridW &grid)
{
if (grid.nvgp) grid.vgpm.resize(grid.nvgp);
grid.mask.resize(grid.size);
grid.cell_center_lon.resize(grid.size);
grid.cell_center_lat.resize(grid.size);
if (mapType == RemapMethod::CONSERV_SCRIP || mapType == RemapMethod::CONSERV)
{
grid.cell_area.resize(grid.size, 0.0);
}
grid.cell_frac.resize(grid.size, 0.0);
if (grid.lneed_cell_corners)
{
if (grid.num_cell_corners > 0)
{
size_t nalloc = grid.num_cell_corners * grid.size;
grid.cell_corner_lon.resize(nalloc, 0.0);
grid.cell_corner_lat.resize(nalloc, 0.0);
}
}
}
struct RemapVarsW
{
bool sort_add;
bool pinit; // true: if the pointers are initialized
RemapMethod mapType; // identifier for remapping method
NormOpt normOpt; // option for normalization (conserv only)
long links_per_value;
size_t max_links; // current size of link arrays
size_t num_links; // actual number of links for remapping
size_t num_wts; // num of weights used in remapping
size_t resize_increment; // default amount to increase array size
std::vector<size_t> src_cell_add; // source grid address for each link
std::vector<size_t> tgt_cell_add; // target grid address for each link
std::vector<double> wts; // map weights for each link [max_links*num_wts]
};
void
remapVarsWInit(RemapMethod mapType, int remapOrder, RemapVarsW &rv)
{
// Initialize all pointer
if (!rv.pinit) rv.pinit = true;
rv.sort_add = (mapType == RemapMethod::CONSERV_SCRIP);
// Determine the number of weights
rv.num_wts = (mapType == RemapMethod::CONSERV_SCRIP) ? 3 : ((mapType == RemapMethod::BICUBIC) ? 4 : 1);
if (mapType == RemapMethod::CONSERV && remapOrder==2) rv.num_wts = 3;
rv.links_per_value = -1;
rv.num_links = 0;
rv.max_links = 0;
rv.resize_increment = 1024;
}
RemapMethod getMapType(int nc_file_id, SubmapType *submapType, int *numNeighbors, int &remapOrder);
void readLinks(int nc_file_id, int nc_add_id, size_t num_links, size_t *cell_add);
static void
remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *submapType,
int *numNeighbors, int &remapOrder, RemapGridW &src_grid, RemapGridW &tgt_grid, RemapVarsW &rv)
{
// The routine reads a NetCDF file to extract remapping info in SCRIP format
bool lgridarea = false;
int status;
int nc_dstgrddims_id; /* id for dest grid dimensions */
int nc_dstgrdimask_id; /* id for dest grid mask */
int nc_srcadd_id; /* id for map source address */
int nc_dstadd_id; /* id for map destination address */
int nc_rmpmatrix_id; /* id for remapping matrix */
size_t attlen, dimlen;
// Open file and read some global information
// nce(nc_open(interp_file, NC_NOWRITE, &nc_file_id));
int nc_file_id = cdf_openread(interp_file);
// Map name
char map_name[1024];
nce(nc_get_att_text(nc_file_id, NC_GLOBAL, "title", map_name));
nce(nc_inq_attlen(nc_file_id, NC_GLOBAL, "title", &attlen));
map_name[attlen] = 0;
if (Options::cdoVerbose)
{
cdoPrint("Reading remapping: %s", map_name);
cdoPrint("From file: %s", interp_file);
}
// Map Type
*mapType = getMapType(nc_file_id, submapType, numNeighbors, remapOrder);
if (*mapType == RemapMethod::CONSERV_SCRIP) lgridarea = true;
remapVarsWInit(*mapType, remapOrder, rv);
rv.mapType = *mapType;
rv.links_per_value = -1;
rv.sort_add = false;
// Normalization option
char normalize_opt[64]; // character string for normalization option
nce(nc_get_att_text(nc_file_id, NC_GLOBAL, "normalization", normalize_opt));
nce(nc_inq_attlen(nc_file_id, NC_GLOBAL, "normalization", &attlen));
normalize_opt[attlen] = 0;
// clang-format off
if (cstrIsEqual(normalize_opt, "none")) rv.normOpt = NormOpt::NONE;
else if (cstrIsEqual(normalize_opt, "fracarea")) rv.normOpt = NormOpt::FRACAREA;
else if (cstrIsEqual(normalize_opt, "destarea")) rv.normOpt = NormOpt::DESTAREA;
else
{
cdoPrint("normalize_opt = %s", normalize_opt);
cdoAbort("Invalid normalization option");
}
// clang-format on
if (Options::cdoVerbose) cdoPrint("normalize_opt = %s", normalize_opt);
// File convention
char convention[64]; // character string for output convention
nce(nc_get_att_text(nc_file_id, NC_GLOBAL, "conventions", convention));
nce(nc_inq_attlen(nc_file_id, NC_GLOBAL, "conventions", &attlen));
convention[attlen] = 0;
if (!cstrIsEqual(convention, "SCRIP"))
{
cdoPrint("convention = %s", convention);
if (cstrIsEqual(convention, "NCAR-CSM"))
cdoAbort("Unsupported file convention: %s!", convention);
else
cdoAbort("Unknown file convention!");
}
// Read some additional global attributes
// Source and destination grid names
char src_grid_name[64]; // grid name for source grid
nce(nc_get_att_text(nc_file_id, NC_GLOBAL, "source_grid", src_grid_name));
nce(nc_inq_attlen(nc_file_id, NC_GLOBAL, "source_grid", &attlen));
src_grid_name[attlen] = 0;
char tgt_grid_name[64]; // grid name for dest grid
nce(nc_get_att_text(nc_file_id, NC_GLOBAL, "dest_grid", tgt_grid_name));
nce(nc_inq_attlen(nc_file_id, NC_GLOBAL, "dest_grid", &attlen));
tgt_grid_name[attlen] = 0;
if (Options::cdoVerbose) cdoPrint("Remapping between: %s and %s", src_grid_name, tgt_grid_name);
// Initialize remapgrid structure
remapGridWInit(src_grid);
remapGridWInit(tgt_grid);
// Read dimension information
int ncdimid;
nce(nc_inq_dimid(nc_file_id, "src_grid_size", &ncdimid));
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
src_grid.size = dimlen;
nce(nc_inq_dimid(nc_file_id, "dst_grid_size", &ncdimid));
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
tgt_grid.size = dimlen;
status = nc_inq_dimid(nc_file_id, "src_grid_corners", &ncdimid);
if (status == NC_NOERR)
{
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
src_grid.num_cell_corners = dimlen;
src_grid.luse_cell_corners = true;
src_grid.lneed_cell_corners = true;
}
status = nc_inq_dimid(nc_file_id, "dst_grid_corners", &ncdimid);
if (status == NC_NOERR)
{
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
tgt_grid.num_cell_corners = dimlen;
tgt_grid.luse_cell_corners = true;
tgt_grid.lneed_cell_corners = true;
}
nce(nc_inq_dimid(nc_file_id, "src_grid_rank", &ncdimid));
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
src_grid.rank = dimlen;
nce(nc_inq_dimid(nc_file_id, "dst_grid_rank", &ncdimid));
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
tgt_grid.rank = dimlen;
nce(nc_inq_dimid(nc_file_id, "num_links", &ncdimid));
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
rv.num_links = dimlen;
/*
if ( rv.num_links == 0 )
cdoAbort("Number of remap links is 0, no remap weights found!");
*/
nce(nc_inq_dimid(nc_file_id, "num_wgts", &ncdimid));
nce(nc_inq_dimlen(nc_file_id, ncdimid, &dimlen));
rv.num_wts = dimlen;
remapGridWAlloc(rv.mapType, src_grid);
remapGridWAlloc(rv.mapType, tgt_grid);
// Allocate address and weight arrays for mapping 1
if (rv.num_links > 0)
{
rv.max_links = rv.num_links;
rv.src_cell_add.resize(rv.num_links);
rv.tgt_cell_add.resize(rv.num_links);
rv.wts.resize(rv.num_wts * rv.num_links);
}
// Get variable ids
nce(nc_inq_varid(nc_file_id, "dst_grid_dims", &nc_dstgrddims_id));
nce(nc_inq_varid(nc_file_id, "dst_grid_imask", &nc_dstgrdimask_id));
nce(nc_inq_varid(nc_file_id, "src_address", &nc_srcadd_id));
nce(nc_inq_varid(nc_file_id, "dst_address", &nc_dstadd_id));
nce(nc_inq_varid(nc_file_id, "remap_matrix", &nc_rmpmatrix_id));
// Read all variables
int ncvarid;
int dims[2];
nce(nc_inq_varid(nc_file_id, "src_grid_dims", &ncvarid));
nce(nc_get_var_int(nc_file_id, ncvarid, dims));
src_grid.dims[0] = dims[0];
src_grid.dims[1] = dims[1];
nce(nc_inq_varid(nc_file_id, "src_grid_imask", &ncvarid));
nce(nc_get_var_int(nc_file_id, ncvarid, &src_grid.mask[0]));
nc_read_coordinate_radian(nc_file_id, "src_grid_center_lat", src_grid.size, src_grid.cell_center_lat.data());
nc_read_coordinate_radian(nc_file_id, "src_grid_center_lon", src_grid.size, src_grid.cell_center_lon.data());
if (src_grid.num_cell_corners)
{
const size_t size = src_grid.num_cell_corners * src_grid.size;
nc_read_coordinate_radian(nc_file_id, "src_grid_corner_lat", size, src_grid.cell_corner_lat.data());
nc_read_coordinate_radian(nc_file_id, "src_grid_corner_lon", size, src_grid.cell_corner_lon.data());
}
if (lgridarea) nc_read_var_double(nc_file_id, "src_grid_area", src_grid.cell_area.data());
nc_read_var_double(nc_file_id, "src_grid_frac", src_grid.cell_frac.data());
nce(nc_get_var_int(nc_file_id, nc_dstgrddims_id, dims));
tgt_grid.dims[0] = dims[0];
tgt_grid.dims[1] = dims[1];
nce(nc_get_var_int(nc_file_id, nc_dstgrdimask_id, &tgt_grid.mask[0]));
nc_read_coordinate_radian(nc_file_id, "dst_grid_center_lat", tgt_grid.size, tgt_grid.cell_center_lat.data());
nc_read_coordinate_radian(nc_file_id, "dst_grid_center_lon", tgt_grid.size, tgt_grid.cell_center_lon.data());
if (tgt_grid.num_cell_corners)
{
const size_t size = tgt_grid.num_cell_corners * tgt_grid.size;
nc_read_coordinate_radian(nc_file_id, "dst_grid_corner_lat", size, tgt_grid.cell_corner_lat.data());
nc_read_coordinate_radian(nc_file_id, "dst_grid_corner_lon", size, tgt_grid.cell_corner_lon.data());
}
if (lgridarea) nc_read_var_double(nc_file_id, "dst_grid_area", tgt_grid.cell_area.data());
nc_read_var_double(nc_file_id, "dst_grid_frac", tgt_grid.cell_frac.data());
if (rv.num_links > 0)
{
readLinks(nc_file_id, nc_srcadd_id, rv.num_links, &rv.src_cell_add[0]);
readLinks(nc_file_id, nc_dstadd_id, rv.num_links, &rv.tgt_cell_add[0]);
for (size_t i = 0; i < rv.num_links; ++i)
{
rv.src_cell_add[i]--;
rv.tgt_cell_add[i]--;
}
nce(nc_get_var_double(nc_file_id, nc_rmpmatrix_id, &rv.wts[0]));
}
// Close input file
nce(nc_close(nc_file_id));
} // remapReadDataScrip
void
remapWeightsCheckAreas(size_t n_a, const std::vector<double> &area_a, const std::vector<double> &area_b,
size_t n_s, const size_t *col, const size_t *row, const double *S)
{
std::vector<double> sum(n_a, 0.0);
// sum weighted ratio of true areas
// ’a’ is source; ’b’ is destination
for (size_t i = 0; i < n_s; ++i) // loop over all elements of S (the weights)
sum[col[i]] = sum[col[i]] + S[i] * area_a[col[i]] / area_b[row[i]];
// check that sums are equal to 1 (within tolerance of 1.e-6)
for (size_t i = 0; i < n_a; ++i) // loop over all source cells
if (std::fabs(sum[i] - 1.0) > 1.e-6)
printf("ERROR\n");
printf("OK\n");
}
static
void verifyWeights(const char *weights_file)
{
// int nc_file_id = cdf_openread(weights_file);
// nce(nc_close(nc_file_id));
//remapReadDataScrip(weights_file, gridID1, gridID2, &mapType, &submapType, &numNeighbors, remapOrder, remaps[0].src_grid,
// remaps[0].tgt_grid, remaps[0].vars);
RemapMethod mapType(RemapMethod::UNDEF);
SubmapType submapType(SubmapType::NONE);
int numNeighbors = 0;
int remapOrder = 0;
RemapGridW src_grid, tgt_grid;
RemapVarsW remapvars;
remapReadDataScrip(weights_file, &mapType, &submapType, &numNeighbors, remapOrder, src_grid, tgt_grid, remapvars);
remapWeightsCheckAreas(src_grid.size, src_grid.cell_area, tgt_grid.cell_area,
remapvars.num_links, remapvars.src_cell_add.data(), remapvars.tgt_cell_add.data(), remapvars.wts.data());
}
#endif
......
......@@ -60,7 +60,7 @@ writeLinks(int nc_file_id, int nc_add_id, nc_type sizetype, size_t num_links, si
#endif
}
static void
void
readLinks(int nc_file_id, int nc_add_id, size_t num_links, size_t *cell_add)
{
if (num_links < 0x7FFFFC00) // 2GB
......@@ -437,7 +437,7 @@ remapWriteDataScrip(const char *interp_file, RemapMethod mapType, SubmapType sub
/*****************************************************************************/
#ifdef HAVE_LIBNETCDF
static RemapMethod
RemapMethod
getMapType(int nc_file_id, SubmapType *submapType, int *numNeighbors, int &remapOrder)
{
// Map method
......@@ -532,10 +532,10 @@ remapReadDataScrip(const char *interp_file, int gridID1, int gridID2, RemapMetho
int nc_rmpmatrix_id; /* id for remapping matrix */
char map_name[1024];
char normalize_opt[64]; /* character string for normalization option */
char convention[64]; /* character string for output convention */
char src_grid_name[64]; /* grid name for source grid */
char tgt_grid_name[64]; /* grid name for dest grid */
char normalize_opt[64]; /* character string for normalization option */
char convention[64]; /* character string for output convention */
char src_grid_name[64]; /* grid name for source grid */
char tgt_grid_name[64]; /* grid name for dest grid */
char src_grid_units[64];
char tgt_grid_units[64];
size_t attlen, dimlen;
......@@ -575,17 +575,16 @@ remapReadDataScrip(const char *interp_file, int gridID1, int gridID2, RemapMetho
nce(nc_inq_attlen(nc_file_id, NC_GLOBAL, "normalization", &attlen));
normalize_opt[attlen] = 0;
if (cstrIsEqual(normalize_opt, "none"))
rv.normOpt = NormOpt::NONE;
else if (cstrIsEqual(normalize_opt, "fracarea"))
rv.normOpt = NormOpt::FRACAREA;
else if (cstrIsEqual(normalize_opt, "destarea"))
rv.normOpt = NormOpt::DESTAREA;
// clang-format off
if (cstrIsEqual(normalize_opt, "none")) rv.normOpt = NormOpt::NONE;
else if (cstrIsEqual(normalize_opt, "fracarea")) rv.normOpt = NormOpt::FRACAREA;
else if (cstrIsEqual(normalize_opt, "destarea")) rv.normOpt = NormOpt::DESTAREA;
else
{
cdoPrint("normalize_opt = %s", normalize_opt);
cdoAbort("Invalid normalization option");
}
// clang-format on
if (Options::cdoVerbose) cdoPrint("normalize_opt = %s", normalize_opt);
......
......@@ -53,9 +53,9 @@ struct remaplink_t
struct RemapVars
{
bool sort_add;
bool pinit; // true: if the pointers are initialized
RemapMethod mapType; // identifier for remapping method
NormOpt normOpt; // option for normalization (conserv only)
bool pinit; // true: if the pointers are initialized
RemapMethod mapType; // identifier for remapping method
NormOpt normOpt; // option for normalization (conserv only)
long links_per_value;
size_t max_links; // current size of link arrays
size_t num_links; // actual number of links for remapping
......@@ -71,10 +71,8 @@ struct RemapVars
void remap(double * dst_array, double missval, size_t dst_size, const RemapVars &rv, const double * src_array,
RemapGradients &gradients);
void remap_laf(double * dst_array, double missval, size_t dst_size, const RemapVars &rv,
const double * src_array);
void remap_sum(double * dst_array, double missval, size_t dst_size, const RemapVars &rv,
const double * src_array);
void remap_laf(double * dst_array, double missval, size_t dst_size, const RemapVars &rv, const double * src_array);
void remap_sum(double * dst_array, double missval, size_t dst_size, const RemapVars &rv, const double * src_array);
void remapVarsInit(RemapMethod mapType, int remapOrder, RemapVars &rv);
void remapVarsEnsureSize(RemapVars &rv, size_t size);
void remapVarsResize(RemapVars &rv, size_t size);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment