Commit 3b3282c7 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

gridsearch: added cdo default search radius

parent 49d45495
......@@ -31,14 +31,72 @@
#include "grid_search.h"
void grid_search_nbr(struct gridsearch *gs, int num_neighbors, int *restrict nbr_add, double *restrict nbr_dist, double plon, double plat);
double nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_mask, int *restrict nbr_mask, const int *restrict nbr_add, double *restrict nbr_dist);
unsigned nbr_normalize_weights(unsigned num_neighbors, double dist_tot, const int *restrict nbr_mask, int *restrict nbr_add, double *restrict nbr_dist);
double intlin(double x, double y1, double x1, double y2, double x2);
double smooth_nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_mask, int *restrict nbr_mask, const int *restrict nbr_add, double *restrict nbr_dist, double search_radius, double weight0, double weightInf)
{
// Compute weights based on inverse distance if mask is false, eliminate those points
double dist_tot = 0.; // sum of neighbor distances (for normalizing)
if ( src_grid_mask )
{
for ( unsigned n = 0; n < num_neighbors; ++n )
{
nbr_mask[n] = FALSE;
if ( nbr_add[n] >= 0 )
if ( src_grid_mask[nbr_add[n]] )
{
nbr_dist[n] = 1./nbr_dist[n];
dist_tot += nbr_dist[n];
nbr_mask[n] = TRUE;
}
}
}
else
{
for ( unsigned n = 0; n < num_neighbors; ++n )
{
nbr_mask[n] = FALSE;
if ( nbr_add[n] >= 0 )
{
nbr_dist[n] = intlin(nbr_dist[n], weight0, 0, weightInf, search_radius);
dist_tot += nbr_dist[n];
nbr_mask[n] = TRUE;
}
}
}
return dist_tot;
}
unsigned smooth_nbr_normalize_weights(unsigned num_neighbors, double dist_tot, const int *restrict nbr_mask, int *restrict nbr_add, double *restrict nbr_dist)
{
// Normalize weights and store the link
unsigned nadds = 0;
for ( unsigned n = 0; n < num_neighbors; ++n )
{
if ( nbr_mask[n] )
{
nbr_dist[nadds] = nbr_dist[n]/dist_tot;
nbr_add[nadds] = nbr_add[n];
nadds++;
}
}
return nadds;
}
static
void smooth(int gridID, double missval, const double *restrict array1, double *restrict array2, int *nmiss)
void smooth(int gridID, double missval, const double *restrict array1, double *restrict array2, int *nmiss,
int num_neighbors, double search_radius, double weight0, double weightInf)
{
*nmiss = 0;
int num_neighbors = 9;
int gridID0 = gridID;
unsigned gridsize = gridInqSize(gridID);
......@@ -73,6 +131,8 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
gs = gridsearch_create_nn(gridsize, xvals, yvals);
else
gs = gridsearch_create(gridsize, xvals, yvals);
gs->search_radius = search_radius;
finish = clock();
......@@ -84,30 +144,35 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
double findex = 0;
/*
#pragma omp parallel for default(none) shared(findex, array1, array2, xvals, yvals, gs, gridsize, num_neighbors) \
private(nbr_mask, nbr_add, nbr_dist)
*/
for ( unsigned i = 0; i < gridsize; ++i )
{
/*
#if defined(_OPENMP)
#include "pragma_omp_atomic_update.h"
#endif
*/
findex++;
if ( cdo_omp_get_thread_num() == 0 ) progressStatus(0, 1, findex/gridsize);
grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, xvals[i], yvals[i]);
/* Compute weights based on inverse distance if mask is false, eliminate those points */
double dist_tot = nbr_compute_weights(num_neighbors, NULL, nbr_mask, nbr_add, nbr_dist);
double dist_tot = smooth_nbr_compute_weights(num_neighbors, NULL, nbr_mask, nbr_add, nbr_dist,
search_radius, weight0, weightInf);
/* Normalize weights and store the link */
unsigned nadds = nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
unsigned nadds = smooth_nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
if ( nadds )
{
/*
printf("n %u %d nadds %u dis %g\n", i, nbr_add[0], nadds, nbr_dist[0]);
for ( unsigned n = 0; n < nadds; ++n )
printf(" n %u add %d dis %g\n", n, nbr_add[n], nbr_dist[n]);
*/
double result = 0;
for ( unsigned n = 0; n < nadds; ++n ) result += array1[nbr_add[n]]*nbr_dist[n];
array2[i] = result;
......@@ -126,12 +191,6 @@ void smooth(int gridID, double missval, const double *restrict array1, double *r
Free(yvals);
}
static inline
void smooth_sum(short m, double sfac, double val, double *avg, double *divavg)
{
if ( m ) { *avg += sfac*val; *divavg += sfac; }
}
static inline
void smooth9_sum(size_t ij, short *mask, double sfac, const double *restrict array, double *avg, double *divavg)
{
......@@ -251,6 +310,10 @@ void *Smooth(void *argument)
int varID, levelID;
int nmiss;
int gridtype;
int max_points = 5;
double search_radius = 60, weight0 = 0.25, weightInf = 0.;
search_radius *= DEG2RAD;
cdoInitialize(argument);
......@@ -291,7 +354,7 @@ void *Smooth(void *argument)
size_t gridsize = vlistGridsizeMax(vlistID1);
double *array1 = (double*) Malloc(gridsize*sizeof(double));
double *array2 = (double*) Malloc(gridsize *sizeof(double));
double *array2 = (double*) Malloc(gridsize*sizeof(double));
int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
......@@ -314,7 +377,7 @@ void *Smooth(void *argument)
gridID = vlistInqVarGrid(vlistID1, varID);
if ( operatorID == SMOOTH )
smooth(gridID, missval, array1, array2, &nmiss);
smooth(gridID, missval, array1, array2, &nmiss, max_points, search_radius, weight0, weightInf);
else if ( operatorID == SMOOTH9 )
smooth9(gridID, missval, array1, array2, &nmiss);
......
......@@ -79,11 +79,11 @@ enum T_EIGEN_MODE {JACOBI, DANIELSON_LANCZOS};
#ifndef M_LN10
#define M_LN10 2.30258509299404568402 /* log_e 10 */
#define M_LN10 2.30258509299404568401799145468436421 /* log_e 10 */
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
#define M_PI 3.14159265358979323846264338327950288 /* pi */
#endif
......
......@@ -6,7 +6,7 @@
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
#define M_PI 3.14159265358979323846264338327950288 /* pi */
#endif
......
......@@ -10,19 +10,31 @@
#include "cdo_int.h"
#include "dmemory.h"
#include "util.h"
#include "grid.h"
#include "grid_search.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288 /* pi */
#endif
#define PI M_PI
#define PI2 (2.0*PI)
static int gridsearch_method_nn = GS_KDTREE;
static
double cdo_default_search_radius(void)
{
extern double gridsearch_radius;
double search_radius = gridsearch_radius;
if ( search_radius < 0. ) search_radius = 0.;
if ( search_radius > 180. ) search_radius = 180.;
search_radius = search_radius*DEG2RAD;
return search_radius;
}
static inline void LLtoXYZ_f(double lon, double lat, float *restrict xyz)
{
double cos_lat = cos(lat);
......@@ -105,6 +117,8 @@ struct gridsearch *gridsearch_create_reg2d(unsigned lcyclic, unsigned nx, unsign
gs->coslat = coslat;
gs->sinlat = sinlat;
gs->search_radius = cdo_default_search_radius();
return gs;
}
......@@ -243,6 +257,8 @@ struct gridsearch *gridsearch_create(unsigned n, const double *restrict lons, co
gs->kdt = gs_create_kdtree(n, lons, lats);
gs->search_radius = cdo_default_search_radius();
return gs;
}
......@@ -259,6 +275,8 @@ struct gridsearch *gridsearch_create_nn(unsigned n, const double *restrict lons,
else if ( gs->method_nn == GS_NEARPT3 ) gs->near = gs_create_nearpt3(n, lons, lats);
else if ( gs->method_nn == GS_FULL ) gs->full = gs_create_full(n, lons, lats);
gs->search_radius = cdo_default_search_radius();
return gs;
}
......
......@@ -36,6 +36,8 @@ struct gridsearch {
struct gsFull *full;
double search_radius;
// reg2d search
double *reg2d_center_lon, *reg2d_center_lat;
double *coslat, *sinlat; // cosine, sine of grid lats (for distance)
......
......@@ -74,8 +74,8 @@ double nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_
if ( nbr_add[n] >= 0 )
if ( src_grid_mask[nbr_add[n]] )
{
nbr_dist[n] = ONE/nbr_dist[n];
dist_tot = dist_tot + nbr_dist[n];
nbr_dist[n] = 1./nbr_dist[n];
dist_tot += nbr_dist[n];
nbr_mask[n] = TRUE;
}
}
......@@ -87,8 +87,8 @@ double nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_
nbr_mask[n] = FALSE;
if ( nbr_add[n] >= 0 )
{
nbr_dist[n] = ONE/nbr_dist[n];
dist_tot = dist_tot + nbr_dist[n];
nbr_dist[n] = 1./nbr_dist[n];
dist_tot += nbr_dist[n];
nbr_mask[n] = TRUE;
}
}
......@@ -117,20 +117,6 @@ unsigned nbr_normalize_weights(unsigned num_neighbors, double dist_tot, const in
return nadds;
}
static
double get_search_radius(void)
{
extern double gridsearch_radius;
double search_radius = gridsearch_radius;
if ( search_radius < 0. ) search_radius = 0.;
if ( search_radius > 180. ) search_radius = 180.;
search_radius = search_radius*DEG2RAD;
return search_radius;
}
// This routine finds the closest num_neighbor points to a search point and computes a distance to each of the neighbors.
......@@ -159,7 +145,7 @@ void grid_search_nbr_reg2d(struct gridsearch *gs, int num_neighbors, remapgrid_t
int *psrc_add = src_add;
int num_add = 0;
double distance; // Angular distance
double cos_search_radius = cos(get_search_radius());
double cos_search_radius = cos(gs->search_radius);
double coslat_dst = cos(plat); // cos(lat) of the search point
double coslon_dst = cos(plon); // cos(lon) of the search point
double sinlat_dst = sin(plat); // sin(lat) of the search point
......@@ -297,7 +283,7 @@ void grid_search_nbr(struct gridsearch *gs, int num_neighbors, int *restrict nbr
double plon, ! longitude of the search point
*/
double search_radius = get_search_radius();
double search_radius = gs->search_radius;
// Initialize distance and address arrays
for ( int n = 0; n < num_neighbors; ++n ) nbr_add[n] = -1;
......
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