Commit 8bbce5d9 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Replaced grid_to_radian() by cdo_grid_to_radian().

parent 72d0a473
Pipeline #2481 passed with stages
in 15 minutes and 5 seconds
......@@ -36,7 +36,6 @@ Arithlat(void *process)
int nrecs;
int varID, levelID;
size_t nmiss;
char units[CDI_MAX_NAME];
Varray<double> scale;
cdoInitialize(process);
......@@ -109,8 +108,7 @@ Arithlat(void *process)
gridInqYvals(gridID, &scale[0]);
// Convert lat/lon units if required
gridInqXunits(gridID, units);
grid_to_radian(units, gridsize, &scale[0], "grid latitudes");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, &scale[0], "grid latitudes");
if (operfunc == func_mul)
for (size_t i = 0; i < gridsize; ++i) scale[i] = std::cos(scale[i]);
......
......@@ -115,12 +115,9 @@ CDIwrite(void *process)
gridInqXvals(gridID2, &xvals[0]);
gridInqYvals(gridID2, &yvals[0]);
/* Convert lat/lon units if required */
char units[CDI_MAX_NAME];
gridInqXunits(gridID2, units);
grid_to_radian(units, gridsize, &xvals[0], "grid center lon");
gridInqYunits(gridID2, units);
grid_to_radian(units, gridsize, &yvals[0], "grid center lat");
// Convert lat/lon units if required
cdo_grid_to_radian(gridID2, CDI_XAXIS, gridsize, &xvals[0], "grid center lon");
cdo_grid_to_radian(gridID2, CDI_YAXIS, gridsize, &yvals[0], "grid center lat");
for (size_t i = 0; i < gridsize; i++) array[i] = 2 - std::cos(std::acos(std::cos(xvals[i]) * std::cos(yvals[i])) / 1.2);
......
......@@ -375,11 +375,8 @@ setmisstodis(Field &field1, Field &field2, int numNeighbors)
gridInqYvals(gridID, yvals.data());
// Convert lat/lon units if required
char units[CDI_MAX_NAME];
gridInqXunits(gridID, units);
grid_to_radian(units, gridsize, &xvals[0], "grid center lon");
gridInqYunits(gridID, units);
grid_to_radian(units, gridsize, &yvals[0], "grid center lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, &xvals[0], "grid center lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize, &yvals[0], "grid center lat");
std::vector<size_t> mindex(nmiss, 1), vindex(nvals, 1);
Varray<double> lons(nvals), lats(nvals);
......
......@@ -212,9 +212,6 @@ Gridcell(void *process)
if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_LCC || projtype == CDI_PROJ_STERE
|| gridtype == GRID_CURVILINEAR)
{
double len1 = 0, len2 = 0;
char units[CDI_MAX_NAME];
if (gridtype != GRID_CURVILINEAR) gridID = gridToCurvilinear(gridID, 1);
gridsize = gridInqSize(gridID);
......@@ -229,13 +226,13 @@ Gridcell(void *process)
gridInqYvals(gridID, yv.data());
// Convert lat/lon units if required
gridInqXunits(gridID, units);
grid_to_radian(units, gridsize, &xv[0], "grid longitudes");
grid_to_radian(units, gridsize, &yv[0], "grid latitudes");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, &xv[0], "grid longitudes");
cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize, &yv[0], "grid latitudes");
const auto planetRadius = getPlanetRadius(gridID);
if (Options::cdoVerbose) cdoPrint("Planet radius: %.2f", planetRadius);
double len1 = 0, len2 = 0;
if (operatorID == GRIDDX)
{
for (size_t j = 0; j < ysize; ++j)
......
......@@ -73,14 +73,8 @@ grid_new(int gridID, const char *txt)
gridInqXbounds(grid->gridID, grid->cell_corner_lon);
gridInqYbounds(grid->gridID, grid->cell_corner_lat);
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
int length = CDI_MAX_NAME;
cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, xunits, &length);
length = CDI_MAX_NAME;
cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, yunits, &length);
grid_to_radian(xunits, grid->num_cell_corners * grid->size, grid->cell_corner_lon, "grid corner lon");
grid_to_radian(yunits, grid->num_cell_corners * grid->size, grid->cell_corner_lat, "grid corner lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, grid->num_cell_corners * grid->size, grid->cell_corner_lon, "grid corner lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, grid->num_cell_corners * grid->size, grid->cell_corner_lat, "grid corner lat");
if (lgrid_destroy) gridDestroy(gridID);
......
......@@ -281,22 +281,17 @@ read_coordinates(const char *filename, long n, double *lon, double *lat, int nv,
gridInqXvals(gridID, lon);
gridInqYvals(gridID, lat);
char units[CDI_MAX_NAME];
/* Convert lat/lon units if required */
gridInqXunits(gridID, units);
grid_to_radian(units, n, lon, "grid center lon");
gridInqYunits(gridID, units);
grid_to_radian(units, n, lat, "grid center lat");
// Convert lat/lon units if required
cdo_grid_to_radian(gridID, CDI_XAXIS, n, lon, "grid center lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, n, lat, "grid center lat");
if (nv == 3 && lon_bnds && lat_bnds)
{
gridInqXbounds(gridID, lon_bnds);
gridInqYbounds(gridID, lat_bnds);
gridInqXunits(gridID, units);
grid_to_radian(units, n * 3, lon_bnds, "grid corner lon");
gridInqYunits(gridID, units);
grid_to_radian(units, n * 3, lat_bnds, "grid corner lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, n * 3, lon_bnds, "grid corner lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, n * 3, lat_bnds, "grid corner lat");
}
streamClose(streamID);
......
......@@ -89,11 +89,8 @@ smooth(int gridID, double missval, const double *restrict array1, double *restri
gridInqYvals(gridID, &yvals[0]);
// Convert lat/lon units if required
char units[CDI_MAX_NAME];
gridInqXunits(gridID, units);
grid_to_radian(units, gridsize, &xvals[0], "grid center lon");
gridInqYunits(gridID, units);
grid_to_radian(units, gridsize, &yvals[0], "grid center lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, &xvals[0], "grid center lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize, &yvals[0], "grid center lat");
std::vector<knnWeightsType> knnWeights;
for (int i = 0; i < Threading::ompNumThreads; ++i) knnWeights.push_back(knnWeightsType(numNeighbors));
......
......@@ -479,11 +479,8 @@ Vargen(void *process)
gridInqYvals(gridID, yvals.data());
// Convert lat/lon units if required
char units[CDI_MAX_NAME];
gridInqXunits(gridID, units);
grid_to_radian(units, gridsize, xvals.data(), "grid center lon");
gridInqYunits(gridID, units);
grid_to_radian(units, gridsize, yvals.data(), "grid center lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, xvals.data(), "grid center lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize, yvals.data(), "grid center lat");
}
if (operatorID == SINCOS)
......
......@@ -306,8 +306,7 @@ gridGenAreaReg2D(const int gridID, double *area, bool lweights)
auto nlat = gridInqYsize(gridID);
if (!nlon) nlon = 1;
if (!nlat) nlat = 1;
std::vector<double> grid_corner_lon(nlon * 2);
std::vector<double> grid_corner_lat(nlat * 2);
std::vector<double> grid_corner_lon(nlon * 2), grid_corner_lat(nlat * 2);
if (gridInqXbounds(gridID, nullptr))
{
......@@ -450,19 +449,11 @@ gridGenAreaUnstruct(int gridID, double *area)
return 1;
}
char xunitstr[CDI_MAX_NAME];
char yunitstr[CDI_MAX_NAME];
gridInqXunits(gridID, xunitstr);
gridInqYunits(gridID, yunitstr);
std::vector<double> grid_center_lon(gridsize);
std::vector<double> grid_center_lat(gridsize);
std::vector<double> grid_center_lon(gridsize), grid_center_lat(gridsize);
gridInqXvals(gridID, grid_center_lon.data());
gridInqYvals(gridID, grid_center_lat.data());
std::vector<double> grid_corner_lon(nv * gridsize);
std::vector<double> grid_corner_lat(nv * gridsize);
std::vector<double> grid_corner_lon(nv * gridsize), grid_corner_lat(nv * gridsize);
if (gridInqYbounds(gridID, nullptr) && gridInqXbounds(gridID, nullptr))
{
......@@ -475,11 +466,11 @@ gridGenAreaUnstruct(int gridID, double *area)
return status;
}
grid_to_radian(xunitstr, gridsize, grid_center_lon.data(), "grid1 center longitudes");
grid_to_radian(xunitstr, gridsize * nv, grid_corner_lon.data(), "grid1 corner longitudes");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, grid_center_lon.data(), "grid1 center longitudes");
cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize * nv, grid_corner_lon.data(), "grid1 corner longitudes");
grid_to_radian(yunitstr, gridsize, grid_center_lat.data(), "grid1 center latitudes");
grid_to_radian(yunitstr, gridsize * nv, grid_corner_lat.data(), "grid1 corner latitudes");
cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize, grid_center_lat.data(), "grid1 center latitudes");
cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize * nv, grid_corner_lat.data(), "grid1 corner latitudes");
if (lgriddestroy) gridDestroy(gridID);
......
......@@ -303,7 +303,6 @@ intlinarr(long nxm, double *ym, double *xm, int nx, double *y, double *x)
void
intgridbil(Field &field1, Field &field2)
{
char xunits[CDI_MAX_NAME];
auto gridID1 = field1.grid;
auto gridID2 = field2.grid;
auto &array1 = field1.vec_d;
......@@ -338,10 +337,8 @@ intgridbil(Field &field1, Field &field2)
{
if (lon_is_circular) lon1[nlon1 - 1] = 0;
gridInqXunits(gridID1, xunits);
grid_to_radian(xunits, nlon1, lon1.data(), "grid1 center lon");
grid_to_radian(xunits, nlat1, lat1.data(), "grid1 center lat");
cdo_grid_to_radian(gridID1, CDI_XAXIS, nlon1, lon1.data(), "grid1 center lon");
cdo_grid_to_radian(gridID1, CDI_YAXIS, nlat1, lat1.data(), "grid1 center lat");
if (lon_is_circular) lon1[nlon1 - 1] = lon1[0] + 2 * M_PI;
}
......@@ -355,10 +352,8 @@ intgridbil(Field &field1, Field &field2)
gridInqXvals(gridID2, &lon2);
gridInqYvals(gridID2, &lat2);
gridInqXunits(gridID2, xunits);
grid_to_radian(xunits, xsize2, &lon2, "grid2 center lon");
grid_to_radian(xunits, ysize2, &lat2, "grid2 center lat");
cdo_grid_to_radian(gridID2, CDI_XAXIS, xsize2, &lon2, "grid2 center lon");
cdo_grid_to_radian(gridID2, CDI_YAXIS, ysize2, &lat2, "grid2 center lat");
if (lon2 < lon1[0]) lon2 += 2 * M_PI;
......@@ -415,10 +410,8 @@ intgridbil(Field &field1, Field &field2)
gridInqXvals(gridID2, xvals2.data());
gridInqYvals(gridID2, yvals2.data());
gridInqXunits(gridID2, xunits);
grid_to_radian(xunits, gridsize2, xvals2.data(), "grid2 center lon");
grid_to_radian(xunits, gridsize2, yvals2.data(), "grid2 center lat");
cdo_grid_to_radian(gridID2, CDI_XAXIS, gridsize2, xvals2.data(), "grid2 center lon");
cdo_grid_to_radian(gridID2, CDI_YAXIS, gridsize2, yvals2.data(), "grid2 center lat");
for (size_t i = 0; i < gridsize2; ++i)
{
......@@ -428,8 +421,7 @@ intgridbil(Field &field1, Field &field2)
}
else
{
Varray<double> xcoord(xsize2);
Varray<double> ycoord(ysize2);
Varray<double> xcoord(xsize2), ycoord(ysize2);
gridInqXvals(gridID2, xcoord.data());
gridInqYvals(gridID2, ycoord.data());
......
......@@ -197,6 +197,15 @@ grid_to_degree(const char *units, size_t nvals, double *restrict values, const c
}
}
void
cdo_grid_to_radian(int gridID, int varID, size_t nvals, double *restrict values, const char *description)
{
char units[CDI_MAX_NAME];
int length = CDI_MAX_NAME;
cdiInqKeyString(gridID, varID, CDI_KEY_UNITS, units, &length);
grid_to_radian(units, nvals, values, description);
}
void
cdo_grid_to_degree(int gridID, int varID, size_t nvals, double *restrict values, const char *description)
{
......
......@@ -49,6 +49,7 @@ void grid_copy_mapping(int gridID1, int gridID2);
bool grid_is_distance_generic(int gridID);
void grid_to_radian(const char *units, size_t nvals, double *values, const char *description);
void cdo_grid_to_radian(int gridID, int varID, size_t nvals, double *values, const char *description);
void cdo_grid_to_degree(int gridID, int varID, size_t nvals, double *values, const char *description);
void grid_gen_corners(size_t n, const double *vals, double *corners);
......
......@@ -302,14 +302,8 @@ remap_define_reg2d(int gridID, RemapGrid &grid)
gridInqYvals(gridID, grid.reg2d_center_lat);
// Convert lat/lon units if required
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
int length = CDI_MAX_NAME;
cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, xunits, &length);
length = CDI_MAX_NAME;
cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, yunits, &length);
grid_to_radian(xunits, nx, grid.reg2d_center_lon, "grid reg2d center lon");
grid_to_radian(yunits, ny, grid.reg2d_center_lat, "grid reg2d center lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, nx, grid.reg2d_center_lon, "grid reg2d center lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, ny, grid.reg2d_center_lat, "grid reg2d center lat");
if (grid.reg2d_center_lon[nx - 1] < grid.reg2d_center_lon[0])
for (size_t i = 1; i < nx; ++i)
......@@ -326,7 +320,7 @@ remap_define_reg2d(int gridID, RemapGrid &grid)
gridInqXbounds(gridID, xbounds.data());
grid.reg2d_corner_lon[0] = xbounds[1];
for (size_t i = 0; i < nx; ++i) grid.reg2d_corner_lon[i + 1] = xbounds[2 * i];
grid_to_radian(xunits, nx + 1, grid.reg2d_corner_lon, "grid reg2d corner lon");
cdo_grid_to_radian(gridID, CDI_XAXIS, nx + 1, grid.reg2d_corner_lon, "grid reg2d corner lon");
}
else
{
......@@ -339,7 +333,7 @@ remap_define_reg2d(int gridID, RemapGrid &grid)
gridInqYbounds(gridID, ybounds.data());
grid.reg2d_corner_lat[0] = ybounds[1];
for (size_t j = 0; j < ny; ++j) grid.reg2d_corner_lat[j + 1] = ybounds[2 * j];
grid_to_radian(yunits, ny + 1, grid.reg2d_corner_lat, "grid reg2d corner lat");
cdo_grid_to_radian(gridID, CDI_YAXIS, ny + 1, grid.reg2d_corner_lat, "grid reg2d corner lat");
}
else
{
......@@ -417,12 +411,6 @@ remapDefineGrid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *tx
gridInqXvals(gridID, grid.cell_center_lon);
gridInqYvals(gridID, grid.cell_center_lat);
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
int length = CDI_MAX_NAME;
cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, xunits, &length);
length = CDI_MAX_NAME;
cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, yunits, &length);
if (grid.lneed_cell_corners)
{
if (gridInqXbounds(gridID, nullptr) && gridInqYbounds(gridID, nullptr))
......@@ -439,14 +427,12 @@ remapDefineGrid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *tx
if (gridInqType(grid.gridID) == GRID_GME) gridInqMaskGME(gridID_gme, &grid.vgpm[0]);
// Convert lat/lon units if required
grid_to_radian(xunits, grid.size, grid.cell_center_lon, "grid center lon");
grid_to_radian(yunits, grid.size, grid.cell_center_lat, "grid center lat");
// Note: using units from cell center instead from bounds
cdo_grid_to_radian(gridID, CDI_XAXIS, grid.size, grid.cell_center_lon, "grid center lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, grid.size, grid.cell_center_lat, "grid center lat");
if (grid.num_cell_corners && grid.lneed_cell_corners)
{
grid_to_radian(xunits, grid.num_cell_corners * grid.size, grid.cell_corner_lon, "grid corner lon");
grid_to_radian(yunits, grid.num_cell_corners * grid.size, grid.cell_corner_lat, "grid corner lat");
cdo_grid_to_radian(gridID, CDI_XAXIS, grid.num_cell_corners * grid.size, grid.cell_corner_lon, "grid corner lon");
cdo_grid_to_radian(gridID, CDI_YAXIS, grid.num_cell_corners * grid.size, grid.cell_corner_lat, "grid corner lat");
}
if (lgrid_destroy) gridDestroy(gridID);
......
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