Commit 0ec5cdbb authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapReadDataScrip: refactor.

parent 76053e6e
......@@ -36,27 +36,69 @@ nce(int istat)
if (istat != NC_NOERR) cdoAbort(nc_strerror(istat));
}
size_t
cdfReadDimlen(int ncfileid, const char *dimname)
{
size_t dimlen = 0;
int ncdimid;
const int status = nc_inq_dimid(ncfileid, dimname, &ncdimid);
if (status == NC_NOERR) nce(nc_inq_dimlen(ncfileid, ncdimid, &dimlen));
return dimlen;
}
void
cdfReadAttText(int ncfileid, int ncvarid, const char *name, char *text)
{
size_t attlen;
nce(nc_get_att_text(ncfileid, ncvarid, name, text));
nce(nc_inq_attlen(ncfileid, ncvarid, name, &attlen));
text[attlen] = 0;
}
void
nc_read_var_double(int nc_file_id, const char *name, double *array)
cdfReadVarInt(int ncfileid, const char *name, int *array)
{
int ncvarid;
nce(nc_inq_varid(nc_file_id, name, &ncvarid));
nce(nc_get_var_double(nc_file_id, ncvarid, array));
nce(nc_inq_varid(ncfileid, name, &ncvarid));
nce(nc_get_var_int(ncfileid, ncvarid, array));
}
void
nc_read_coordinate_radian(int nc_file_id, const char *name, size_t size, double *array)
cdfReadVarSize(int ncfileid, const char *name, size_t len, size_t *array)
{
if (len < 0x7FFFFC00) // 2GB
{
std::vector<int> iarray(len);
cdfReadVarInt(ncfileid, name, iarray.data());
for (size_t i = 0; i < len; ++i) array[i] = (size_t) iarray[i];
}
#ifdef HAVE_NETCDF4
else
{
int ncvarid;
nce(nc_inq_varid(ncfileid, name, &ncvarid));
nce(nc_get_var_ulonglong(ncfileid, ncvarid, (unsigned long long *) array));
}
#endif
}
void
cdfReadVarDouble(int ncfileid, 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));
nce(nc_inq_varid(ncfileid, name, &ncvarid));
nce(nc_get_var_double(ncfileid, 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;
void
cdfReadCoordinateRadian(int ncfileid, const char *name, size_t size, double *array)
{
int ncvarid;
nce(nc_inq_varid(ncfileid, name, &ncvarid));
nce(nc_get_var_double(ncfileid, ncvarid, array));
char grid_units[64];
cdfReadAttText(ncfileid, ncvarid, "units", grid_units);
grid_to_radian(grid_units, size, array, name);
}
......@@ -145,6 +187,25 @@ struct RemapVarsW
std::vector<double> wts; // map weights for each link [max_links*num_wts]
};
NormOpt
getNormOpt(const char *normOptStr)
{
NormOpt normOpt(NormOpt::NONE);
// clang-format off
if (cstrIsEqual(normOptStr, "none")) normOpt = NormOpt::NONE;
else if (cstrIsEqual(normOptStr, "fracarea")) normOpt = NormOpt::FRACAREA;
else if (cstrIsEqual(normOptStr, "destarea")) normOpt = NormOpt::DESTAREA;
else
{
cdoPrint("normalize_opt = %s", normOptStr);
cdoAbort("Invalid normalization option");
}
// clang-format on
if (Options::cdoVerbose) cdoPrint("normalize_opt = %s", normOptStr);
return normOpt;
}
void
remapVarsWInit(RemapMethod mapType, int remapOrder, RemapVarsW &rv)
......@@ -164,8 +225,7 @@ remapVarsWInit(RemapMethod mapType, int remapOrder, RemapVarsW &rv)
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);
RemapMethod getMapType(int ncfileid, SubmapType *submapType, int *numNeighbors, int &remapOrder);
static void
remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *submapType,
......@@ -174,24 +234,15 @@ remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *su
// 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;
size_t 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);
// nce(nc_open(interp_file, NC_NOWRITE, &ncfileid));
int ncfileid = 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;
cdfReadAttText(ncfileid, NC_GLOBAL, "title", map_name);
if (Options::cdoVerbose)
{
......@@ -200,7 +251,7 @@ remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *su
}
// Map Type
*mapType = getMapType(nc_file_id, submapType, numNeighbors, remapOrder);
*mapType = getMapType(ncfileid, submapType, numNeighbors, remapOrder);
if (*mapType == RemapMethod::CONSERV_SCRIP) lgridarea = true;
remapVarsWInit(*mapType, remapOrder, rv);
......@@ -211,28 +262,12 @@ remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *su
// 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);
cdfReadAttText(ncfileid, NC_GLOBAL, "normalization", normalize_opt);
rv.normOpt = getNormOpt(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;
cdfReadAttText(ncfileid, NC_GLOBAL, "conventions", convention);
if (!cstrIsEqual(convention, "SCRIP"))
{
......@@ -247,14 +282,10 @@ remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *su
// 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;
cdfReadAttText(ncfileid, NC_GLOBAL, "source_grid", src_grid_name);
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;
cdfReadAttText(ncfileid, NC_GLOBAL, "dest_grid", tgt_grid_name);
if (Options::cdoVerbose) cdoPrint("Remapping between: %s and %s", src_grid_name, tgt_grid_name);
......@@ -263,51 +294,30 @@ remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *su
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;
src_grid.size = cdfReadDimlen(ncfileid, "src_grid_size");
tgt_grid.size = cdfReadDimlen(ncfileid, "dst_grid_size");
status = nc_inq_dimid(nc_file_id, "src_grid_corners", &ncdimid);
if (status == NC_NOERR)
dimlen = cdfReadDimlen(ncfileid, "src_grid_corners");
if (dimlen > 0)
{
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)
dimlen = cdfReadDimlen(ncfileid, "dst_grid_corners");
if (dimlen > 0)
{
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;
src_grid.rank = cdfReadDimlen(ncfileid, "src_grid_rank");;
tgt_grid.rank = cdfReadDimlen(ncfileid, "dst_grid_rank");;
rv.num_links = cdfReadDimlen(ncfileid, "num_links");
// if ( rv.num_links == 0 ) cdoAbort("Number of remap links is 0, no remap weights found!");
rv.num_wts = cdfReadDimlen(ncfileid, "num_wgts");
remapGridWAlloc(rv.mapType, src_grid);
remapGridWAlloc(rv.mapType, tgt_grid);
......@@ -321,74 +331,55 @@ remapReadDataScrip(const char *interp_file, RemapMethod *mapType, SubmapType *su
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));
// Read all variables of source grid
cdfReadVarSize(ncfileid, "src_grid_dims", 2, src_grid.dims);
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));
cdfReadVarInt(ncfileid, "src_grid_imask", src_grid.mask.data());
// 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());
cdfReadCoordinateRadian(ncfileid, "src_grid_center_lat", src_grid.size, src_grid.cell_center_lat.data());
cdfReadCoordinateRadian(ncfileid, "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());
cdfReadCoordinateRadian(ncfileid, "src_grid_corner_lat", size, src_grid.cell_corner_lat.data());
cdfReadCoordinateRadian(ncfileid, "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());
if (lgridarea) cdfReadVarDouble(ncfileid, "src_grid_area", src_grid.cell_area.data());
cdfReadVarDouble(ncfileid, "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];
// Read all variables of target grid
cdfReadVarSize(ncfileid, "dst_grid_dims", 2, tgt_grid.dims);
nce(nc_get_var_int(nc_file_id, nc_dstgrdimask_id, &tgt_grid.mask[0]));
cdfReadVarInt(ncfileid, "dst_grid_imask", tgt_grid.mask.data());
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());
cdfReadCoordinateRadian(ncfileid, "dst_grid_center_lat", tgt_grid.size, tgt_grid.cell_center_lat.data());
cdfReadCoordinateRadian(ncfileid, "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());
cdfReadCoordinateRadian(ncfileid, "dst_grid_corner_lat", size, tgt_grid.cell_corner_lat.data());
cdfReadCoordinateRadian(ncfileid, "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 (lgridarea) cdfReadVarDouble(ncfileid, "dst_grid_area", tgt_grid.cell_area.data());
cdfReadVarDouble(ncfileid, "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]);
cdfReadVarSize(ncfileid, "src_address", rv.num_links, rv.src_cell_add.data());
cdfReadVarSize(ncfileid, "dst_address", rv.num_links, rv.tgt_cell_add.data());
for (size_t i = 0; i < rv.num_links; ++i)
{
rv.src_cell_add[i]--;
rv.tgt_cell_add[i]--;
}
for (size_t i = 0; i < rv.num_links; ++i) rv.src_cell_add[i]--;
for (size_t i = 0; i < rv.num_links; ++i) rv.tgt_cell_add[i]--;
nce(nc_get_var_double(nc_file_id, nc_rmpmatrix_id, &rv.wts[0]));
cdfReadVarDouble(ncfileid, "remap_matrix", rv.wts.data());
}
// Close input file
nce(nc_close(nc_file_id));
nce(nc_close(ncfileid));
} // remapReadDataScrip
void
......@@ -412,8 +403,8 @@ remapWeightsCheckAreas(size_t n_a, const std::vector<double> &area_a, const std:
static
void verifyWeights(const char *weights_file)
{
// int nc_file_id = cdf_openread(weights_file);
// nce(nc_close(nc_file_id));
// int ncfileid = cdf_openread(weights_file);
// nce(nc_close(ncfileid));
RemapMethod mapType(RemapMethod::UNDEF);
SubmapType submapType(SubmapType::NONE);
......
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