Commit 0804ebf5 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

grid_gen_bounds: changed type of parameter to vector.

parent e0b72162
......@@ -328,7 +328,7 @@ gridGenAreaReg2D(const int gridID, double *area, bool lweights)
{
std::vector<double> grid_center_lon(nlon);
gridInqXvals(gridID, grid_center_lon.data());
grid_gen_bounds(nlon, grid_center_lon.data(), grid_corner_lon.data());
grid_gen_bounds(nlon, grid_center_lon, grid_corner_lon);
gridInqXunits(gridID, unitstr);
}
......@@ -351,7 +351,7 @@ gridGenAreaReg2D(const int gridID, double *area, bool lweights)
{
std::vector<double> grid_center_lat(nlat);
gridInqYvals(gridID, grid_center_lat.data());
grid_gen_bounds(nlat, grid_center_lat.data(), grid_corner_lat.data());
grid_gen_bounds(nlat, grid_center_lat, grid_corner_lat);
grid_check_lat_borders(nlat * 2, grid_corner_lat.data());
gridInqYunits(gridID, unitstr);
}
......
......@@ -200,17 +200,16 @@ gen_grid_lonlat(GridDesciption &grid, const char *pline, double inc, double lon1
if (lbounds && nlon > 1 && nlat > 1)
{
std::vector<double> xbounds(2 * nlon);
std::vector<double> ybounds(2 * nlat);
std::vector<double> xbounds(2 * nlon), ybounds(2 * nlat);
grid_gen_bounds(nlon, grid.xvals.data(), xbounds.data());
grid_gen_bounds(nlat, yvals.data(), ybounds.data());
grid_gen_bounds(nlon, grid.xvals, xbounds);
grid_gen_bounds(nlat, yvals, ybounds);
grid_check_lat_borders(2 * nlat, ybounds.data());
grid.xbounds.resize(4 * gridsize);
grid.ybounds.resize(4 * gridsize);
grid_gen_xbounds2D(nlon, nlat, xbounds.data(), grid.xbounds.data());
grid_gen_ybounds2D(nlon, nlat, ybounds.data(), grid.ybounds.data());
grid_gen_xbounds2D(nlon, nlat, xbounds, grid.xbounds);
grid_gen_ybounds2D(nlon, nlat, ybounds, grid.ybounds);
}
}
}
......
......@@ -319,7 +319,7 @@ grid_gen_corners(size_t n, const double *restrict vals, double *restrict corners
}
void
grid_gen_bounds(size_t n, const double *restrict vals, double *restrict bounds)
grid_gen_bounds(size_t n, const std::vector<double> &vals, std::vector<double> &bounds)
{
for (size_t i = 0; i < n - 1; ++i)
{
......@@ -462,8 +462,8 @@ grid_cell_center_to_bounds_Y2D(const char *yunitstr, size_t xsize, size_t ysize,
}
static void
gridGenRotBounds(double xpole, double ypole, double angle, size_t nx, size_t ny, double *xbounds, double *ybounds,
double *xbounds2D, double *ybounds2D)
gridGenRotBounds(double xpole, double ypole, double angle, size_t nx, size_t ny, const std::vector<double> &xbounds,
const std::vector<double> &ybounds, std::vector<double> &xbounds2D, std::vector<double> &ybounds2D)
{
double minlon, maxlon;
double minlat, maxlat;
......@@ -501,7 +501,7 @@ gridGenRotBounds(double xpole, double ypole, double angle, size_t nx, size_t ny,
}
void
grid_gen_xbounds2D(size_t nx, size_t ny, const double *restrict xbounds, double *restrict xbounds2D)
grid_gen_xbounds2D(size_t nx, size_t ny, const std::vector<double> &xbounds, std::vector<double> &xbounds2D)
{
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(nx, ny, xbounds, xbounds2D)
......@@ -513,7 +513,7 @@ grid_gen_xbounds2D(size_t nx, size_t ny, const double *restrict xbounds, double
for (size_t j = 0; j < ny; ++j)
{
const size_t index = j * 4 * nx + 4 * i;
const size_t index = 4 * (j * nx + i);
xbounds2D[index] = minlon;
xbounds2D[index + 1] = maxlon;
xbounds2D[index + 2] = maxlon;
......@@ -523,7 +523,7 @@ grid_gen_xbounds2D(size_t nx, size_t ny, const double *restrict xbounds, double
}
void
grid_gen_ybounds2D(size_t nx, size_t ny, const double *restrict ybounds, double *restrict ybounds2D)
grid_gen_ybounds2D(size_t nx, size_t ny, const std::vector<double> &ybounds, std::vector<double> &ybounds2D)
{
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(nx, ny, ybounds, ybounds2D)
......@@ -544,7 +544,7 @@ grid_gen_ybounds2D(size_t nx, size_t ny, const double *restrict ybounds, double
for (size_t i = 0; i < nx; ++i)
{
const size_t index = j * 4 * nx + 4 * i;
const size_t index = 4 * (j * nx + i);
ybounds2D[index] = minlat;
ybounds2D[index + 1] = minlat;
ybounds2D[index + 2] = maxlat;
......@@ -957,20 +957,52 @@ check_units(const char *name, const char *units)
len ? "Invalid" : "Missing", name, len ? "\"" : "", units, len ? "\"" : "");
}
static void
center_1D_to_2D(size_t nx, size_t ny, const std::vector<double> &xvals, const std::vector<double> &yvals,
std::vector<double> &xvals2D, std::vector<double> &yvals2D, double xscale, double yscale)
{
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
xvals2D[j * nx + i] = xscale * xvals[i];
yvals2D[j * nx + i] = yscale * yvals[j];
}
}
static void
bounds_1D_to_2D(size_t nx, size_t ny, const std::vector<double> &xbounds, const std::vector<double> &ybounds,
std::vector<double> &xbounds2D, std::vector<double> &ybounds2D, double xscale, double yscale)
{
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
const size_t index = 4 * (j * nx + i);
xbounds2D[index + 0] = xscale * xbounds[2 * i];
ybounds2D[index + 0] = yscale * ybounds[2 * j];
xbounds2D[index + 1] = xscale * xbounds[2 * i];
ybounds2D[index + 1] = yscale * ybounds[2 * j + 1];
xbounds2D[index + 2] = xscale * xbounds[2 * i + 1];
ybounds2D[index + 2] = yscale * ybounds[2 * j + 1];
xbounds2D[index + 3] = xscale * xbounds[2 * i + 1];
ybounds2D[index + 3] = yscale * ybounds[2 * j];
}
}
int
gridToCurvilinear(int gridID1, int lbounds)
{
size_t index;
int gridtype = gridInqType(gridID1);
if (!(gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || gridtype == GRID_PROJECTION))
cdoAbort("Grid type >%s< unsupported!", gridNamePtr(gridtype));
size_t nx = gridInqXsize(gridID1);
size_t ny = gridInqYsize(gridID1);
bool lxyvals = gridInqXvals(gridID1, nullptr) && gridInqYvals(gridID1, nullptr);
const bool lxyvals = gridInqXvals(gridID1, nullptr) && gridInqYvals(gridID1, nullptr);
if (!lxyvals) cdoAbort("Grid coordinates missing!");
size_t gridsize = gridInqSize(gridID1);
int gridID2 = gridCreate(GRID_CURVILINEAR, gridsize);
const size_t gridsize = gridInqSize(gridID1);
const int gridID2 = gridCreate(GRID_CURVILINEAR, gridsize);
gridDefDatatype(gridID2, CDI_DATATYPE_FLT32);
char *proj4param = nullptr;
......@@ -990,7 +1022,7 @@ gridToCurvilinear(int gridID1, int lbounds)
}
else
{
int projtype = gridInqProjType(gridID1);
const int projtype = gridInqProjType(gridID1);
// clang-format off
if (projtype == CDI_PROJ_RLL) lproj_rll = true;
else if (projtype == CDI_PROJ_LAEA) lproj_laea = true;
......@@ -1009,182 +1041,151 @@ gridToCurvilinear(int gridID1, int lbounds)
}
}
if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN)
double xpole = 0, ypole = 0, angle = 0;
double xscale = 1, yscale = 1;
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
const size_t nvertex = (size_t) gridInqNvertex(gridID1);
gridInqXunits(gridID1, xunits);
gridInqYunits(gridID1, yunits);
if (lproj_laea || lproj_lcc || lproj_sinu || lproj_stere)
{
double xpole = 0, ypole = 0, angle = 0;
double xscale = 1, yscale = 1;
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
check_units("x", xunits);
check_units("y", yunits);
}
size_t nvertex = (size_t) gridInqNvertex(gridID1);
gridInqXunits(gridID1, xunits);
gridInqYunits(gridID1, yunits);
if (memcmp(xunits, "km", 2) == 0) xscale = 1000;
if (memcmp(yunits, "km", 2) == 0) yscale = 1000;
if (lproj_laea || lproj_lcc || lproj_sinu || lproj_stere)
{
check_units("x", xunits);
check_units("y", yunits);
}
gridDefXsize(gridID2, nx);
gridDefYsize(gridID2, ny);
if (memcmp(xunits, "km", 2) == 0) xscale = 1000;
if (memcmp(yunits, "km", 2) == 0) yscale = 1000;
std::vector<double> xvals2D(gridsize);
std::vector<double> yvals2D(gridsize);
gridDefXsize(gridID2, nx);
gridDefYsize(gridID2, ny);
if (nx == 0) nx = 1;
if (ny == 0) ny = 1;
std::vector<double> xvals2D(gridsize);
std::vector<double> yvals2D(gridsize);
std::vector<double> xvals(nx, 0);
std::vector<double> yvals(ny, 0);
if (nx == 0) nx = 1;
if (ny == 0) ny = 1;
if (gridInqXvals(gridID1, nullptr)) gridInqXvals(gridID1, xvals.data());
if (gridInqYvals(gridID1, nullptr)) gridInqYvals(gridID1, yvals.data());
std::vector<double> xvals(nx);
std::vector<double> yvals(ny);
if (lproj_rll)
{
gridInqParamRLL(gridID1, &xpole, &ypole, &angle);
gridDefProj(gridID2, gridID1);
if (gridInqXvals(gridID1, nullptr))
gridInqXvals(gridID1, xvals.data());
else
for (size_t i = 0; i < nx; ++i) xvals[i] = 0;
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
xvals2D[j * nx + i] = lamrot_to_lam(yvals[j], xvals[i], ypole, xpole, angle);
yvals2D[j * nx + i] = phirot_to_phi(yvals[j], xvals[i], ypole, angle);
}
}
else
{
center_1D_to_2D(nx, ny, xvals, yvals, xvals2D, yvals2D, xscale, yscale);
if (lproj_sinu)
cdo_sinu_to_lonlat(gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_laea)
cdo_laea_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_lcc)
cdo_lcc_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_stere)
cdo_stere_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj4)
cdo_proj_to_lonlat(proj4param, gridsize, xvals2D.data(), yvals2D.data());
}
if (gridInqYvals(gridID1, nullptr))
gridInqYvals(gridID1, yvals.data());
else
for (size_t i = 0; i < ny; ++i) yvals[i] = 0;
gridDefXvals(gridID2, xvals2D.data());
gridDefYvals(gridID2, yvals2D.data());
if (lproj_rll)
{
gridInqParamRLL(gridID1, &xpole, &ypole, &angle);
gridDefProj(gridID2, gridID1);
if (lbounds)
{
std::vector<double> xbounds, ybounds;
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
xvals2D[j * nx + i] = lamrot_to_lam(yvals[j], xvals[i], ypole, xpole, angle);
yvals2D[j * nx + i] = phirot_to_phi(yvals[j], xvals[i], ypole, angle);
}
if (nvertex == 2 && gridInqXbounds(gridID1, nullptr))
{
xbounds.resize(2 * nx);
gridInqXbounds(gridID1, xbounds.data());
if (check_range(2 * nx, xbounds.data(), -720, 720))
{
cdoWarning("longitude bounds out of range, skipped!");
xbounds.clear();
}
}
else
else if (nx > 1)
{
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
xvals2D[j * nx + i] = xscale * xvals[i];
yvals2D[j * nx + i] = yscale * yvals[j];
}
if (lproj_sinu)
cdo_sinu_to_lonlat(gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_laea)
cdo_laea_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_lcc)
cdo_lcc_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_stere)
cdo_stere_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj4)
cdo_proj_to_lonlat(proj4param, gridsize, xvals2D.data(), yvals2D.data());
xbounds.resize(2 * nx);
grid_gen_bounds(nx, xvals, xbounds);
}
gridDefXvals(gridID2, xvals2D.data());
gridDefYvals(gridID2, yvals2D.data());
if (lbounds)
if (nvertex == 2 && gridInqYbounds(gridID1, nullptr))
{
std::vector<double> xbounds, ybounds;
if (nvertex == 2 && gridInqXbounds(gridID1, nullptr))
ybounds.resize(2 * ny);
gridInqYbounds(gridID1, ybounds.data());
if (check_range(2 * ny, ybounds.data(), -180, 180))
{
xbounds.resize(2 * nx);
gridInqXbounds(gridID1, xbounds.data());
if (check_range(2 * nx, xbounds.data(), -720, 720))
{
cdoWarning("longitude bounds out of range, skipped!");
xbounds.clear();
}
cdoWarning("latitude bounds out of range, skipped!");
ybounds.clear();
}
else if (nx > 1)
}
else if (ny > 1)
{
ybounds.resize(2 * ny);
if (lproj_sinu || lproj_laea || lproj_lcc || lproj_stere || lproj4)
grid_gen_bounds(ny, yvals, ybounds);
else
{
xbounds.resize(2 * nx);
grid_gen_bounds(nx, xvals.data(), xbounds.data());
grid_gen_bounds(ny, yvals, ybounds);
grid_check_lat_borders(2 * ny, ybounds.data());
}
}
if (xbounds.size() && ybounds.size())
{
std::vector<double> xbounds2D(4 * gridsize);
std::vector<double> ybounds2D(4 * gridsize);
if (nvertex == 2 && gridInqYbounds(gridID1, nullptr))
if (lproj_rll)
{
ybounds.resize(2 * ny);
gridInqYbounds(gridID1, ybounds.data());
if (check_range(2 * ny, ybounds.data(), -180, 180))
{
cdoWarning("latitude bounds out of range, skipped!");
ybounds.clear();
}
gridGenRotBounds(xpole, ypole, angle, nx, ny, xbounds, ybounds, xbounds2D, ybounds2D);
}
else if (ny > 1)
else
{
ybounds.resize(2 * ny);
if (lproj_sinu || lproj_laea || lproj_lcc || lproj_stere || lproj4)
grid_gen_bounds(ny, yvals.data(), ybounds.data());
else
{
grid_gen_bounds(ny, yvals.data(), ybounds.data());
grid_check_lat_borders(2 * ny, ybounds.data());
}
}
if (xbounds.size() && ybounds.size())
{
std::vector<double> xbounds2D(4 * gridsize);
std::vector<double> ybounds2D(4 * gridsize);
if (lproj_rll)
{
gridGenRotBounds(xpole, ypole, angle, nx, ny, xbounds.data(), ybounds.data(), xbounds2D.data(), ybounds2D.data());
bounds_1D_to_2D(nx, ny, xbounds, ybounds, xbounds2D, ybounds2D, xscale, yscale);
if (lproj_sinu)
cdo_sinu_to_lonlat(4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_laea)
cdo_laea_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_lcc)
cdo_lcc_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_stere)
cdo_stere_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj4)
cdo_proj_to_lonlat(proj4param, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
}
else
{
if (lproj_sinu || lproj_laea || lproj_lcc || lproj_stere || lproj4)
{
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
index = j * 4 * nx + 4 * i;
xbounds2D[index + 0] = xscale * xbounds[2 * i];
ybounds2D[index + 0] = yscale * ybounds[2 * j];
xbounds2D[index + 1] = xscale * xbounds[2 * i];
ybounds2D[index + 1] = yscale * ybounds[2 * j + 1];
xbounds2D[index + 2] = xscale * xbounds[2 * i + 1];
ybounds2D[index + 2] = yscale * ybounds[2 * j + 1];
xbounds2D[index + 3] = xscale * xbounds[2 * i + 1];
ybounds2D[index + 3] = yscale * ybounds[2 * j];
}
if (lproj_sinu)
cdo_sinu_to_lonlat(4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_laea)
cdo_laea_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_lcc)
cdo_lcc_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_stere)
cdo_stere_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj4)
cdo_proj_to_lonlat(proj4param, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
}
else
{
grid_gen_xbounds2D(nx, ny, xbounds.data(), xbounds2D.data());
grid_gen_ybounds2D(nx, ny, ybounds.data(), ybounds2D.data());
}
grid_gen_xbounds2D(nx, ny, xbounds, xbounds2D);
grid_gen_ybounds2D(nx, ny, ybounds, ybounds2D);
}
gridDefXbounds(gridID2, xbounds2D.data());
gridDefYbounds(gridID2, ybounds2D.data());
}
}
gridCopyMask(gridID1, gridID2, gridsize);
}
else
{
cdoAbort("Grid type >%s< unsupported!", gridNamePtr(gridtype));
gridDefXbounds(gridID2, xbounds2D.data());
gridDefYbounds(gridID2, ybounds2D.data());
}
}
gridCopyMask(gridID1, gridID2, gridsize);
if (proj4param) Free(proj4param);
return gridID2;
......@@ -1341,7 +1342,7 @@ gridToUnstructured(int gridID1, int lbounds)
else if (nx > 1)
{
xbounds.resize(2 * nx);
grid_gen_bounds(nx, xvals.data(), xbounds.data());
grid_gen_bounds(nx, xvals, xbounds);
}
if (nvertex == 2 && gridInqYbounds(gridID1, nullptr))
......@@ -1352,7 +1353,7 @@ gridToUnstructured(int gridID1, int lbounds)
else if (ny > 1)
{
ybounds.resize(2 * ny);
grid_gen_bounds(ny, yvals.data(), ybounds.data());
grid_gen_bounds(ny, yvals, ybounds);
grid_check_lat_borders(2 * ny, ybounds.data());
}
......@@ -1362,13 +1363,12 @@ gridToUnstructured(int gridID1, int lbounds)
if (lproj_rll)
{
gridGenRotBounds(xpole, ypole, angle, nx, ny, xbounds.data(), ybounds.data(), xbounds2D.data(),
ybounds2D.data());
gridGenRotBounds(xpole, ypole, angle, nx, ny, xbounds, ybounds, xbounds2D, ybounds2D);
}
else
{
grid_gen_xbounds2D(nx, ny, xbounds.data(), xbounds2D.data());
grid_gen_ybounds2D(nx, ny, ybounds.data(), ybounds2D.data());
grid_gen_xbounds2D(nx, ny, xbounds, xbounds2D);
grid_gen_ybounds2D(nx, ny, ybounds, ybounds2D);
}
gridDefXbounds(gridID2, xbounds2D.data());
......@@ -1422,19 +1422,7 @@ gridToUnstructured(int gridID1, int lbounds)
if (lbounds)
{
size_t nvertex = (size_t) gridInqNvertex(gridID1);
std::vector<double> xbounds, ybounds;
/*
if (nvertex == 2 && gridInqXbounds(gridID1, nullptr))
{
xbounds.resize(2 * nlon);
gridInqXbounds(gridID1, xbounds.data());
}
else if (nlon > 1)
{
xbounds.resize(2 * nlon);
grid_gen_bounds(nx, xvals.data(), xbounds.data());
}
*/
std::vector<double> ybounds;
if (nvertex == 2 && gridInqYbounds(gridID1, nullptr))
{
ybounds.resize(2 * nlat);
......
......@@ -21,18 +21,11 @@
#include "grid_rot.h"
#include "grid_convert.h"
#include <stdio.h>
#include <cstdio>
#include <vector>
#ifdef __cplusplus
extern "C"
{
#endif
extern int qu2reg3_double(double *pfield, int *kpoint, int klat, int klon, double msval, int *kret, int omisng, int operio,
int oveggy);
#ifdef __cplusplus
}
#endif
extern "C" int qu2reg3_double(double *pfield, int *kpoint, int klat, int klon, double msval, int *kret, int omisng,
int operio, int oveggy);
extern bool gridVerbose;
extern char *gridSearchDir;
......@@ -59,11 +52,11 @@ void grid_to_radian(const char *units, size_t nvals, double *values, const char
void grid_to_degree(const char *units, size_t nvals, double *values, const char *description);
void grid_gen_corners(size_t n, const double *vals, double *corners);
void grid_gen_bounds(size_t n, const double *vals, double *bounds);
void grid_gen_bounds(size_t n, const std::vector<double> &vals, std::vector<double> &bounds);
void grid_check_lat_borders(int n, double *ybounds);
void grid_gen_xbounds2D(size_t nx, size_t ny, const double *xbounds, double *xbounds2D);
void grid_gen_ybounds2D(size_t nx, size_t ny, const double *ybounds, double *ybounds2D);
void grid_gen_xbounds2D(size_t nx, size_t ny, const std::vector<double> &xbounds, std::vector<double> &xbounds2D);
void grid_gen_ybounds2D(size_t nx, size_t ny, const std::vector<double> &ybounds, std::vector<double> &ybounds2D);
void grid_cell_center_to_bounds_X2D(const char *xunitstr, size_t xsize, size_t ysize, const double *grid_center_lon,
double *grid_corner_lon, double dlon);
......
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