Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
mpim-sw
cdo
Commits
2d8ae0a0
Commit
2d8ae0a0
authored
Nov 04, 2017
by
Uwe Schulzweida
Browse files
Merge declaration and initialization.
parent
b992f802
Changes
2
Hide whitespace changes
Inline
Side-by-side
src/remap_bicubic_scrip.cc
View file @
2d8ae0a0
...
...
@@ -32,12 +32,12 @@ void set_bicubic_weights(double iw, double jw, double wgts[4][4])
wgts
[
3
][
3
]
=
iw
*
(
iw
-
1.
)
*
(
iw
-
1.
)
*
jw
*
jw
*
(
jw
-
1.
);
}
int
num_src_points
(
const
int
*
restrict
mask
,
const
size_t
src_add
[
4
],
double
src_lats
[
4
]);
unsigned
num_src_points
(
const
int
*
restrict
mask
,
const
size_t
src_add
[
4
],
double
src_lats
[
4
]);
static
void
renormalize_weights
(
const
double
src_lats
[
4
],
double
wgts
[
4
][
4
])
{
int
n
;
unsigned
n
;
double
sum_wgts
=
0.0
;
/* sum of weights for normalization */
/* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
for
(
n
=
0
;
n
<
4
;
++
n
)
sum_wgts
+=
fabs
(
src_lats
[
n
]);
...
...
@@ -50,11 +50,11 @@ void renormalize_weights(const double src_lats[4], double wgts[4][4])
static
void
bicubic_warning
(
void
)
{
static
int
lwarn
=
TRUE
;
static
bool
lwarn
=
true
;
if
(
cdoVerbose
||
lwarn
)
{
lwarn
=
FALSE
;
lwarn
=
false
;
// cdoWarning("Iteration for iw,jw exceed max iteration count of %d!", remap_max_iter);
cdoWarning
(
"Bicubic interpolation failed for some grid points - used a distance-weighted average instead!"
);
}
...
...
@@ -65,7 +65,7 @@ void bicubic_remap(double* restrict tgt_point, const double* restrict src_array,
const
double
*
restrict
grad1
,
const
double
*
restrict
grad2
,
const
double
*
restrict
grad3
)
{
*
tgt_point
=
0.
;
for
(
int
n
=
0
;
n
<
4
;
++
n
)
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
*
tgt_point
+=
src_array
[
src_add
[
n
]]
*
wgts
[
n
][
0
]
+
grad1
[
src_add
[
n
]]
*
wgts
[
n
][
1
]
+
grad2
[
src_add
[
n
]]
*
wgts
[
n
][
2
]
+
...
...
@@ -81,13 +81,6 @@ void bicubic_remap(double* restrict tgt_point, const double* restrict src_array,
*/
void
scrip_remap_bicubic_weights
(
remapgrid_t
*
src_grid
,
remapgrid_t
*
tgt_grid
,
remapvars_t
*
rv
)
{
/* Local variables */
int
search_result
;
size_t
src_add
[
4
];
/* address for the four source points */
double
src_lats
[
4
];
/* latitudes of four bilinear corners */
double
src_lons
[
4
];
/* longitudes of four bilinear corners */
double
wgts
[
4
][
4
];
/* bicubic weights for four corners */
double
plat
,
plon
;
/* lat/lon coords of destination point */
extern
int
timer_remap_bic
;
int
remap_grid_type
=
src_grid
->
remap_grid_type
;
...
...
@@ -102,7 +95,7 @@ void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
if
(
src_grid
->
rank
!=
2
)
cdoAbort
(
"Can not do bicubic interpolation when source grid rank != 2"
);
long
tgt_grid_size
=
tgt_grid
->
size
;
size_t
tgt_grid_size
=
tgt_grid
->
size
;
weightlinks4_t
*
weightlinks
=
(
weightlinks4_t
*
)
Malloc
(
tgt_grid_size
*
sizeof
(
weightlinks4_t
));
weightlinks
[
0
].
addweights
=
(
addweight4_t
*
)
Malloc
(
4
*
tgt_grid_size
*
sizeof
(
addweight4_t
));
...
...
@@ -115,10 +108,9 @@ void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
private(src_add, src_lats, src_lons, wgts, plat, plon, search_result)
shared(weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex)
#endif
for
(
long
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
for
(
size_t
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
{
#if defined(_OPENMP)
#include
"pragma_omp_atomic_update.h"
...
...
@@ -130,10 +122,16 @@ void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
if
(
!
tgt_grid
->
mask
[
tgt_cell_add
]
)
continue
;
plat
=
tgt_grid
->
cell_center_lat
[
tgt_cell_add
];
plon
=
tgt_grid
->
cell_center_lon
[
tgt_cell_add
];
double
plat
=
tgt_grid
->
cell_center_lat
[
tgt_cell_add
];
double
plon
=
tgt_grid
->
cell_center_lon
[
tgt_cell_add
];
size_t
src_add
[
4
];
// address for the four source points
double
src_lats
[
4
];
// latitudes of four bilinear corners
double
src_lons
[
4
];
// longitudes of four bilinear corners
double
wgts
[
4
][
4
];
// bicubic weights for four corners
/* Find nearest square of grid points on source grid */
int
search_result
;
if
(
remap_grid_type
==
REMAP_GRID_TYPE_REG2D
)
search_result
=
grid_search_reg2d
(
src_grid
,
src_add
,
src_lats
,
src_lons
,
plat
,
plon
,
src_grid
->
dims
,
...
...
@@ -147,17 +145,16 @@ void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
/* Check to see if points are land points */
if
(
search_result
>
0
)
{
for
(
int
n
=
0
;
n
<
4
;
++
n
)
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
if
(
!
src_grid
->
mask
[
src_add
[
n
]]
)
search_result
=
0
;
}
/* If point found, find local iw,jw coordinates for weights */
if
(
search_result
>
0
)
{
double
iw
,
jw
;
/* current guess for bilinear coordinate */
tgt_grid
->
cell_frac
[
tgt_cell_add
]
=
1.
;
double
iw
,
jw
;
/* current guess for bilinear coordinate */
if
(
find_ij_weights
(
plon
,
plat
,
src_lats
,
src_lons
,
&
iw
,
&
jw
)
)
{
/* Successfully found iw,jw - compute weights */
...
...
@@ -207,20 +204,13 @@ void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, r
*/
void
scrip_remap_bicubic
(
remapgrid_t
*
src_grid
,
remapgrid_t
*
tgt_grid
,
const
double
*
restrict
src_array
,
double
*
restrict
tgt_array
,
double
missval
)
{
/* Local variables */
int
search_result
;
size_t
src_add
[
4
];
/* address for the four source points */
double
src_lats
[
4
];
/* latitudes of four bilinear corners */
double
src_lons
[
4
];
/* longitudes of four bilinear corners */
double
wgts
[
4
][
4
];
/* bicubic weights for four corners */
double
plat
,
plon
;
/* lat/lon coords of destination point */
int
remap_grid_type
=
src_grid
->
remap_grid_type
;
if
(
cdoVerbose
)
cdoPrint
(
"Called %s()"
,
__func__
);
progressInit
();
long
tgt_grid_size
=
tgt_grid
->
size
;
size_t
tgt_grid_size
=
tgt_grid
->
size
;
/* Compute mappings from source to target grid */
...
...
@@ -239,10 +229,9 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, grad1_lat, grad1_lon, grad1_latlon, findex) \
private(src_add, src_lats, src_lons, wgts, plat, plon, search_result)
shared(remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, grad1_lat, grad1_lon, grad1_latlon, findex)
#endif
for
(
long
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
for
(
size_t
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
{
#if defined(_OPENMP)
#include
"pragma_omp_atomic_update.h"
...
...
@@ -254,10 +243,16 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
if
(
!
tgt_grid
->
mask
[
tgt_cell_add
]
)
continue
;
plat
=
tgt_grid
->
cell_center_lat
[
tgt_cell_add
];
plon
=
tgt_grid
->
cell_center_lon
[
tgt_cell_add
];
double
plat
=
tgt_grid
->
cell_center_lat
[
tgt_cell_add
];
double
plon
=
tgt_grid
->
cell_center_lon
[
tgt_cell_add
];
size_t
src_add
[
4
];
// address for the four source points
double
src_lats
[
4
];
// latitudes of four bilinear corners
double
src_lons
[
4
];
// longitudes of four bilinear corners
double
wgts
[
4
][
4
];
// bicubic weights for four corners
/* Find nearest square of grid points on source grid */
int
search_result
;
if
(
remap_grid_type
==
REMAP_GRID_TYPE_REG2D
)
search_result
=
grid_search_reg2d
(
src_grid
,
src_add
,
src_lats
,
src_lons
,
plat
,
plon
,
src_grid
->
dims
,
...
...
@@ -271,17 +266,16 @@ void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const dou
/* Check to see if points are land points */
if
(
search_result
>
0
)
{
for
(
int
n
=
0
;
n
<
4
;
++
n
)
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
if
(
!
src_grid
->
mask
[
src_add
[
n
]]
)
search_result
=
0
;
}
/* If point found, find local iw,jw coordinates for weights */
if
(
search_result
>
0
)
{
double
iw
,
jw
;
/* current guess for bilinear coordinate */
tgt_grid
->
cell_frac
[
tgt_cell_add
]
=
1.
;
double
iw
,
jw
;
/* current guess for bilinear coordinate */
if
(
find_ij_weights
(
plon
,
plat
,
src_lats
,
src_lons
,
&
iw
,
&
jw
)
)
{
/* Successfully found iw,jw - compute weights */
...
...
src/remap_bilinear_scrip.cc
View file @
2d8ae0a0
...
...
@@ -92,11 +92,11 @@ void set_bilinear_weights(double iw, double jw, double wgts[4])
}
int
num_src_points
(
const
int
*
restrict
mask
,
const
size_t
src_add
[
4
],
double
src_lats
[
4
])
unsigned
num_src_points
(
const
int
*
restrict
mask
,
const
size_t
src_add
[
4
],
double
src_lats
[
4
])
{
int
icount
=
0
;
unsigned
icount
=
0
;
for
(
int
n
=
0
;
n
<
4
;
++
n
)
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
{
if
(
mask
[
src_add
[
n
]]
)
icount
++
;
...
...
@@ -112,8 +112,8 @@ void renormalize_weights(const double src_lats[4], double wgts[4])
{
double
sum_wgts
=
0.0
;
/* sum of weights for normalization */
/* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
for
(
int
n
=
0
;
n
<
4
;
++
n
)
sum_wgts
+=
fabs
(
src_lats
[
n
]);
for
(
int
n
=
0
;
n
<
4
;
++
n
)
wgts
[
n
]
=
fabs
(
src_lats
[
n
])
/
sum_wgts
;
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
sum_wgts
+=
fabs
(
src_lats
[
n
]);
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
wgts
[
n
]
=
fabs
(
src_lats
[
n
])
/
sum_wgts
;
}
static
...
...
@@ -148,7 +148,7 @@ static
void
bilinear_remap
(
double
*
restrict
tgt_point
,
const
double
*
restrict
src_array
,
const
double
wgts
[
4
],
const
size_t
src_add
[
4
])
{
// *tgt_point = 0.;
// for (
int
n = 0; n < 4; ++n ) *tgt_point += src_array[src_add[n]]*wgts[n];
// for (
unsigned
n = 0; n < 4; ++n ) *tgt_point += src_array[src_add[n]]*wgts[n];
*
tgt_point
=
src_array
[
src_add
[
0
]]
*
wgts
[
0
]
+
src_array
[
src_add
[
1
]]
*
wgts
[
1
]
+
src_array
[
src_add
[
2
]]
*
wgts
[
2
]
+
src_array
[
src_add
[
3
]]
*
wgts
[
3
];
}
...
...
@@ -162,12 +162,6 @@ void bilinear_remap(double* restrict tgt_point, const double *restrict src_array
*/
void
scrip_remap_bilinear_weights
(
remapgrid_t
*
src_grid
,
remapgrid_t
*
tgt_grid
,
remapvars_t
*
rv
)
{
/* Local variables */
int
search_result
;
size_t
src_add
[
4
];
/* address for the four source points */
double
src_lats
[
4
];
/* latitudes of four bilinear corners */
double
src_lons
[
4
];
/* longitudes of four bilinear corners */
double
wgts
[
4
];
/* bilinear weights for four corners */
extern
int
timer_remap_bil
;
int
remap_grid_type
=
src_grid
->
remap_grid_type
;
...
...
@@ -182,11 +176,11 @@ void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid,
if
(
src_grid
->
rank
!=
2
)
cdoAbort
(
"Can not do bilinear interpolation when source grid rank != 2"
);
long
tgt_grid_size
=
tgt_grid
->
size
;
size_t
tgt_grid_size
=
tgt_grid
->
size
;
weightlinks_t
*
weightlinks
=
(
weightlinks_t
*
)
Malloc
(
tgt_grid_size
*
sizeof
(
weightlinks_t
));
weightlinks
[
0
].
addweights
=
(
addweight_t
*
)
Malloc
(
4
*
tgt_grid_size
*
sizeof
(
addweight_t
));
for
(
unsigned
tgt_cell_add
=
1
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
for
(
size_t
tgt_cell_add
=
1
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
weightlinks
[
tgt_cell_add
].
addweights
=
weightlinks
[
0
].
addweights
+
4
*
tgt_cell_add
;
double
findex
=
0
;
...
...
@@ -194,12 +188,10 @@ void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid,
/* Loop over destination grid */
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
private(src_add, src_lats, src_lons, wgts, search_result) \
schedule(static)
#pragma omp parallel for default(none) schedule(static) \
shared(weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex)
#endif
for
(
long
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
for
(
size_t
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
{
#if defined(_OPENMP)
#include
"pragma_omp_atomic_update.h"
...
...
@@ -214,7 +206,13 @@ void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid,
double
plon
=
0
,
plat
=
0
;
remapgrid_get_lonlat
(
tgt_grid
,
tgt_cell_add
,
&
plon
,
&
plat
);
size_t
src_add
[
4
];
// address for the four source points
double
src_lats
[
4
];
// latitudes of four bilinear corners
double
src_lons
[
4
];
// longitudes of four bilinear corners
double
wgts
[
4
];
// bilinear weights for four corners
// Find nearest square of grid points on source grid
int
search_result
;
if
(
remap_grid_type
==
REMAP_GRID_TYPE_REG2D
)
search_result
=
grid_search_reg2d
(
src_grid
,
src_add
,
src_lats
,
src_lons
,
plat
,
plon
,
src_grid
->
dims
,
...
...
@@ -228,17 +226,16 @@ void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid,
// Check to see if points are mask points
if
(
search_result
>
0
)
{
for
(
int
n
=
0
;
n
<
4
;
++
n
)
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
if
(
!
src_grid
->
mask
[
src_add
[
n
]]
)
search_result
=
0
;
}
// If point found, find local iw,jw coordinates for weights
if
(
search_result
>
0
)
{
double
iw
,
jw
;
// current guess for bilinear coordinate
tgt_grid
->
cell_frac
[
tgt_cell_add
]
=
1.
;
double
iw
,
jw
;
// current guess for bilinear coordinate
if
(
find_ij_weights
(
plon
,
plat
,
src_lats
,
src_lons
,
&
iw
,
&
jw
)
)
{
// Successfully found iw,jw - compute weights
...
...
@@ -287,12 +284,6 @@ void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid,
*/
void
scrip_remap_bilinear
(
remapgrid_t
*
src_grid
,
remapgrid_t
*
tgt_grid
,
const
double
*
restrict
src_array
,
double
*
restrict
tgt_array
,
double
missval
)
{
/* Local variables */
int
search_result
;
size_t
src_add
[
4
];
/* address for the four source points */
double
src_lats
[
4
];
/* latitudes of four bilinear corners */
double
src_lons
[
4
];
/* longitudes of four bilinear corners */
double
wgts
[
4
];
/* bilinear weights for four corners */
extern
int
timer_remap_bil
;
int
remap_grid_type
=
src_grid
->
remap_grid_type
;
...
...
@@ -302,7 +293,7 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
progressInit
();
long
tgt_grid_size
=
tgt_grid
->
size
;
size_t
tgt_grid_size
=
tgt_grid
->
size
;
/* Compute mappings from source to target grid */
...
...
@@ -314,12 +305,10 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
/* Loop over destination grid */
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex) \
private(src_add, src_lats, src_lons, wgts, search_result) \
schedule(static)
#pragma omp parallel for default(none) schedule(static) \
shared(cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex)
#endif
for
(
long
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
for
(
size_t
tgt_cell_add
=
0
;
tgt_cell_add
<
tgt_grid_size
;
++
tgt_cell_add
)
{
#if defined(_OPENMP)
#include
"pragma_omp_atomic_update.h"
...
...
@@ -334,7 +323,13 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
double
plon
=
0
,
plat
=
0
;
remapgrid_get_lonlat
(
tgt_grid
,
tgt_cell_add
,
&
plon
,
&
plat
);
size_t
src_add
[
4
];
// address for the four source points
double
src_lats
[
4
];
// latitudes of four bilinear corners
double
src_lons
[
4
];
// longitudes of four bilinear corners
double
wgts
[
4
];
// bilinear weights for four corners
// Find nearest square of grid points on source grid
int
search_result
;
if
(
remap_grid_type
==
REMAP_GRID_TYPE_REG2D
)
search_result
=
grid_search_reg2d
(
src_grid
,
src_add
,
src_lats
,
src_lons
,
plat
,
plon
,
src_grid
->
dims
,
...
...
@@ -348,17 +343,16 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
// Check to see if points are mask points
if
(
search_result
>
0
)
{
for
(
int
n
=
0
;
n
<
4
;
++
n
)
for
(
unsigned
n
=
0
;
n
<
4
;
++
n
)
if
(
!
src_grid
->
mask
[
src_add
[
n
]]
)
search_result
=
0
;
}
// If point found, find local iw,jw coordinates for weights
if
(
search_result
>
0
)
{
double
iw
,
jw
;
// current guess for bilinear coordinate
tgt_grid
->
cell_frac
[
tgt_cell_add
]
=
1.
;
double
iw
,
jw
;
// current guess for bilinear coordinate
if
(
find_ij_weights
(
plon
,
plat
,
src_lats
,
src_lons
,
&
iw
,
&
jw
)
)
{
// Successfully found iw,jw - compute weights
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment