Skip to content
GitLab
Commit 67c308a5 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Remapstat: test implementation for unstructured target grids.

parent 257074b2
Pipeline #4539 passed with stages
in 16 minutes and 44 seconds
Hide whitespace changes
Inline Side-by-side
src/Remapstat.cc
View file @ 67c308a5
......@@ -30,6 +30,7 @@
#include "grid_point_search.h"
#include <mpim_grid.h>
#include "cimdOmp.h"
#include "verifygrid.h"
//#define TESTIMPL 1
......@@ -71,15 +72,15 @@ struct StatInfo
static
size_t read_target_cell_bounds(int gridID, Varray<double> &xbounds, Varray<double> &ybounds)
{
if (!gridHasCoordinates(gridID)) cdoAbort("Target cell corner coordinates missing!");
if (!gridHasBounds(gridID)) cdoAbort("Target cell corner coordinates missing!");
const auto nv = gridInqNvertex(gridID);
const auto gridsize = gridInqSize(gridID);
xbounds.resize(nv * gridsize);
ybounds.resize(nv * gridsize);
gridInqXvals(gridID, xbounds.data());
gridInqYvals(gridID, xbounds.data());
gridInqXbounds(gridID, xbounds.data());
gridInqYbounds(gridID, ybounds.data());
// Convert lat/lon units if required
cdo_grid_to_radian(gridID, CDI_XAXIS, nv * gridsize, xbounds.data(), "grid corner lon");
......@@ -195,6 +196,7 @@ double calc_maxdist(size_t i, size_t nv, double plon, double plat, const Varray<
const auto sdist = squareDistance(p1, p2);
maxdist = std::max(maxdist, sdist);
}
maxdist = 1.01 * std::sqrt(maxdist);
return maxdist;
......@@ -232,16 +234,83 @@ double calc_maxdist_rec2d(size_t i, size_t nlon, double plon, double plat, const
}
static
size_t find_points(size_t i, size_t nv, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
size_t find_points(std::vector<char> &vmask, size_t cell_no, size_t ncorner, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
const Varray<double> &yvals1, const Varray<double> &xbounds, const Varray<double> &ybounds)
{
#ifndef TESTIMPL
cdoAbort("Internal error: find_points() not implemented!");
#endif
Point3D center_point_xyz;
Point center_point_plane_projection;
Varray<Point> cell_corners_plane_projection(ncorner + 1);
Varray<Point3D> cell_corners_xyz(ncorner + 1);
Varray<Point3D> cell_corners_xyz_open_cell(ncorner);
std::vector<bool> marked_duplicate_indices(ncorner);
set_cell_corners_xyz(ncorner, &xbounds[cell_no * ncorner], &ybounds[cell_no * ncorner], cell_corners_xyz_open_cell);
auto actual_number_of_corners = get_actual_number_of_corners(ncorner, cell_corners_xyz_open_cell);
auto no_duplicates = get_no_duplicates(actual_number_of_corners, cell_corners_xyz_open_cell, marked_duplicate_indices);
copy_unique_corners(actual_number_of_corners, cell_corners_xyz_open_cell, marked_duplicate_indices, cell_corners_xyz);
actual_number_of_corners -= no_duplicates;
cell_corners_xyz[actual_number_of_corners] = cell_corners_xyz[0];
if (actual_number_of_corners < 3) return 0;
const auto coordinate_to_ignore = find_coordinate_to_ignore(cell_corners_xyz);
bool invert_result = false;
const auto cval = (coordinate_to_ignore == 1) ? cell_corners_xyz[0].X : ((coordinate_to_ignore == 2) ? cell_corners_xyz[0].Y : cell_corners_xyz[0].Z);
if (cval < 0.0) invert_result = true;
set_cell_corners_plane_projection(coordinate_to_ignore, actual_number_of_corners, cell_corners_xyz, cell_corners_plane_projection);
const auto polygon_area = calculate_the_polygon_area(cell_corners_plane_projection, actual_number_of_corners + 1);
auto is_clockwise = are_polygon_vertices_arranged_in_clockwise_order(polygon_area);
if (invert_result) is_clockwise = !is_clockwise;
if (is_clockwise) return 0;
double center_coordinates[3];
size_t nvalues = 0;
for (size_t k = 0; k < nadds; ++k)
{
const auto index1 = adds[k];
const auto lon = xvals1[index1];
const auto lat = yvals1[index1];
gcLLtoXYZ(lon, lat, center_coordinates);
center_point_xyz.X = center_coordinates[0];
center_point_xyz.Y = center_coordinates[1];
center_point_xyz.Z = center_coordinates[2];
set_center_point_plane_projection(coordinate_to_ignore, center_point_xyz, center_point_plane_projection);
auto winding_number
= winding_numbers_algorithm(cell_corners_plane_projection, actual_number_of_corners + 1, center_point_plane_projection);
if (winding_number != 0)
{
if (vmask[index1] == 0)
{
adds[nvalues] = adds[k];
nvalues++;
}
if (vmask[index1] < 127) vmask[index1]++;
}
}
return nvalues;
}
static
size_t find_points_rec2d(size_t i, size_t nlon2, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
size_t find_points_rec2d(std::vector<char> &vmask, size_t i, size_t nlon2, size_t nadds, std::vector<size_t> &adds, const Varray<double> &xvals1,
const Varray<double> &yvals1, const Varray<double> &xbounds2d, const Varray<double> &ybounds2d)
{
const auto iy = i / nlon2;
......@@ -250,13 +319,15 @@ size_t find_points_rec2d(size_t i, size_t nlon2, size_t nadds, std::vector<size_
size_t nvalues = 0;
for (size_t k = 0; k < nadds; ++k)
{
const auto x = xvals1[adds[k]];
const auto y = yvals1[adds[k]];
const auto index1 = adds[k];
const auto x = xvals1[index1];
const auto y = yvals1[index1];
if (y >= ybounds2d[2*iy] && y < ybounds2d[2*iy + 1])
if ((x >= xbounds2d[2*ix] && x < xbounds2d[2*ix + 1]) ||
((x-PI2) >= xbounds2d[2*ix] && (x-PI2) < xbounds2d[2*ix + 1]))
{
adds[nvalues] = adds[k];
if (vmask[index1] < 127) vmask[index1]++;
adds[nvalues] = index1;
nvalues++;
}
}
......@@ -264,6 +335,19 @@ size_t find_points_rec2d(size_t i, size_t nlon2, size_t nadds, std::vector<size_
return nvalues;
}
static
void check_vmask(std::vector<char> vmask)
{
constexpr int max_vals = 128;
size_t vm[max_vals] = { 0 };
const auto size = vmask.size();
for (size_t i = 0; i < size; ++i) if (vmask[i] >= 0 && vmask[i] < max_vals) vm[(int)vmask[i]]++;
for (int i = 0; i < max_vals; ++i) if (vm[i]) cdoPrint("Number of source values: %d --> %zu/%zu", i, vm[i], size);
size_t sum = 0;
for (int i = 1; i < max_vals; ++i) sum += vm[i];
cdoPrint("Sum of used source values: %zu/%zu", sum, size);
}
static
Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
{
......@@ -322,7 +406,7 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
start = Options::cdoVerbose ? cdo_get_wtime() : 0;
size_t ndist = gridsize1;
std::vector<double> values(gridsize1);
std::vector<char> vmask(gridsize1, 0);
std::vector<size_t> adds(gridsize1);
std::vector<double> dist(gridsize1);
......@@ -341,8 +425,8 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
// printf("%zu nadds %zu\n", i+1, nadds);
const auto nvalues = grid2_is_reg2d
? find_points_rec2d(i, nlon2, nadds, adds, xvals1, yvals1, xbounds2d, ybounds2d)
: find_points(i, nv, nadds, adds, xvals1, yvals1, xbounds, ybounds);
? find_points_rec2d(vmask, i, nlon2, nadds, adds, xvals1, yvals1, xbounds2d, ybounds2d)
: find_points(vmask, i, nv, nadds, adds, xvals1, yvals1, xbounds, ybounds);
if (nvalues)
{
......@@ -356,6 +440,8 @@ Varray2D<size_t> gen_mapdata(int gridID1, int gridID2)
if (Options::cdoVerbose) statInfo.print();
if (Options::cdoVerbose) check_vmask(vmask);
if (Options::cdoVerbose) cdoPrint("Point search qnearest: %.2f seconds", cdo_get_wtime() - start);
gridPointSearchDelete(gps);
......@@ -468,14 +554,17 @@ Remapstat(void *process)
operatorInputArg("grid description file or name");
const auto gridID2 = cdoDefineGrid(cdoOperatorArgv(0));
auto gridtype2 = gridInqType(gridID2);
const bool lprojparams = (gridtype2 == GRID_PROJECTION) && grid_has_proj_params(gridID2);
{
#ifdef TESTIMPL
if (!gridProjIsSupported(gridID2) && !lprojparams && gridtype2 != GRID_LONLAT && gridtype2 != GRID_GAUSSIAN
&& gridtype2 != GRID_CURVILINEAR && gridtype2 != GRID_UNSTRUCTURED)
const bool lprojparams = (gridtype2 == GRID_PROJECTION) && grid_has_proj_params(gridID2);
if (!gridProjIsSupported(gridID2) && !lprojparams && gridtype2 != GRID_LONLAT && gridtype2 != GRID_GAUSSIAN
&& gridtype2 != GRID_CURVILINEAR && gridtype2 != GRID_UNSTRUCTURED)
#else
if (gridtype2 != GRID_LONLAT && gridtype2 != GRID_GAUSSIAN)
if (gridtype2 != GRID_LONLAT && gridtype2 != GRID_GAUSSIAN)
#endif
cdoAbort("Remapping to %s data unsupported!", gridNamePtr(gridtype2));
cdoAbort("Remapping to %s data unsupported!", gridNamePtr(gridtype2));
}
const auto streamID1 = cdoOpenRead(0);
......
src/Verifygrid.cc
View file @ 67c308a5
......@@ -204,23 +204,15 @@ winding_numbers_algorithm(const Varray<Point> &cell_corners, int number_corners,
{
if (cell_corners[k + 1].y > point.y)
{
const auto &point1 = cell_corners[k];
const auto &point2 = cell_corners[k + 1];
if (is_point_left_of_edge(point1, point2, point) > 0) winding_number++;
//printf(" 1: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
if (is_point_left_of_edge(cell_corners[k], cell_corners[k + 1], point) > 0) winding_number++;
}
}
else
{
if (cell_corners[k + 1].y <= point.y)
{
const auto &point1 = cell_corners[k];
const auto &point2 = cell_corners[k + 1];
if (is_point_left_of_edge(point1, point2, point) < 0) winding_number--;
//printf(" 2: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
}
if (is_point_left_of_edge(cell_corners[k], cell_corners[k + 1], point) < 0) winding_number--;
}
}
}
......@@ -691,21 +683,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
{
no_counterclockwise_cells += 1;
}
/*
printf("ncorner actual_number_of_corners %d %d\n", ncorner, actual_number_of_corners);
printf("cell_no, is_clockwise, polygon_area %zu %d, %g\n", cell_no, is_clockwise, polygon_area);
printf("bounds %g %g %g %g %g %g %g %g\n",
grid_corner_lon[cell_no * ncorner], grid_corner_lat[cell_no * ncorner],
grid_corner_lon[cell_no * ncorner+1], grid_corner_lat[cell_no * ncorner+1],
grid_corner_lon[cell_no * ncorner+2], grid_corner_lat[cell_no * ncorner+2],
grid_corner_lon[cell_no * ncorner+3], grid_corner_lat[cell_no * ncorner+3]);
printf("plane %g %g %g %g %g %g %g %g %g %g\n",
cell_corners_plane_projection[0].x, cell_corners_plane_projection[0].y,
cell_corners_plane_projection[1].x, cell_corners_plane_projection[1].y,
cell_corners_plane_projection[2].x, cell_corners_plane_projection[2].y,
cell_corners_plane_projection[3].x, cell_corners_plane_projection[3].y,
center_point_plane_projection.x, center_point_plane_projection.y);
*/
// The winding numbers algorithm is used to test whether the presumed center point is within the bounds of the cell.
auto winding_number = winding_numbers_algorithm(cell_corners_plane_projection, actual_number_of_corners + 1, center_point_plane_projection);
......
src/verifygrid.h
View file @ 67c308a5
......@@ -20,6 +20,11 @@
struct Point { double x, y; };
struct Point3D { double X, Y, Z; };
int get_actual_number_of_corners(int ncorner, const Varray<Point3D> &cell_corners_xyz_open_cell);
int get_no_duplicates(int actual_number_of_corners, const Varray<Point3D> &cell_corners_xyz_open_cell, std::vector<bool> &marked_duplicate_indices);
void copy_unique_corners(int actual_number_of_corners, const Varray<Point3D> &cell_corners_xyz_open_cell, std::vector<bool> &marked_duplicate_indices,
Varray<Point3D> &cell_corners_xyz_without_duplicates);
void set_cell_corners_xyz(int ncorner, const double *cell_corners_lon, const double *cell_corners_lat, Varray<Point3D> &cell_corners_xyz);
void set_center_point_plane_projection(int coordinate_to_ignore, const Point3D &center_point_xyz, Point &center_point_plane_projection);
void set_cell_corners_plane_projection(int coordinate_to_ignore, int ncorner, const Varray<Point3D> &cell_corners_xyz, Varray<Point> &cell_corners_plane_projection);
......
Supports Markdown
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