Skip to content
GitLab
Commit 63fec1bc authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Removed operator intgridcon.

parent 9af4fe20
Hide whitespace changes
Inline Side-by-side
src/Intgrid.c
View file @ 63fec1bc
......@@ -291,12 +291,10 @@ void *Intgrid(void *argument)
int nmiss;
int xinc = 0, yinc = 0;
double missval;
double slon, slat;
cdoInitialize(argument);
int INTGRIDBIL = cdoOperatorAdd("intgridbil", 0, 0, NULL);
int INTGRIDCON = cdoOperatorAdd("intgridcon", 0, 0, NULL);
int INTPOINT = cdoOperatorAdd("intpoint", 0, 0, NULL);
int INTERPOLATE = cdoOperatorAdd("interpolate", 0, 0, NULL);
int BOXAVG = cdoOperatorAdd("boxavg", 0, 0, NULL);
......@@ -306,7 +304,7 @@ void *Intgrid(void *argument)
int streamID1 = streamOpenRead(cdoStreamName(0));
if ( operatorID == INTGRIDBIL || operatorID == INTGRIDCON || operatorID == INTERPOLATE )
if ( operatorID == INTGRIDBIL || operatorID == INTERPOLATE )
{
operatorInputArg("grid description file or name");
gridID2 = cdoDefineGrid(operatorArgv()[0]);
......@@ -315,8 +313,8 @@ void *Intgrid(void *argument)
{
operatorInputArg("longitude and latitude");
operatorCheckArgc(2);
slon = parameter2double(operatorArgv()[0]);
slat = parameter2double(operatorArgv()[1]);
double slon = parameter2double(operatorArgv()[0]);
double slat = parameter2double(operatorArgv()[1]);
gridID2 = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID2, 1);
gridDefYsize(gridID2, 1);
......@@ -410,8 +408,6 @@ void *Intgrid(void *argument)
if ( operatorID == INTGRIDBIL || operatorID == INTPOINT )
intgridbil(&field1, &field2);
if ( operatorID == INTGRIDCON )
intgridcon(&field1, &field2);
else if ( operatorID == INTERPOLATE )
interpolate(&field1, &field2);
else if ( operatorID == BOXAVG )
......
src/interpol.c
View file @ 63fec1bc
......@@ -202,26 +202,6 @@ double intlinarr2p(long nxm, long nym, double **fieldm, const double *xm, const
return value;
}
static
void rect_find_ij_weights(double plon, double plat, long ii, long jj, const double *restrict xm, const double *restrict ym, double *ig, double *jg)
{
/*
wgts[0] = (plon-xm[ii]) * (plat-ym[jj]) / ((xm[ii-1]-xm[ii]) * (ym[jj-1]-ym[jj]));
wgts[1] = (plon-xm[ii-1]) * (plat-ym[jj]) / ((xm[ii]-xm[ii-1]) * (ym[jj-1]-ym[jj]));
wgts[2] = (plon-xm[ii-1]) * (plat-ym[jj-1]) / ((xm[ii]-xm[ii-1]) * (ym[jj]-ym[jj-1]));
wgts[3] = (plon-xm[ii]) * (plat-ym[jj-1]) / ((xm[ii-1]-xm[ii]) * (ym[jj]-ym[jj-1]));
*/
double iw, jw;
double wgts0 = (plon-xm[ii]) * (plat-ym[jj]) / ((xm[ii-1]-xm[ii]) * (ym[jj-1]-ym[jj]));
double wgts1 = (plon-xm[ii-1]) * (plat-ym[jj]) / ((xm[ii]-xm[ii-1]) * (ym[jj-1]-ym[jj]));
iw = 1./(wgts0/wgts1 + 1.);
jw = 1. - wgts1/iw;
*ig = iw;
*jg = jw;
}
static
void intlinarr2(double missval, int lon_is_circular,
long nxm, long nym, const double *restrict fieldm, const double *restrict xm, const double *restrict ym,
......@@ -319,342 +299,6 @@ void intlinarr2(double missval, int lon_is_circular,
if ( grid1_mask ) Free(grid1_mask);
}
static
void restrict_boundbox(const double *restrict grid_bound_box, double *restrict bound_box)
{
if ( bound_box[0] < grid_bound_box[0] && bound_box[1] > grid_bound_box[0] ) bound_box[0] = grid_bound_box[0];
if ( bound_box[1] > grid_bound_box[1] && bound_box[0] < grid_bound_box[1] ) bound_box[1] = grid_bound_box[1];
// if ( bound_box[2] < grid_bound_box[2] && bound_box[3] > grid_bound_box[2] ) bound_box[2] = grid_bound_box[2];
// if ( bound_box[3] > grid_bound_box[3] && bound_box[2] < grid_bound_box[3] ) bound_box[3] = grid_bound_box[3];
}
static
void boundbox_from_corners(long ic, long nc, const double *restrict corner_lon,
const double *restrict corner_lat, double *restrict bound_box)
{
long inc, j;
double clon, clat;
inc = ic*nc;
clat = corner_lat[inc];
clon = corner_lon[inc];
bound_box[0] = clat;
bound_box[1] = clat;
bound_box[2] = clon;
bound_box[3] = clon;
for ( j = 1; j < nc; ++j )
{
clat = corner_lat[inc+j];
clon = corner_lon[inc+j];
if ( clat < bound_box[0] ) bound_box[0] = clat;
if ( clat > bound_box[1] ) bound_box[1] = clat;
if ( clon < bound_box[2] ) bound_box[2] = clon;
if ( clon > bound_box[3] ) bound_box[3] = clon;
}
}
#if defined(HAVE_LIBYAC)
#include "clipping.h"
#include "area.h"
#endif
static
void intconarr2(double missval, int lon_is_circular,
long nxm, long nym, const double *restrict fieldm, const double *restrict xm, const double *restrict ym,
long nc2, long gridsize2, double *field, const double *restrict x, const double *restrict y)
{
long ndeps;
int *deps;
long ii = -1, jj = -1;
long gridsize1;
long nlon1 = nxm;
long nx, ny;
double findex = 0;
int *grid1_mask = NULL;
nx = nxm - 1;
ny = nym - 1;
printf(" nx, ny %ld %ld\n", nx, ny);
//if ( lon_is_circular ) nlon1--;
gridsize1 = nx*ny;
deps = (int*) Malloc(gridsize1*sizeof(int));
grid1_mask = (int*) Calloc(1, gridsize1*sizeof(int));
for ( jj = 0; jj < ny; ++jj )
for ( ii = 0; ii < nx; ++ii )
{
if ( !DBL_IS_EQUAL(fieldm[jj*nx+ii], missval) ) grid1_mask[jj*nx+ii] = 1;
}
double grid1_bound_box[4];
grid1_bound_box[0] = ym[0];
grid1_bound_box[1] = ym[ny];
if ( ym[0] > ym[ny] )
{
grid1_bound_box[0] = ym[ny];
grid1_bound_box[1] = ym[0];
}
grid1_bound_box[2] = xm[0];
grid1_bound_box[3] = xm[nx];
progressInit();
#if defined(HAVE_LIBYAC)
enum edge_type quad_type[] = {GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE}; // not used !
// enum edge_type quad_type[] = {LON_CIRCLE, LON_CIRCLE, LON_CIRCLE, LON_CIRCLE};
int n;
double weight_sum;
double *weight;
weight = (double*) Malloc(gridsize1*sizeof(double));
double tgt_area;
double *area;
area = (double*) Malloc(gridsize1*sizeof(double));
struct grid_cell *SourceCell;
SourceCell = (struct grid_cell*) Malloc(gridsize1 * sizeof(struct grid_cell));
for ( int n = 0; n < gridsize1; n++ ) {
SourceCell[n].num_corners = 4;
SourceCell[n].edge_type = quad_type;
SourceCell[n].coordinates_x = (double*) Malloc(4 * sizeof(double));
SourceCell[n].coordinates_y = (double*) Malloc(4 * sizeof(double));
}
struct grid_cell TargetCell;
TargetCell.num_corners = nc2;
TargetCell.edge_type = quad_type;
TargetCell.coordinates_x = (double*) Malloc(nc2 * sizeof(double));
TargetCell.coordinates_y = (double*) Malloc(nc2 * sizeof(double));
unsigned const * curr_deps;
#endif
/*
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, field, fieldm, x, y, xm, ym, nxm, nym, gridsize2, missval, findex, nlon1, lon_is_circular, grid1_mask, nc2) \
private(jj, ii)
#endif
*/
for ( int i = 0; i < gridsize2; ++i )
{
int src_add[4]; /* address for the four source points */
long n;
long iix;
int lfound;
int lprogress = 1;
ndeps = 0;
/*
#if defined(_OPENMP)
if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;
#endif
*/
field[i] = missval;
// printf("%ld lonb: %g %g %g %g latb: %g %g %g %g\n", i+1, x[i*nc2]*RAD2DEG, x[i*nc2+1]*RAD2DEG, x[i*nc2+2]*RAD2DEG, x[i*nc2+3]*RAD2DEG, y[i*nc2]*RAD2DEG, y[i*nc2+1]*RAD2DEG, y[i*nc2+2]*RAD2DEG, y[i*nc2+3]*RAD2DEG);
double bound_lon1, bound_lon2;
double bound_box[4];
boundbox_from_corners(i, nc2, x, y, bound_box);
restrict_boundbox(grid1_bound_box, bound_box);
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
//printf("bound_box %ld lon: %g %g lat: %g %g\n", i, bound_box[2]*RAD2DEG, bound_box[3]*RAD2DEG, bound_box[0]*RAD2DEG, bound_box[1]*RAD2DEG);
/*
#if defined(_OPENMP)
#include "pragma_omp_atomic_update.h"
#endif
*/
findex++;
if ( lprogress ) progressStatus(0, 1, findex/gridsize2);
// lfound = rect_grid_search(&ii, &jj, x[i], y[i], nxm, nym, xm, ym);
// for ( int k = 0; k < nxm; ++k ) printf("x: %d %g\n", k+1, xm[k]);
// for ( int k = 0; k < nym; ++k ) printf("y: %d %g\n", k+1, ym[k]*RAD2DEG);
long imin = nxm, imax = -1, jmin = nym, jmax = -1;
long im, jm;
lfound = rect_grid_search2(&jmin, &jmax, bound_box[0], bound_box[1], nym, ym);
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
if ( bound_lon1 <= grid1_bound_box[3] && bound_lon2 >= grid1_bound_box[2] )
{
//printf("b1 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < grid1_bound_box[2] && bound_lon2 > grid1_bound_box[2] ) bound_lon1 = grid1_bound_box[2];
if ( bound_lon2 > grid1_bound_box[3] && bound_lon1 < grid1_bound_box[3] ) bound_lon2 = grid1_bound_box[3];
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
deps[ndeps++] = jm*nx + im;
}
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
if ( bound_lon1 <= grid1_bound_box[2] && bound_lon2 > grid1_bound_box[2] )
{
bound_lon1 += 2*M_PI;
bound_lon2 += 2*M_PI;
//printf("b2 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < grid1_bound_box[2] && bound_lon2 > grid1_bound_box[2] ) bound_lon1 = grid1_bound_box[2];
if ( bound_lon2 > grid1_bound_box[3] && bound_lon1 < grid1_bound_box[3] ) bound_lon2 = grid1_bound_box[3];
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
deps[ndeps++] = jm*nx + im;
}
bound_lon1 = bound_box[2];
bound_lon2 = bound_box[3];
if ( bound_lon1 < grid1_bound_box[3] && bound_lon2 >= grid1_bound_box[3] )
{
bound_lon1 -= 2*M_PI;
bound_lon2 -= 2*M_PI;
//printf("b3 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < grid1_bound_box[2] && bound_lon2 > grid1_bound_box[2] ) bound_lon1 = grid1_bound_box[2];
if ( bound_lon2 > grid1_bound_box[3] && bound_lon1 < grid1_bound_box[3] ) bound_lon2 = grid1_bound_box[3];
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxm, xm);
//printf("imin %ld imax %ld jmin %ld jmax %ld\n", imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
deps[ndeps++] = jm*nx + im;
}
//for ( long id = 0; id < ndeps; ++id )
// printf("dep %ld %d\n", id+1, deps[id]);
/*
if ( bound_lon1 < grid1_bound_box[2] && bound_box[3] >= grid1_bound_box[2] )
lfound = rect_grid_search2(&imin, &imax, bound_box[2], bound_box[3], nxm, xm);
*/
#if defined(HAVE_LIBYAC)
int index2 = i;
/*
int ilat2 = index2/nlonOut;
int ilon2 = index2 - ilat2*nlonOut;
*/
double addtest = 0;
//if ( i == 296 ) addtest = 360*DEG2RAD;
for ( int ic = 0; ic < nc2; ++ic )
{
TargetCell.coordinates_x[ic] = addtest+x[index2*nc2+ic];
TargetCell.coordinates_y[ic] = y[index2*nc2+ic];
}
if ( cdoVerbose )
{
printf("target: ");
for ( int n = 0; n < nc2; ++n )
printf(" %g %g", TargetCell.coordinates_x[n]/DEG2RAD, TargetCell.coordinates_y[n]/DEG2RAD);
printf("\n");
}
/*
if ( cdoVerbose )
printf("num_deps_per_element %d %d\n", i, ndeps);
*/
int num_deps = ndeps;
int nSourceCells = num_deps;
for ( int k = 0; k < num_deps; ++k )
{
int index1 = deps[k];
int ilat1 = index1/nx;
int ilon1 = index1 - ilat1*nx;
/*
if ( cdoVerbose )
printf(" dep: %d %d %d %d %d %d\n", k, nlonOut, nlatOut, index1, ilon1, ilat1);
*/
if ( ym[ilat1] < ym[ilat1+1] )
{
SourceCell[k].coordinates_x[0] = addtest+xm[ilon1];
SourceCell[k].coordinates_y[0] = ym[ilat1];
SourceCell[k].coordinates_x[1] = addtest+xm[ilon1+1];
SourceCell[k].coordinates_y[1] = ym[ilat1];
SourceCell[k].coordinates_x[2] = addtest+xm[ilon1+1];
SourceCell[k].coordinates_y[2] = ym[ilat1+1];
SourceCell[k].coordinates_x[3] = addtest+xm[ilon1];
SourceCell[k].coordinates_y[3] = ym[ilat1+1];
}
else
{
SourceCell[k].coordinates_x[0] = addtest+xm[ilon1];
SourceCell[k].coordinates_y[0] = ym[ilat1+1];
SourceCell[k].coordinates_x[1] = addtest+xm[ilon1+1];
SourceCell[k].coordinates_y[1] = ym[ilat1+1];
SourceCell[k].coordinates_x[2] = addtest+xm[ilon1+1];
SourceCell[k].coordinates_y[2] = ym[ilat1];
SourceCell[k].coordinates_x[3] = addtest+xm[ilon1];
SourceCell[k].coordinates_y[3] = ym[ilat1];
}
if ( cdoVerbose )
{
printf("source: %d %d", num_deps, k);
for ( int n = 0; n < 4; ++n )
printf(" %g %g", SourceCell[k].coordinates_x[n]/DEG2RAD, SourceCell[k].coordinates_y[n]/DEG2RAD);
printf("\n");
}
}
compute_overlap_areas ( nSourceCells, SourceCell, TargetCell, area);
tgt_area = huiliers_area(TargetCell);
// tgt_area = cell_area(TargetCell);
for (n = 0; n < nSourceCells; ++n)
weight[n] = area[n] / tgt_area;
correct_weights ( nSourceCells, weight );
field[i] = missval;
if ( num_deps ) field[i] = 0;
for ( int k = 0; k < num_deps; ++k )
{
int index1 = deps[k];
/*
int ilat1 = index1/nx;
int ilon1 = index1 - ilat1*nx;
long add1, add2;
add1 = index1;
add2 = index2;
yar_store_link_cnsrv(&remap.vars, add1, add2, weight[k]);
*/
if ( weight[k] > 0 )
{
if ( cdoVerbose )
printf("tgt_cell_add %ld, src_cell_add %ld, weight[n] %g, tgt_area %g\n", i, index1, weight[k], tgt_area);
//printf("tgt_cell_add %ld, n %ld, src_cell_add %ld, weight[n] %g, tgt_area %g\n", i, k, index1, weight[k], tgt_area);
field[i] += fieldm[index1] * weight[k];
}
/*
if ( cdoVerbose )
printf(" result dep: %ld %d %d %g\n", i, k, index1, weight[k]);
*/
}
#endif
}
#if defined(HAVE_LIBYAC)
Free(weight);
Free(area);
#endif
if ( findex < gridsize2 ) progressStatus(0, 1, 1);
if ( deps ) Free(deps);
if ( grid1_mask ) Free(grid1_mask);
}
double intlin(double x, double y1, double x1, double y2, double x2)
{
......@@ -855,133 +499,6 @@ void intgridbil(field_t *field1, field_t *field2)
Free(array1_2D);
}
void intgridcon(field_t *field1, field_t *field2)
{
int nlon1, nlat1;
int nlon1b, nlat1b;
int nlon2, nlat2;
int ilat;
int gridID1, gridID2;
int i, nmiss;
int lon_is_circular;
double *lon1, *lat1;
double *lon1bounds, *lat1bounds;
double **field;
double *array = NULL;
double *array1, *array2;
double missval;
char units[CDI_MAX_NAME];
/* static int index = 0; */
gridID1 = field1->grid;
gridID2 = field2->grid;
array1 = field1->ptr;
array2 = field2->ptr;
missval = field1->missval;
if ( ! (gridInqXvals(gridID1, NULL) && gridInqYvals(gridID1, NULL)) )
cdoAbort("Source grid has no values");
lon_is_circular = gridIsCircular(gridID1);
nlon1 = gridInqXsize(gridID1);
nlat1 = gridInqYsize(gridID1);
lon1 = (double*) Malloc(nlon1*sizeof(double));
lat1 = (double*) Malloc(nlat1*sizeof(double));
gridInqXvals(gridID1, lon1);
gridInqYvals(gridID1, lat1);
gridInqXunits(gridID1, units);
grid_to_radian(units, nlon1, lon1, "grid1 center lon");
grid_to_radian(units, nlat1, lat1, "grid1 center lat");
nlon1b = nlon1 + 1;
nlat1b = nlat1 + 1;
lon1bounds = (double*) Malloc(nlon1b*sizeof(double));
lat1bounds = (double*) Malloc(nlat1b*sizeof(double));
grid_gen_corners(nlon1, lon1, lon1bounds);
grid_gen_corners(nlat1, lat1, lat1bounds);
if ( cdoVerbose )
{
for ( int i = 0; i < nlat1; ++i )
printf("lat1 %d %g\n", i+1, lat1[i]*RAD2DEG);
printf("lat1bounds: %g %g %g ... %g %g %g\n", lat1bounds[0]*RAD2DEG, lat1bounds[1]*RAD2DEG, lat1bounds[2]*RAD2DEG, lat1bounds[nlat1b-3]*RAD2DEG, lat1bounds[nlat1b-2]*RAD2DEG, lat1bounds[nlat1b-1]*RAD2DEG);
printf("lon1bounds: %g %g %g ... %g %g %g\n", lon1bounds[0]*RAD2DEG, lon1bounds[1]*RAD2DEG, lon1bounds[2]*RAD2DEG, lon1bounds[nlon1b-3]*RAD2DEG, lon1bounds[nlon1b-2]*RAD2DEG, lon1bounds[nlon1b-1]*RAD2DEG);
}
nlon2 = gridInqXsize(gridID2);
nlat2 = gridInqYsize(gridID2);
int gridsize2;
if ( gridInqType(gridID2) == GRID_GME ) gridID2 = gridToUnstructured(gridID2, 0);
if ( gridInqType(gridID2) != GRID_UNSTRUCTURED && gridInqType(gridID2) != GRID_CURVILINEAR )
gridID2 = gridToCurvilinear(gridID2, 0);
if ( ! (gridInqXvals(gridID2, NULL) && gridInqYvals(gridID2, NULL)) )
cdoAbort("Target grid has no values");
gridsize2 = gridInqSize(gridID2);
int nc2;
if ( gridInqType(gridID2) == GRID_UNSTRUCTURED )
nc2 = gridInqNvertex(gridID2);
else
nc2 = 4;
double *grid2_corner_lon = NULL, *grid2_corner_lat = NULL;
if ( gridInqYbounds(gridID2, NULL) && gridInqXbounds(gridID2, NULL) )
{
grid2_corner_lon = (double*) Malloc(nc2*gridsize2*sizeof(double));
grid2_corner_lat = (double*) Malloc(nc2*gridsize2*sizeof(double));
gridInqXbounds(gridID2, grid2_corner_lon);
gridInqYbounds(gridID2, grid2_corner_lat);
}
else
{
cdoAbort("grid2 corner missing!");
}
gridInqXunits(gridID2, units);
grid_to_radian(units, nc2*gridsize2, grid2_corner_lon, "grid2 corner lon");
grid_to_radian(units, nc2*gridsize2, grid2_corner_lat, "grid2 corner lat");
/*
for ( int i = 0; i < gridsize2; ++i )
{
if ( lon2[i] < lon1[0] ) lon2[i] += 2*M_PI;
if ( lon2[i] > lon1[nlon1-1] ) lon2[i] -= 2*M_PI;
}
*/
for ( i = 0; i < gridsize2; ++i ) array2[i] = missval;
intconarr2(missval, lon_is_circular,
nlon1b, nlat1b, array1, lon1bounds, lat1bounds,
nc2, gridsize2, array2, grid2_corner_lon, grid2_corner_lat);
nmiss = 0;
for ( i = 0; i < gridsize2; ++i )
if ( DBL_IS_EQUAL(array2[i], missval) ) nmiss++;
field2->nmiss = nmiss;
if (grid2_corner_lon) Free(grid2_corner_lon);
if (grid2_corner_lat) Free(grid2_corner_lat);
if (array) Free(array);
Free(lon1);
Free(lat1);
}
/* source code from pingo */
void interpolate(field_t *field1, field_t *field2)
{
......
src/modules.c
View file @ 63fec1bc
......@@ -342,7 +342,7 @@ void *Maggraph(void *argument);
#define ImportobsOperators {"import_obs"}
#define InfoOperators {"info", "infop", "infon", "infoc", "map"}
#define InputOperators {"input", "inputsrv", "inputext"}
#define IntgridOperators {"intgridbil", "intgridcon", "intpoint", "interpolate", "boxavg", "thinout"}
#define IntgridOperators {"intgridbil", "intpoint", "interpolate", "boxavg", "thinout"}
#define IntgridtrajOperators {"intgridtraj"}
#define IntlevelOperators {"intlevel", "intlevelx"}
#define Intlevel3dOperators {"intlevel3d", "intlevelx3d"}
......
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