Commit fdfff5ae authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapconclip update

parent aa2b879b
......@@ -2,7 +2,7 @@
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2013 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
Copyright (C) 2003-2014 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
......@@ -233,7 +233,7 @@ void usage(void)
operatorPrintAll();
fprintf(stderr, "\n");
fprintf(stderr, " CDO version %s, Copyright (C) 2003-2013 Uwe Schulzweida\n", VERSION);
fprintf(stderr, " CDO version %s, Copyright (C) 2003-2014 Uwe Schulzweida\n", VERSION);
// fprintf(stderr, " Available from <http://code.zmaw.de/projects/cdo>\n");
fprintf(stderr, " This is free software and comes with ABSOLUTELY NO WARRANTY\n");
fprintf(stderr, " Report bugs to <http://code.zmaw.de/projects/cdo>\n");
......
......@@ -4321,11 +4321,11 @@ long get_srch_cells_reg2d(const int *restrict src_grid_dims,
bound_lon2 = tgt_cell_bound_box[3];
if ( bound_lon1 <= src_lon_max && bound_lon2 >= src_lon_min )
{
printf(" b1 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
//printf(" b1 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxp1, src_corner_lon);
printf(" %g %g imin %ld imax %ld jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin], RAD2DEG*src_corner_lon[imax+1], imin, imax, jmin, jmax);
//printf(" %g %g imin %ld imax %ld jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin], RAD2DEG*src_corner_lon[imax+1], imin, imax, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
srch_add[num_srch_cells++] = jm*nx + im;
......@@ -4335,15 +4335,18 @@ long get_srch_cells_reg2d(const int *restrict src_grid_dims,
bound_lon2 = tgt_cell_bound_box[3];
if ( bound_lon1 <= src_lon_min && bound_lon2 >= src_lon_min )
{
long imin2 = nxp1, imax2 = -1;
bound_lon1 += 2*M_PI;
bound_lon2 += 2*M_PI;
printf(" b2 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
//printf(" b2 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxp1, src_corner_lon);
printf(" %g %g imin %ld imax %ld jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin], RAD2DEG*src_corner_lon[imax+1], imin, imax, jmin, jmax);
lfound = rect_grid_search2(&imin2, &imax2, bound_lon1, bound_lon2, nxp1, src_corner_lon);
if ( imin2 != -1 && imin2 == imax ) imin2 += 1;
if ( imax2 != -1 && imax2 == imax ) imax2 += 1;
//printf(" %g %g imin %ld imax %ld jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin2], RAD2DEG*src_corner_lon[imax2+1], imin2, imax2, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
for ( im = imin2; im <= imax2; ++im )
srch_add[num_srch_cells++] = jm*nx + im;
}
......@@ -4351,15 +4354,18 @@ long get_srch_cells_reg2d(const int *restrict src_grid_dims,
bound_lon2 = tgt_cell_bound_box[3];
if ( bound_lon1 <= src_lon_max && bound_lon2 >= src_lon_max )
{
long imin3 = nxp1, imax3 = -1;
bound_lon1 -= 2*M_PI;
bound_lon2 -= 2*M_PI;
printf(" b3 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
//printf(" b3 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxp1, src_corner_lon);
printf(" %g %g imin %ld imax %ld jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin], RAD2DEG*src_corner_lon[imax+1], imin, imax, jmin, jmax);
lfound = rect_grid_search2(&imin3, &imax3, bound_lon1, bound_lon2, nxp1, src_corner_lon);
if ( imin3 != -1 && imin3 == imin ) imin3 -= 1;
if ( imax3 != -1 && imax3 == imin ) imax3 -= 1;
//printf(" %g %g imin %ld imax %ld jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin3], RAD2DEG*src_corner_lon[imax3+1], imin3, imax3, jmin, jmax);
for ( jm = jmin; jm <= jmax; ++jm )
for ( im = imin; im <= imax; ++im )
for ( im = imin3; im <= imax3; ++im )
srch_add[num_srch_cells++] = jm*nx + im;
}
......@@ -5400,8 +5406,46 @@ void boundbox_from_corners1(long ic, long nc, const double *restrict corner_lon,
{
long inc, j;
double clon, clat;
double clonx[nc], clatx[nc];
inc = ic*nc;
for ( j = 0; j < nc; ++j )
{
clonx[j] = corner_lon[inc+j];
clatx[j] = corner_lat[inc+j];
// if ( ic == 14 ) printf(" clon: %d %g\n", j ,clonx[j]);
}
for ( j = 1; j < nc; ++j )
{
if ( fabs(clonx[j] - clonx[j-1]) > PI )
{
if ( clonx[j] > clonx[j-1] )
clonx[j] -= PI2;
else
clonx[j-1] -= PI2;
}
}
clat = clatx[0];
clon = clonx[0];
bound_box[0] = clat;
bound_box[1] = clat;
bound_box[2] = clon;
bound_box[3] = clon;
for ( j = 1; j < nc; ++j )
{
clat = clatx[j];
clon = clonx[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;
}
/*
clat = corner_lat[inc];
clon = corner_lon[inc];
......@@ -5425,6 +5469,7 @@ void boundbox_from_corners1(long ic, long nc, const double *restrict corner_lon,
bound_box[2] = 0;
bound_box[3] = PI2;
}
*/
}
#if defined(HAVE_LIBYAC)
......@@ -5627,11 +5672,11 @@ void remap_conclip(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
{
boundbox_from_corners1(tgt_grid_add, grid2_corners, tgt_grid->cell_corner_lon, tgt_grid->cell_corner_lat, tgt_cell_bound_box);
printf("bound_box %ld lon: %g %g lat: %g %g\n", tgt_grid_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
// printf("bound_box %ld lon: %g %g lat: %g %g\n", tgt_grid_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
restrict_boundbox(src_grid_bound_box, tgt_cell_bound_box);
if ( cdoVerbose )
printf("bound_box %ld lon: %g %g lat: %g %g\n", tgt_grid_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
printf("bound_box %ld lon: %g %g lat: %g %g\n", tgt_grid_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
// printf("bound_box %ld lon: %g %g lat: %g %g\n", tgt_grid_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
}
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
......@@ -5643,7 +5688,7 @@ void remap_conclip(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
if ( cdoVerbose )
printf("tgt_grid_add %ld num_srch_cells %ld\n", tgt_grid_add, num_srch_cells);
printf("tgt_grid_add %ld num_srch_cells %ld\n", tgt_grid_add, num_srch_cells);
// printf("tgt_grid_add %ld num_srch_cells %ld\n", tgt_grid_add, num_srch_cells);
if ( num_srch_cells == 0 ) continue;
......@@ -5730,6 +5775,7 @@ void remap_conclip(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv
{
if ( area[n] > 0 )
{
//printf(">>>> %d %d %g %g\n", n, srch_add[n], tgt_area, area[n]);
weight[num_weights] = area[n] / tgt_area;
srch_add[num_weights] = srch_add[n];
num_weights++;
......
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