Commit 2bff97c2 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added function get_projection().

parent 0804ebf5
......@@ -988,6 +988,38 @@ bounds_1D_to_2D(size_t nx, size_t ny, const std::vector<double> &xbounds, const
}
}
enum struct Projection
{
none,
proj4,
proj_rll,
proj_laea,
proj_lcc,
proj_sinu,
proj_stere
};
Projection get_projection(int gridID1)
{
Projection projection(Projection::none);
const int projtype = gridInqProjType(gridID1);
// clang-format off
if (projtype == CDI_PROJ_RLL) projection = Projection::proj_rll;
else if (projtype == CDI_PROJ_LAEA) projection = Projection::proj_laea;
else if (projtype == CDI_PROJ_LCC) projection = Projection::proj_lcc;
else if (projtype == CDI_PROJ_SINU) projection = Projection::proj_sinu;
else if (projtype == CDI_PROJ_STERE) projection = Projection::proj_stere;
else
{
char mapname[CDI_MAX_NAME];
cdiGridInqKeyStr(gridID1, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapname);
cdoAbort("Projection type >%s< unsupported!", mapname);
}
// clang-format on
return projection;
}
int
gridToCurvilinear(int gridID1, int lbounds)
{
......@@ -1005,42 +1037,18 @@ gridToCurvilinear(int gridID1, int lbounds)
const int gridID2 = gridCreate(GRID_CURVILINEAR, gridsize);
gridDefDatatype(gridID2, CDI_DATATYPE_FLT32);
Projection projection(Projection::none);
char *proj4param = nullptr;
bool lproj4 = false;
bool lproj_rll = false;
bool lproj_laea = false;
bool lproj_lcc = false;
bool lproj_sinu = false;
bool lproj_stere = false;
if (gridtype == GRID_PROJECTION && gridsize == nx * ny)
{
proj4param = grid_get_proj4param(gridID1);
if (proj4param)
{
lproj4 = true;
gridtype = GRID_LONLAT;
}
else
{
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;
else if (projtype == CDI_PROJ_LCC) lproj_lcc = true;
else if (projtype == CDI_PROJ_SINU) lproj_sinu = true;
else if (projtype == CDI_PROJ_STERE) lproj_stere = true;
else
{
char mapname[CDI_MAX_NAME];
cdiGridInqKeyStr(gridID1, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapname);
cdoAbort("Projection type >%s< unsupported!", mapname);
}
// clang-format on
if (lproj_rll || lproj_laea || lproj_lcc || lproj_sinu || lproj_stere) gridtype = GRID_LONLAT;
}
projection = proj4param ? Projection::proj4 : get_projection(gridID1);
if (projection != Projection::none) gridtype = GRID_LONLAT;
}
const bool lprojection = projection == Projection::proj_laea || projection == Projection::proj_lcc ||
projection == Projection::proj_sinu || projection == Projection::proj_stere;
double xpole = 0, ypole = 0, angle = 0;
double xscale = 1, yscale = 1;
char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];
......@@ -1049,7 +1057,7 @@ gridToCurvilinear(int gridID1, int lbounds)
gridInqXunits(gridID1, xunits);
gridInqYunits(gridID1, yunits);
if (lproj_laea || lproj_lcc || lproj_sinu || lproj_stere)
if (lprojection)
{
check_units("x", xunits);
check_units("y", yunits);
......@@ -1061,22 +1069,19 @@ gridToCurvilinear(int gridID1, int lbounds)
gridDefXsize(gridID2, nx);
gridDefYsize(gridID2, ny);
std::vector<double> xvals2D(gridsize);
std::vector<double> yvals2D(gridsize);
std::vector<double> xvals2D(gridsize), yvals2D(gridsize);
if (nx == 0) nx = 1;
if (ny == 0) ny = 1;
std::vector<double> xvals(nx, 0);
std::vector<double> yvals(ny, 0);
std::vector<double> xvals(nx, 0), yvals(ny, 0);
if (gridInqXvals(gridID1, nullptr)) gridInqXvals(gridID1, xvals.data());
if (gridInqYvals(gridID1, nullptr)) gridInqYvals(gridID1, yvals.data());
if (lproj_rll)
if (projection == Projection::proj_rll)
{
gridInqParamRLL(gridID1, &xpole, &ypole, &angle);
gridDefProj(gridID2, gridID1);
gridInqParamRLL(gridID1, &xpole, &ypole, &angle);
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
......@@ -1089,15 +1094,15 @@ gridToCurvilinear(int gridID1, int lbounds)
{
center_1D_to_2D(nx, ny, xvals, yvals, xvals2D, yvals2D, xscale, yscale);
if (lproj_sinu)
if (projection == Projection::proj_sinu)
cdo_sinu_to_lonlat(gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_laea)
else if (projection == Projection::proj_laea)
cdo_laea_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_lcc)
else if (projection == Projection::proj_lcc)
cdo_lcc_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj_stere)
else if (projection == Projection::proj_stere)
cdo_stere_to_lonlat(gridID1, gridsize, xvals2D.data(), yvals2D.data());
else if (lproj4)
else if (projection == Projection::proj4)
cdo_proj_to_lonlat(proj4param, gridsize, xvals2D.data(), yvals2D.data());
}
......@@ -1137,7 +1142,7 @@ gridToCurvilinear(int gridID1, int lbounds)
else if (ny > 1)
{
ybounds.resize(2 * ny);
if (lproj_sinu || lproj_laea || lproj_lcc || lproj_stere || lproj4)
if (lprojection || projection == Projection::proj4)
grid_gen_bounds(ny, yvals, ybounds);
else
{
......@@ -1148,28 +1153,27 @@ gridToCurvilinear(int gridID1, int lbounds)
if (xbounds.size() && ybounds.size())
{
std::vector<double> xbounds2D(4 * gridsize);
std::vector<double> ybounds2D(4 * gridsize);
std::vector<double> xbounds2D(4 * gridsize), ybounds2D(4 * gridsize);
if (lproj_rll)
if (projection == Projection::proj_rll)
{
gridGenRotBounds(xpole, ypole, angle, nx, ny, xbounds, ybounds, xbounds2D, ybounds2D);
}
else
{
if (lproj_sinu || lproj_laea || lproj_lcc || lproj_stere || lproj4)
if (lprojection || projection == Projection::proj4)
{
bounds_1D_to_2D(nx, ny, xbounds, ybounds, xbounds2D, ybounds2D, xscale, yscale);
if (lproj_sinu)
if (projection == Projection::proj_sinu)
cdo_sinu_to_lonlat(4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_laea)
else if (projection == Projection::proj_laea)
cdo_laea_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_lcc)
else if (projection == Projection::proj_lcc)
cdo_lcc_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj_stere)
else if (projection == Projection::proj_stere)
cdo_stere_to_lonlat(gridID1, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
else if (lproj4)
else if (projection == Projection::proj4)
cdo_proj_to_lonlat(proj4param, 4 * gridsize, xbounds2D.data(), ybounds2D.data());
}
else
......@@ -1297,41 +1301,39 @@ gridToUnstructured(int gridID1, int lbounds)
gridDefXsize(gridID2, gridsize);
gridDefYsize(gridID2, gridsize);
std::vector<double> xvals(nx), yvals(ny);
std::vector<double> xvals2D(gridsize), yvals2D(gridsize);
{
std::vector<double> xvals2D(gridsize), yvals2D(gridsize);
gridInqXvals(gridID1, xvals.data());
gridInqYvals(gridID1, yvals.data());
std::vector<double> xvals(nx, 0), yvals(ny, 0);
if (gridInqXvals(gridID1, nullptr)) gridInqXvals(gridID1, xvals.data());
if (gridInqXvals(gridID1, nullptr)) gridInqYvals(gridID1, yvals.data());
if (lproj_rll)
{
gridInqParamRLL(gridID1, &xpole, &ypole, &angle);
if (lproj_rll)
{
gridInqParamRLL(gridID1, &xpole, &ypole, &angle);
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
{
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
xvals2D[j * nx + i] = xvals[i];
yvals2D[j * nx + i] = yvals[j];
}
}
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
{
for (size_t j = 0; j < ny; j++)
for (size_t i = 0; i < nx; i++)
{
xvals2D[j * nx + i] = xvals[i];
yvals2D[j * nx + i] = yvals[j];
}
}
gridDefXvals(gridID2, xvals2D.data());
gridDefYvals(gridID2, yvals2D.data());
}
gridDefXvals(gridID2, xvals2D.data());
gridDefYvals(gridID2, yvals2D.data());
if (lbounds)
{
size_t nvertex = (size_t) gridInqNvertex(gridID1);
const size_t nvertex = (size_t) gridInqNvertex(gridID1);
std::vector<double> xbounds, ybounds;
if (nvertex == 2 && gridInqXbounds(gridID1, nullptr))
......
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