Commit f31a5421 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

fieldzon: refactor and added support for GRID_GAUSSIAN_REDUCED.

parent a29e625a
......@@ -39,17 +39,9 @@
#include <mpim_grid.h>
#include "percentiles.h"
void *
Zonstat(void *process)
static void
addOperators(void)
{
int gridID1 = -1, gridID2 = -1;
int zongridID = -1;
int index;
int nrecs;
int varID, levelID;
cdoInitialize(process);
// clang-format off
cdoOperatorAdd("zonmin", func_min, 0, nullptr);
cdoOperatorAdd("zonmax", func_max, 0, nullptr);
......@@ -63,6 +55,20 @@ Zonstat(void *process)
cdoOperatorAdd("zonstd1", func_std1, 0, nullptr);
cdoOperatorAdd("zonpctl", func_pctl, 0, nullptr);
// clang-format on
}
void *
Zonstat(void *process)
{
int gridID1 = -1, gridID2 = -1;
int zongridID = -1;
int index;
int nrecs;
int varID, levelID;
cdoInitialize(process);
addOperators();
const int operatorID = cdoOperatorID();
const int operfunc = cdoOperatorF1(operatorID);
......
......@@ -19,30 +19,40 @@
#include "cdo_int.h"
#include "percentiles.h"
static
size_t fillRowlon(int gridID, size_t ny, std::vector<int> &rowlon, std::vector<int> &cumrowlon)
{
rowlon.resize(ny);
cumrowlon.resize(ny);
gridInqRowlon(gridID, rowlon.data());
cumrowlon[0] = 0;
for (size_t j = 1; j < ny; j++) cumrowlon[j] = cumrowlon[j-1] + rowlon[j-1];
size_t nx = rowlon[0];
for (size_t j = 1; j < ny; j++) if (rowlon[j] > (int)nx) nx = rowlon[j];
return nx;
}
void
zonmin(const Field &field1, Field &field2)
{
size_t rnmiss = 0;
double rmin;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (field1.nmiss)
{
rmin = arrayMinMV(nx, v, field1.missval);
if (DBL_IS_EQUAL(rmin, field1.missval)) rnmiss++;
}
else
{
rmin = arrayMin(nx, v);
}
const double rmin = field1.nmiss ? arrayMinMV(nx, v, field1.missval) : arrayMin(nx, v);
if (DBL_IS_EQUAL(rmin, field1.missval)) rnmiss++;
field2.vec[j] = rmin;
}
......@@ -54,26 +64,23 @@ void
zonmax(const Field &field1, Field &field2)
{
size_t rnmiss = 0;
double rmax;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (field1.nmiss)
{
rmax = arrayMaxMV(nx, v, field1.missval);
if (DBL_IS_EQUAL(rmax, field1.missval)) rnmiss++;
}
else
{
rmax = arrayMax(nx, v);
}
const double rmax = field1.nmiss ? arrayMaxMV(nx, v, field1.missval) : arrayMax(nx, v);
if (DBL_IS_EQUAL(rmax, field1.missval)) rnmiss++;
field2.vec[j] = rmax;
}
......@@ -85,26 +92,23 @@ void
zonrange(const Field &field1, Field &field2)
{
size_t rnmiss = 0;
double range;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (field1.nmiss)
{
range = arrayRangeMV(nx, v, field1.missval);
if (DBL_IS_EQUAL(range, field1.missval)) rnmiss++;
}
else
{
range = arrayRange(nx, v);
}
const double range = field1.nmiss ? arrayRangeMV(nx, v, field1.missval) : arrayRange(nx, v);
if (DBL_IS_EQUAL(range, field1.missval)) rnmiss++;
field2.vec[j] = range;
}
......@@ -116,26 +120,25 @@ void
zonsum(const Field &field1, Field &field2)
{
size_t rnmiss = 0;
double rsum = 0;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (field1.nmiss)
{
rsum = arraySumMV(nx, v, field1.missval);
if (DBL_IS_EQUAL(rsum, field1.missval)) rnmiss++;
}
else
{
rsum = arraySum(nx, v);
}
const double rsum = field1.nmiss ? arraySumMV(nx, v, field1.missval) : arraySum(nx, v);
if (DBL_IS_EQUAL(rsum, field1.missval)) rnmiss++;
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
field2.vec[j] = rsum;
}
......@@ -147,46 +150,21 @@ void
zonmean(const Field &field1, Field &field2)
{
size_t rnmiss = 0;
double rmean = 0;
size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
if (lreduced)
{
rowlon.resize(ny);
cumrowlon.resize(ny);
gridInqRowlon(field1.grid, rowlon.data());
cumrowlon[0] = 0;
for (size_t j = 1; j < ny; j++) cumrowlon[j] = cumrowlon[j-1] + rowlon[j-1];
nx = rowlon[0];
for (size_t j = 1; j < ny; j++) if (rowlon[j] > (int)nx) nx = rowlon[j];
}
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
if (lreduced)
{
nx = rowlon[j];
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[cumrowlon[j] + i];
}
else
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
}
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (field1.nmiss)
{
rmean = arrayMeanMV(nx, v, field1.missval);
}
else
{
rmean = arrayMean(nx, v);
}
const double rmean = field1.nmiss ? arrayMeanMV(nx, v, field1.missval) : arrayMean(nx, v);
if (DBL_IS_EQUAL(rmean, field1.missval)) rnmiss++;
......@@ -200,25 +178,21 @@ void
zonavg(const Field &field1, Field &field2)
{
size_t rnmiss = 0;
double ravg = 0;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (field1.nmiss)
{
ravg = arrayAvgMV(nx, v, field1.missval);
}
else
{
ravg = arrayMean(nx, v);
}
const double ravg = field1.nmiss ? arrayAvgMV(nx, v, field1.missval) : arrayMean(nx, v);
if (DBL_IS_EQUAL(ravg, field1.missval)) rnmiss++;
......@@ -271,14 +245,18 @@ zonvar(const Field &field1, Field &field2)
double rsum = 0, rsumw = 0;
double rsumq = 0, rsumwq = 0;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; j++)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
prevarsum_zon(v, nx, nmiss, missval1, rsum, rsumw, rsumq, rsumwq);
......@@ -300,14 +278,18 @@ zonvar1(const Field &field1, Field &field2)
double rsum = 0, rsumw = 0;
double rsumq = 0, rsumwq = 0;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; j++)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
prevarsum_zon(v, nx, nmiss, missval1, rsum, rsumw, rsumq, rsumwq);
......@@ -363,8 +345,11 @@ zonpctl(Field &field1, Field &field2, int p)
{
size_t rnmiss = 0;
const double missval = field1.missval;
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
const bool lreduced = gridInqType(field1.grid) == GRID_GAUSSIAN_REDUCED;
std::vector<int> rowlon, cumrowlon;
size_t nx = lreduced ? fillRowlon(field1.grid, ny, rowlon, cumrowlon) : gridInqXsize(field1.grid);
std::vector<double> v(nx);
......@@ -374,7 +359,9 @@ zonpctl(Field &field1, Field &field2, int p)
for (size_t j = 0; j < ny; j++)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
size_t l = 0;
for (size_t i = 0; i < nx; i++)
......@@ -395,7 +382,9 @@ zonpctl(Field &field1, Field &field2, int p)
{
for (size_t j = 0; j < ny; j++)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (lreduced) nx = rowlon[j];
const size_t offset = lreduced ? cumrowlon[j] : j * nx;
std::copy(field1.vec.begin()+offset, field1.vec.begin()+offset+nx, v.begin());
if (nx > 0)
{
......
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