Skip to content
Snippets Groups Projects
Commit 96e5ae58 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

thinout: added support for curvilinear grids.

parent 221f5d83
No related branches found
No related tags found
No related merge requests found
Pipeline #20686 passed
Stage: build
Stage: check
Hide whitespace changes
Inline Side-by-side
Showing
with 85 additions and 60 deletions
src/Intgrid.cc
+ 85
− 60
View file @ 96e5ae58
......@@ -21,8 +21,48 @@
#include "griddes.h"
#include "matrix_view.h"
template <typename T>
static void
thinout(const Varray<T> &varray1, Varray<T> &varray2, int gridID1, int gridID2, size_t xinc, size_t yinc)
{
const auto nlon1 = gridInqXsize(gridID1);
const auto nlat1 = gridInqYsize(gridID1);
const auto nlon2 = gridInqXsize(gridID2);
const auto nlat2 = gridInqYsize(gridID2);
MatrixView<const T> xfield1(varray1.data(), nlat1, nlon1);
MatrixView<T> xfield2(varray2.data(), nlat2, nlon2);
size_t olat = 0;
for (size_t ilat = 0; ilat < nlat1; ilat += yinc)
{
size_t olon = 0;
for (size_t ilon = 0; ilon < nlon1; ilon += xinc)
{
xfield2[olat][olon] = xfield1[ilat][ilon];
olon++;
}
olat++;
}
}
static void
thinout(const Field &field1, Field &field2, size_t xinc, size_t yinc)
{
if (field1.memType != field2.memType) cdo_abort("Interal error, memType of field1 and field2 differ!");
if (field1.memType == MemType::Float)
thinout(field1.vec_f, field2.vec_f, field1.grid, field2.grid, xinc, yinc);
else
thinout(field1.vec_d, field2.vec_d, field1.grid, field2.grid, xinc, yinc);
field_num_mv(field2);
}
static int
genThinoutGrid(int gridID1, size_t xinc, size_t yinc)
gen_thinout_grid(int gridID1, size_t xinc, size_t yinc)
{
const auto nlon1 = gridInqXsize(gridID1);
const auto nlat1 = gridInqYsize(gridID1);
......@@ -33,17 +73,19 @@ genThinoutGrid(int gridID1, size_t xinc, size_t yinc)
if (nlat1 % yinc) nlat2++;
const auto gridsize2 = nlon2 * nlat2;
const auto gridID2 = gridCreate(GRID_LONLAT, gridsize2);
const auto gridtype = gridInqType(gridID1);
const auto gridID2 = gridCreate(gridtype, gridsize2);
gridDefXsize(gridID2, nlon2);
gridDefYsize(gridID2, nlat2);
const auto gridtype = gridInqType(gridID1);
grid_copy_keys(gridID1, gridID2);
if (gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT)
{
Varray<double> xvals1(nlon1), yvals1(nlat1);
Varray<double> xvals2(nlon2), yvals2(nlat2);
gridInqXvals(gridID1, &xvals1[0]);
gridInqYvals(gridID1, &yvals1[0]);
gridInqXvals(gridID1, xvals1.data());
gridInqYvals(gridID1, yvals1.data());
size_t olat = 0;
for (size_t ilat = 0; ilat < nlat1; ilat += yinc) yvals2[olat++] = yvals1[ilat];
......@@ -51,8 +93,21 @@ genThinoutGrid(int gridID1, size_t xinc, size_t yinc)
size_t olon = 0;
for (size_t ilon = 0; ilon < nlon1; ilon += xinc) xvals2[olon++] = xvals1[ilon];
gridDefXvals(gridID2, &xvals2[0]);
gridDefYvals(gridID2, &yvals2[0]);
gridDefXvals(gridID2, xvals2.data());
gridDefYvals(gridID2, yvals2.data());
}
else if (gridtype == GRID_CURVILINEAR)
{
Varray<double> xvals1(nlon1*nlat1), yvals1(nlon1*nlat1);
Varray<double> xvals2(nlon2*nlat2), yvals2(nlon2*nlat2);
gridInqXvals(gridID1, xvals1.data());
gridInqYvals(gridID1, yvals1.data());
thinout(xvals1, xvals2, gridID1, gridID2, xinc, yinc);
thinout(yvals1, yvals2, gridID1, gridID2, xinc, yinc);
gridDefXvals(gridID2, xvals2.data());
gridDefYvals(gridID2, yvals2.data());
}
else
{
......@@ -63,7 +118,7 @@ genThinoutGrid(int gridID1, size_t xinc, size_t yinc)
}
static int
genBoxavgGrid(int gridID1, size_t xinc, size_t yinc)
gen_boxavg_grid(int gridID1, size_t xinc, size_t yinc)
{
const auto nlon1 = gridInqXsize(gridID1);
const auto nlat1 = gridInqYsize(gridID1);
......@@ -192,49 +247,10 @@ boxavg(const Field &field1, Field &field2, size_t xinc, size_t yinc)
field_num_mv(field2);
}
template <typename T>
static void
thinout(const Varray<T> &varray1, Varray<T> &varray2, int gridID1, int gridID2, size_t xinc, size_t yinc)
{
const auto nlon1 = gridInqXsize(gridID1);
const auto nlat1 = gridInqYsize(gridID1);
const auto nlon2 = gridInqXsize(gridID2);
const auto nlat2 = gridInqYsize(gridID2);
MatrixView<const T> xfield1(varray1.data(), nlat1, nlon1);
MatrixView<T> xfield2(varray2.data(), nlat2, nlon2);
size_t olat = 0;
for (size_t ilat = 0; ilat < nlat1; ilat += yinc)
{
size_t olon = 0;
for (size_t ilon = 0; ilon < nlon1; ilon += xinc)
{
xfield2[olat][olon] = xfield1[ilat][ilon];
olon++;
}
olat++;
}
}
static void
thinout(const Field &field1, Field &field2, size_t xinc, size_t yinc)
{
if (field1.memType != field2.memType) cdo_abort("Interal error, memType of field1 and field2 differ!");
if (field1.memType == MemType::Float)
thinout(field1.vec_f, field2.vec_f, field1.grid, field2.grid, xinc, yinc);
else
thinout(field1.vec_d, field2.vec_d, field1.grid, field2.grid, xinc, yinc);
field_num_mv(field2);
}
void *
Intgrid(void *process)
{
int gridID1 = -1, gridID2 = -1;
int gridID2 = -1;
int xinc = 1, yinc = 1;
cdo_initialize(process);
......@@ -277,30 +293,42 @@ Intgrid(void *process)
const auto ngrids = vlistNgrids(vlistID1);
for (int index = 0; index < ngrids; ++index)
{
gridID1 = vlistGrid(vlistID1, index);
const auto gridID1 = vlistGrid(vlistID1, index);
const auto gridtype = gridInqType(gridID1);
if (operatorID == BOXAVG || operatorID == THINOUT)
if (operatorID == BOXAVG)
{
if (index == 0)
{
if (gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN)
cdo_abort("Interpolation of %s data unsupported!", gridNamePtr(gridtype));
gridID2 = gen_boxavg_grid(gridID1, xinc, yinc);
}
else
cdo_abort("Too many different grids!");
}
if (operatorID == THINOUT)
{
if (index == 0)
{
if (gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN /* && gridtype != GRID_CURVILINEAR */)
if (gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN && gridtype != GRID_CURVILINEAR)
cdo_abort("Interpolation of %s data unsupported!", gridNamePtr(gridtype));
gridID2 = (operatorID == BOXAVG) ? genBoxavgGrid(gridID1, xinc, yinc) : genThinoutGrid(gridID1, xinc, yinc);
gridID2 = gen_thinout_grid(gridID1, xinc, yinc);
}
else
cdo_abort("Too many different grids!");
}
else if (operatorID == INTGRIDBIL || operatorID == INTERPOLATE)
{
const bool ldistgen = (grid_is_distance_generic(gridID1) && grid_is_distance_generic(gridID2));
const auto ldistgen = (grid_is_distance_generic(gridID1) && grid_is_distance_generic(gridID2));
if (!ldistgen && gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN)
cdo_abort("Interpolation of %s data unsupported!", gridNamePtr(gridtype));
}
else if (operatorID == INTGRIDNN || operatorID == INTGRIDDIS)
{
const bool lprojparams = (gridtype == GRID_PROJECTION) && grid_has_proj_params(gridID1);
const auto lprojparams = (gridtype == GRID_PROJECTION) && grid_has_proj_params(gridID1);
if (!gridProjIsSupported(gridID1) && !lprojparams && gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN
&& gridtype != GRID_GME && gridtype != GRID_CURVILINEAR && gridtype != GRID_UNSTRUCTURED)
cdo_abort("Interpolation of %s data unsupported!", gridNamePtr(gridtype));
......@@ -347,13 +375,10 @@ Intgrid(void *process)
{
cdo_read_record(streamID1, field1.vec_d.data(), &field1.nmiss);
gridID1 = varList1[varID].gridID;
const auto missval = varList1[varID].missval;
field1.grid = gridID1;
field1.missval = missval;
field1.grid = varList1[varID].gridID;
field1.missval = varList1[varID].missval;
field2.grid = gridID2;
field2.missval = missval;
field2.missval = field1.missval;
field2.nmiss = 0;
}
......
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

Imprint | Privacy Policy