From 1738fbacca4665d76f5f6f133f8199438c84c1ac Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Sat, 7 Jul 2018 18:10:08 +0200
Subject: [PATCH] Added remap_cell_search.cc.
---
src/Makefile.am | 3 ++-
src/grid_cell_search.cc | 47 ++++++++++++++++++++++++++++++++++++++--
src/grid_cell_search.h | 7 ++++++
src/grid_point_search.cc | 8 +++----
src/remap.h | 1 +
src/remap_cell_search.cc | 43 ++++++++++++++++++++++++++++++++++++
src/remap_conserv.cc | 22 +++----------------
7 files changed, 105 insertions(+), 26 deletions(-)
create mode 100644 src/remap_cell_search.cc
diff --git a/src/Makefile.am b/src/Makefile.am
index 17e88ce14..aec2f6865 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -167,7 +167,8 @@ libcdo_la_SOURCES = \
remap_scrip_io.cc \
remap_search_reg2d.cc \
remap_search_latbins.cc \
- remap_search.cc \
+ remap_point_search.cc \
+ remap_cell_search.cc \
remap_store_link.cc \
remap_store_link.h \
remap_store_link_cnsrv.cc \
diff --git a/src/grid_cell_search.cc b/src/grid_cell_search.cc
index 3890193da..f1442e351 100644
--- a/src/grid_cell_search.cc
+++ b/src/grid_cell_search.cc
@@ -46,8 +46,6 @@ gridCellSearchCreateReg2d(size_t dims[2], const double *reg2d_corner_lon, const
size_t nx = dims[0];
size_t ny = dims[1];
- printf("nx, ny %zu %zu\n", nx, ny);
-
gcs->reg2d_corner_lon = (double *) Malloc(nx * sizeof(double));
gcs->reg2d_corner_lat = (double *) Malloc(ny * sizeof(double));
@@ -62,6 +60,8 @@ gridCellSearchCreate(size_t numCells, size_t numCellCorners, double *cellCornerL
{
GridCellSearch *gcs = (GridCellSearch *) Calloc(1, sizeof(GridCellSearch));
+ gcs->method = cellSearchMethod;
+
struct grid *yacSrcGrid = yac_scrip_grid_new(cellCornerLon, cellCornerLat, numCells, numCellCorners);
gcs->yacSearch = yac_sphere_part_search_new(yacSrcGrid);
gcs->yacSrcGrid = yacSrcGrid;
@@ -83,3 +83,46 @@ gridCellSearchDelete(GridCellSearch *gcs)
Free(gcs);
}
}
+
+static size_t
+doGridCellSearchYac(grid_search *yacSearch, bool isLonLatCell, size_t numCellCorners, grid_cell &gridCell, size_t *srchAddr)
+{
+ struct bounding_circle bndCircle;
+ double *xyz = gridCell.coordinates_xyz;
+ if (numCellCorners == 4 && isLonLatCell)
+ yac_get_cell_bounding_circle_reg_quad(xyz, xyz+3, xyz+6, xyz+9, &bndCircle);
+ else if (numCellCorners == 3)
+ yac_get_cell_bounding_circle_unstruct_triangle(xyz, xyz+3, xyz+6, &bndCircle);
+ else
+ yac_get_cell_bounding_circle(gridCell, &bndCircle);
+
+ struct dep_list resultList;
+ yac_do_bnd_circle_search(yacSearch, &bndCircle, 1, &resultList);
+
+ size_t numSrchCells = resultList.num_deps_per_element[0];
+ unsigned const *currNeighs = yac_get_dependencies_of_element(resultList, 0);
+ for (size_t i = 0; i < numSrchCells; ++i) srchAddr[i] = currNeighs[i];
+
+ yac_free_dep_list(&resultList);
+
+ return numSrchCells;
+}
+
+size_t
+doGridCellSearch(GridCellSearch *gcs, bool isLonLatCell, size_t numCellCorners, grid_cell &gridCell, size_t *srchAddr)
+{
+ if (gcs)
+ {
+ size_t nadds = 0;
+ // clang-format off
+ if (gcs->method == CellSearchMethod::spherepart)
+ nadds = doGridCellSearchYac((grid_search*)gcs->yacSearch, isLonLatCell, numCellCorners, gridCell, srchAddr);
+ else
+ cdoAbort("%s::method undefined!", __func__);
+ // clang-format on
+
+ return nadds;
+ }
+
+ return 0;
+}
diff --git a/src/grid_cell_search.h b/src/grid_cell_search.h
index f97f7d181..1a3b39592 100644
--- a/src/grid_cell_search.h
+++ b/src/grid_cell_search.h
@@ -18,6 +18,10 @@
#ifndef GRID_CELL_SEARCH_H
#define GRID_CELL_SEARCH_H
+extern "C" {
+#include "lib/yac/grid_cell.h"
+}
+
enum struct CellSearchMethod
{
spherepart,
@@ -26,6 +30,8 @@ enum struct CellSearchMethod
struct GridCellSearch
{
+ CellSearchMethod method;
+
bool is_reg2d;
size_t dims[2];
@@ -40,5 +46,6 @@ struct GridCellSearch
GridCellSearch *gridCellSearchCreateReg2d(size_t dims[2], const double *reg2d_corner_lon, const double *reg2d_corner_lat);
GridCellSearch *gridCellSearchCreate(size_t numCells, size_t numCellCorners, double *cellCornerLon, double *cellCornerLat);
void gridCellSearchDelete(GridCellSearch *gcs);
+size_t doGridCellSearch(GridCellSearch *gcs, bool isLonLatCell, size_t numCellCorners, grid_cell &gridCell, size_t *srchAddr);
#endif
diff --git a/src/grid_point_search.cc b/src/grid_point_search.cc
index f8ed42fe3..7887d70f6 100644
--- a/src/grid_point_search.cc
+++ b/src/grid_point_search.cc
@@ -616,10 +616,10 @@ gridPointSearchNearest(GridPointSearch *gps, double lon, double lat, size_t *add
double searchRadius = gps->searchRadius;
void *sc = gps->search_container;
// clang-format off
- if ( gps->method == PointSearchMethod::kdtree ) nadds = gps_nearest_kdtree(sc, lon, lat, searchRadius, addr, dist, gps);
- else if ( gps->method == PointSearchMethod::nanoflann ) nadds = gps_nearest_nanoflann(sc, lon, lat, searchRadius, addr, dist, gps);
- else if ( gps->method == PointSearchMethod::spherepart ) nadds = gps_nearest_spherepart(sc, lon, lat, searchRadius, addr, dist, gps);
- else if ( gps->method == PointSearchMethod::full ) nadds = gps_nearest_full(sc, lon, lat, searchRadius, addr, dist);
+ if (gps->method == PointSearchMethod::kdtree) nadds = gps_nearest_kdtree(sc, lon, lat, searchRadius, addr, dist, gps);
+ else if (gps->method == PointSearchMethod::nanoflann) nadds = gps_nearest_nanoflann(sc, lon, lat, searchRadius, addr, dist, gps);
+ else if (gps->method == PointSearchMethod::spherepart) nadds = gps_nearest_spherepart(sc, lon, lat, searchRadius, addr, dist, gps);
+ else if (gps->method == PointSearchMethod::full) nadds = gps_nearest_full(sc, lon, lat, searchRadius, addr, dist);
else cdoAbort("%s::method undefined!", __func__);
// clang-format on
diff --git a/src/remap.h b/src/remap.h
index 54a83b9dc..9ddfa7ab5 100644
--- a/src/remap.h
+++ b/src/remap.h
@@ -138,6 +138,7 @@ void remapSearchFree(RemapSearch &search);
void remapSearchPoints(RemapSearch &rsearch, double plon, double plat, knnWeightsType &knnWeights);
int remapSearchSquare(RemapSearch &rsearch, double plon, double plat, size_t *src_add, double *src_lats, double *src_lons);
+size_t remapSearchCells(RemapSearch &rsearch, bool isLonLatCell, size_t numCellCorners, grid_cell &gridCell, size_t *srchAddr);
void remapBilinearWeights(RemapSearch &rsearch, RemapVars &rv);
void remapBicubicWeights(RemapSearch &rsearch, RemapVars &rv);
diff --git a/src/remap_cell_search.cc b/src/remap_cell_search.cc
new file mode 100644
index 000000000..24bfd121a
--- /dev/null
+++ b/src/remap_cell_search.cc
@@ -0,0 +1,43 @@
+/*
+ This file is part of CDO. CDO is a collection of Operators to
+ manipulate and analyse Climate model Data.
+
+ Copyright (C) 2003-2018 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
+ See COPYING file for copying and redistribution conditions.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; version 2 of the License.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+*/
+
+#include "cdo_int.h"
+#include "remap.h"
+
+static void
+gridSearchCellsReg2d(GridCellSearch *gcs)
+{
+}
+
+static size_t
+gridSearchCells(GridCellSearch *gcs, bool isLonLatCell, size_t numCellCorners, grid_cell &gridCell, size_t *srchAddr)
+{
+ return doGridCellSearch(gcs, isLonLatCell, numCellCorners, gridCell, srchAddr);
+}
+
+size_t
+remapSearchCells(RemapSearch &rsearch, bool isLonLatCell, size_t numCellCorners, grid_cell &gridCell, size_t *srchAddr)
+{
+ int nadds = 0;
+
+ if (rsearch.srcGrid->remap_grid_type == REMAP_GRID_TYPE_REG2D)
+ gridSearchCellsReg2d(rsearch.gcs);
+ else
+ nadds = gridSearchCells(rsearch.gcs, isLonLatCell, numCellCorners, gridCell, srchAddr);
+
+ return nadds;
+}
diff --git a/src/remap_conserv.cc b/src/remap_conserv.cc
index bcedc0d24..7638c3221 100644
--- a/src/remap_conserv.cc
+++ b/src/remap_conserv.cc
@@ -753,26 +753,10 @@ remapConservWeights(RemapSearch &rsearch, RemapVars &rv)
// Get search cells
- if ( useCellsearch )
+ if (useCellsearch)
{
- struct bounding_circle bnd_circle;
- double *xyz = tgt_grid_cell[ompthID].coordinates_xyz;
- if (tgt_num_cell_corners == 4 && target_cell_type == LON_LAT_CELL)
- yac_get_cell_bounding_circle_reg_quad(xyz, xyz+3, xyz+6, xyz+9, &bnd_circle);
- else if (tgt_num_cell_corners == 3)
- yac_get_cell_bounding_circle_unstruct_triangle(xyz, xyz+3, xyz+6, &bnd_circle);
- else
- yac_get_cell_bounding_circle(tgt_grid_cell[ompthID], &bnd_circle);
-
- struct dep_list result_list;
- yac_do_bnd_circle_search((struct grid_search *)rsearch.gcs->yacSearch, &bnd_circle, 1, &result_list);
-
- num_srch_cells = result_list.num_deps_per_element[0];
- unsigned const *curr_neighs = yac_get_dependencies_of_element(result_list,0);
- for (size_t i = 0; i < num_srch_cells; ++i)
- srch_add[ompthID][i] = curr_neighs[i];
-
- yac_free_dep_list(&result_list);
+ num_srch_cells = remapSearchCells(rsearch, target_cell_type == LON_LAT_CELL, tgt_num_cell_corners,
+ tgt_grid_cell[ompthID], srch_add[ompthID].data());
}
else if (src_remap_grid_type == REMAP_GRID_TYPE_REG2D)
{
--
GitLab