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

Collgrid: Add support for variables on different grids.

parent a5d06bf5
No related branches found
No related tags found
No related merge requests found
Pipeline #8758 failed
......@@ -4,6 +4,11 @@
* Using spherepart and clipping from YAC2
* Version 2.0.0 release
2021-06-20 Uwe Schulzweida
* Collgrid: Add support for variables on different grids
* Collgrid: Add support for global cell indices from variable global_cell_indices
2021-06-16 Uwe Schulzweida
* Select: Add parameter levrange (level range)
......
......@@ -33,8 +33,6 @@ static int globalGridType = CDI_UNDEFID;
struct GridInfo
{
int globalIndicesID;
bool hasCenter;
bool hasBounds;
bool needed;
};
......@@ -44,7 +42,7 @@ struct CollgridInfo
int vlistID;
VarList varList;
Field field;
std::vector<long> cellIndex;
std::vector<std::vector<long>> cellIndex;
};
struct xyinfoType
......@@ -72,7 +70,7 @@ cmpxy_gt(const xyinfoType &a, const xyinfoType &b)
}
static int
gen_coll_grid(int ngrids, int nfiles, std::vector<CollgridInfo> &collgridInfo, bool initGrid, int igrid, long nxblocks)
gen_coll_grid(int ngrids, int nfiles, std::vector<CollgridInfo> &collgridInfo, int gindex, long nxblocks)
{
auto lsouthnorth = true;
auto lregular = false;
......@@ -80,7 +78,7 @@ gen_coll_grid(int ngrids, int nfiles, std::vector<CollgridInfo> &collgridInfo, b
long nx = (nxblocks != -1) ? nxblocks : -1;
auto gridID = vlistGrid(collgridInfo[0].vlistID, igrid);
auto gridID = vlistGrid(collgridInfo[0].vlistID, gindex);
const auto gridtype0 = (globalGridType != CDI_UNDEFID) ? globalGridType : gridInqType(gridID);
if (ngrids > 1 && gridtype0 == GRID_GENERIC && gridInqXsize(gridID) == 0 && gridInqYsize(gridID) == 0) return -1;
......@@ -96,7 +94,7 @@ gen_coll_grid(int ngrids, int nfiles, std::vector<CollgridInfo> &collgridInfo, b
for (int fileID = 0; fileID < nfiles; fileID++)
{
gridID = vlistGrid(collgridInfo[fileID].vlistID, igrid);
gridID = vlistGrid(collgridInfo[fileID].vlistID, gindex);
const auto gridtype = (globalGridType != CDI_UNDEFID) ? globalGridType : gridInqType(gridID);
if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || gridtype == GRID_PROJECTION)
lregular = true;
......@@ -231,31 +229,28 @@ gen_coll_grid(int ngrids, int nfiles, std::vector<CollgridInfo> &collgridInfo, b
yoff[j + 1] = yoff[j] + ysize[idx];
}
if (initGrid)
for (int fileID = 0; fileID < nfiles; fileID++)
{
for (int fileID = 0; fileID < nfiles; fileID++)
{
const long idx = xyinfo[fileID].id;
const long iy = fileID / nx;
const long ix = fileID - iy * nx;
const long offset = yoff[iy] * xsize2 + xoff[ix];
// printf("fileID %d %d, iy %d, ix %d, offset %d\n", fileID, xyinfo[fileID].id, iy, ix, offset);
long ij = 0;
for (long j = 0; j < ysize[idx]; ++j)
for (long i = 0; i < xsize[idx]; ++i)
const long idx = xyinfo[fileID].id;
const long iy = fileID / nx;
const long ix = fileID - iy * nx;
const long offset = yoff[iy] * xsize2 + xoff[ix];
// printf("fileID %d %d, iy %d, ix %d, offset %d\n", fileID, xyinfo[fileID].id, iy, ix, offset);
long ij = 0;
for (long j = 0; j < ysize[idx]; ++j)
for (long i = 0; i < xsize[idx]; ++i)
{
if (lcurvilinear || lunstructured)
{
if (lcurvilinear || lunstructured)
{
if (lcenter) xvals2[offset + j * xsize2 + i] = xvals[idx][ij];
if (lcenter) yvals2[offset + j * xsize2 + i] = yvals[idx][ij];
if (lbounds) for (long k = 0; k < nv; ++k) xbounds2[(offset + j * xsize2 + i) * nv + k] = xbounds[idx][ij * nv + k];
if (lbounds) for (long k = 0; k < nv; ++k) ybounds2[(offset + j * xsize2 + i) * nv + k] = ybounds[idx][ij * nv + k];
}
collgridInfo[idx].cellIndex[ij++] = offset + j * xsize2 + i;
if (lcenter) xvals2[offset + j * xsize2 + i] = xvals[idx][ij];
if (lcenter) yvals2[offset + j * xsize2 + i] = yvals[idx][ij];
if (lbounds) for (long k = 0; k < nv; ++k) xbounds2[(offset + j * xsize2 + i) * nv + k] = xbounds[idx][ij * nv + k];
if (lbounds) for (long k = 0; k < nv; ++k) ybounds2[(offset + j * xsize2 + i) * nv + k] = ybounds[idx][ij * nv + k];
}
}
collgridInfo[idx].cellIndex[gindex][ij++] = offset + j * xsize2 + i;
}
}
const auto gridID2 = gridCreate(gridtype0, xsize2 * ysize2);
......@@ -277,7 +272,7 @@ gen_coll_grid(int ngrids, int nfiles, std::vector<CollgridInfo> &collgridInfo, b
if (lbounds) gridDefYbounds(gridID2, ybounds2.data());
}
gridID = vlistGrid(collgridInfo[0].vlistID, igrid);
gridID = vlistGrid(collgridInfo[0].vlistID, gindex);
grid_copy_keys(gridID, gridID2);
......@@ -512,43 +507,36 @@ Collgrid(void *process)
// if ( Options::cdoVerbose ) vlistPrint(vlistID1);
// if ( Options::cdoVerbose ) vlistPrint(vlistID2);
const auto ngrids1 = vlistNgrids(vlistID1);
const auto ngrids2 = vlistNgrids(vlistID2);
for (int fileID = 0; fileID < nfiles; fileID++)
{
const auto cellIndexSize = vlistGridsizeMax(collgridInfo[fileID].vlistID);
collgridInfo[fileID].cellIndex.resize(cellIndexSize);
collgridInfo[fileID].cellIndex.resize(ngrids1);
for (int gindex = 0; gindex < ngrids1; ++gindex)
{
if (gridInfo[gindex].needed)
{
const auto patchSize = gridInqSize(vlistGrid(collgridInfo[fileID].vlistID, gindex));
collgridInfo[fileID].cellIndex[gindex].resize(patchSize);
}
}
}
const auto ngrids1 = vlistNgrids(vlistID1);
const auto ngrids2 = vlistNgrids(vlistID2);
std::vector<int> gridID2s(ngrids2);
auto initGrid = true;
for (int i2 = 0; i2 < ngrids2; ++i2)
{
int i1;
for (i1 = 0; i1 < ngrids1; ++i1)
if (vlistGrid(vlistID1, i1) == vlistGrid(vlistID2, i2)) break;
if (initGrid)
{
gridID2s[i2] = gen_coll_grid(ngrids2, nfiles, collgridInfo, initGrid, i1, nxblocks);
if (gridID2s[i2] != -1) initGrid = false;
}
else
gridID2s[i2] = gen_coll_grid(ngrids2, nfiles, collgridInfo, initGrid, i1, nxblocks);
gridID2s[i2] = gen_coll_grid(ngrids2, nfiles, collgridInfo, i1, nxblocks);
}
size_t gridsize2 = 0;
for (int i = 0; i < ngrids2; ++i)
{
if (gridID2s[i] != -1)
{
if (gridsize2 == 0) gridsize2 = gridInqSize(gridID2s[i]);
if (gridsize2 != gridInqSize(gridID2s[i]))
cdo_abort("Target gridsize differ (gridID=%d gridsize=%zu)!", i + 1, gridInqSize(gridID2s[i]));
vlistChangeGridIndex(vlistID2, i, gridID2s[i]);
}
if (gridID2s[i] != -1) vlistChangeGridIndex(vlistID2, i, gridID2s[i]);
}
VarList varList2;
......@@ -561,9 +549,9 @@ Collgrid(void *process)
const auto gridID = varList2[varID].gridID;
for (int i = 0; i < ngrids2; ++i)
{
if (gridID2s[i] != -1)
if (gridID2s[i] != -1 && gridID == vlistGrid(vlistID2, i))
{
if (gridID == vlistGrid(vlistID2, i) || gridsize2 == gridInqSize(gridID)) collectVars2[varID] = true;
collectVars2[varID] = true;
break;
}
}
......@@ -610,7 +598,7 @@ Collgrid(void *process)
size_t nmiss;
cdo_read_record(collgridInfo[fileID].streamID, cellIndex.data(), &nmiss);
for (size_t i = 0; i < patchSize; ++i)
collgridInfo[fileID].cellIndex[i] = std::lround(cellIndex[i]) - 1;
collgridInfo[fileID].cellIndex[gindex][i] = std::lround(cellIndex[i]) - 1;
}
}
......@@ -622,11 +610,12 @@ Collgrid(void *process)
field2.init(varList2[varID2]);
const auto missval = varList2[varID2].missval;
const auto gridsize = field2.size;
const auto missval = field2.missval;
if (field2.memType == MemType::Float)
for (size_t i = 0; i < gridsize2; i++) field2.vec_f[i] = missval;
for (size_t i = 0; i < gridsize; i++) field2.vec_f[i] = missval;
else
for (size_t i = 0; i < gridsize2; i++) field2.vec_d[i] = missval;
for (size_t i = 0; i < gridsize; i++) field2.vec_d[i] = missval;
#ifdef _OPENMP
#pragma omp parallel for default(shared)
......@@ -637,7 +626,7 @@ Collgrid(void *process)
field1.init(collgridInfo[fileID].varList[varID]);
cdo_read_record(collgridInfo[fileID].streamID, field1);
if (collectVars2[varID2]) collect_cells(field1, field2, collgridInfo[fileID].cellIndex);
if (collectVars2[varID2]) collect_cells(field1, field2, collgridInfo[fileID].cellIndex[gindex]);
}
cdo_def_record(streamID2, varID2, levelID2);
......@@ -648,7 +637,9 @@ Collgrid(void *process)
cdo_write_record(streamID2, field2);
}
else
cdo_write_record(streamID2, collgridInfo[0].field);
{
cdo_write_record(streamID2, collgridInfo[0].field);
}
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment