diff --git a/raster/r.resamp.interp/Makefile b/raster/r.resamp.interp/Makefile index ba9472ca18a..59546ad50ba 100644 --- a/raster/r.resamp.interp/Makefile +++ b/raster/r.resamp.interp/Makefile @@ -2,8 +2,9 @@ MODULE_TOPDIR = ../.. PGM = r.resamp.interp -LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) +LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(OMPLIB) DEPENDENCIES = $(RASTERDEP) $(GISDEP) +EXTRA_CFLAGS = $(OMPCFLAGS) include $(MODULE_TOPDIR)/include/Make/Module.make diff --git a/raster/r.resamp.interp/benchmark/benchmark_r_resamp_interp.py b/raster/r.resamp.interp/benchmark/benchmark_r_resamp_interp.py new file mode 100644 index 00000000000..20196bc44bc --- /dev/null +++ b/raster/r.resamp.interp/benchmark/benchmark_r_resamp_interp.py @@ -0,0 +1,57 @@ +"""Benchmarking of r.resamp.interp + +@author Aaron Saw Min Sern +""" + +from grass.exceptions import CalledModuleError +from grass.pygrass.modules import Module +from subprocess import DEVNULL + +import grass.benchmark as bm + + +def main(): + results = [] + + # Users can add more or modify existing reference maps + benchmark(7071, "r.resamp.interp_50M", results) + benchmark(14142, "r.resamp.interp_200M", results) + benchmark(10000, "r.resamp.interp_100M", results) + benchmark(20000, "r.resamp.interp_400M", results) + + bm.nprocs_plot(results) + + +def benchmark(size, label, results): + reference = "r_resamp_interp_reference_map" + output = "benchmark_r_resamp_interp" + + generate_map(rows=size, cols=size, fname=reference) + module = Module( + "r.resamp.interp", + input=reference, + output=output, + method="lanczos", + memory=500, + run_=False, + stdout_=DEVNULL, + overwrite=True, + ) + results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=16, repeat=3)) + Module("g.remove", quiet=True, flags="f", type="raster", name=reference) + Module("g.remove", quiet=True, flags="f", type="raster", name=output) + + +def generate_map(rows, cols, fname): + Module("g.region", flags="p", s=0, n=rows, w=0, e=cols, res=1) + # Generate using r.random.surface if r.surf.fractal fails + try: + print("Generating reference map using r.surf.fractal...") + Module("r.surf.fractal", output=fname) + except CalledModuleError: + print("r.surf.fractal fails, using r.random.surface instead...") + Module("r.random.surface", output=fname) + + +if __name__ == "__main__": + main() diff --git a/raster/r.resamp.interp/main.c b/raster/r.resamp.interp/main.c index e57c388590c..1c1d2360e1c 100644 --- a/raster/r.resamp.interp/main.c +++ b/raster/r.resamp.interp/main.c @@ -6,8 +6,9 @@ * Paul Kelly , * Dylan Beaudette, Hamish Bowman , * Markus Metz (lanczos) + * Aaron Saw Min Sern (OpenMP parallelization) * PURPOSE: - * COPYRIGHT: (C) 2006-2007 by the GRASS Development Team + * COPYRIGHT: (C) 2006-2022 by the GRASS Development Team * * This program is free software under the GNU General Public * License (>=v2). Read the file COPYING that comes with GRASS @@ -17,23 +18,26 @@ /* TODO: add fallback methods, see e.g. r.proj */ +#if defined(_OPENMP) +#include +#endif #include #include #include #include #include #include +#define FIRST_THREAD 0 static int neighbors; -static DCELL *bufs[5]; -static int cur_row; +static DCELL *(*bufs)[5]; -static void read_rows(int infile, int row) +static void read_rows(int infile, int row, int t_id, int *cur_row) { - int end_row = cur_row + neighbors; + int end_row = *cur_row + neighbors; int first_row = row; int last_row = row + neighbors; - int offset = first_row - cur_row; + int offset = first_row - *cur_row; int keep = end_row - first_row; int i; @@ -42,32 +46,35 @@ static void read_rows(int infile, int row) if (keep > 0 && offset > 0) for (i = 0; i < keep; i++) { - DCELL *tmp = bufs[i]; + DCELL *tmp = bufs[t_id][i]; - bufs[i] = bufs[i + offset]; - bufs[i + offset] = tmp; + bufs[t_id][i] = bufs[t_id][i + offset]; + bufs[t_id][i + offset] = tmp; } if (keep < 0) keep = 0; for (i = keep; i < neighbors; i++) - Rast_get_d_row(infile, bufs[i], first_row + i); + Rast_get_d_row(infile, bufs[t_id][i], first_row + i); - cur_row = first_row; + *cur_row = first_row; } int main(int argc, char *argv[]) { struct GModule *module; - struct Option *rastin, *rastout, *method; + struct Option *rastin, *rastout, *method, *nprocs, *memory; struct History history; char title[64]; char buf_nsres[100], buf_ewres[100]; struct Colors colors; - int infile, outfile; + int *infile; + int outfile; DCELL *outbuf; - int row, col; + int bufrows; + int threads; + int t, row, col; struct Cell_head dst_w, src_w; G_gisinit(argv[0]); @@ -80,6 +87,7 @@ int main(int argc, char *argv[]) G_add_keyword(_("bilinear")); G_add_keyword(_("bicubic")); G_add_keyword(_("lanczos")); + G_add_keyword(_("parallel")); module->description = _("Resamples raster map to a finer grid using interpolation."); @@ -91,6 +99,9 @@ int main(int argc, char *argv[]) method->answer = "bilinear"; method->guisection = _("Method"); + nprocs = G_define_standard_option(G_OPT_M_NPROCS); + memory = G_define_standard_option(G_OPT_MEMORYMB); + if (G_parser(argc, argv)) exit(EXIT_FAILURE); @@ -107,6 +118,30 @@ int main(int argc, char *argv[]) G_get_set_window(&dst_w); + sscanf(nprocs->answer, "%d", &threads); + if (threads < 1) { + G_fatal_error(_("<%d> is not valid number of threads."), threads); + } +#if defined(_OPENMP) + omp_set_num_threads(threads); +#else + if (threads != 1) + G_warning(_("GRASS is compiled without OpenMP support. Ignoring " + "threads setting.")); + threads = 1; +#endif + + bufrows = + atoi(memory->answer) * (((1 << 20) / sizeof(DCELL)) / dst_w.cols); + /* set the output buffer rows to be at most covering the entire map */ + if (bufrows > dst_w.rows) { + bufrows = dst_w.rows; + } + /* but at least the number of threads */ + if (bufrows < threads) { + bufrows = threads; + } + /* set window to old map */ Rast_get_cellhd(rastin->answer, "", &src_w); @@ -150,212 +185,275 @@ int main(int argc, char *argv[]) Rast_set_input_window(&src_w); /* allocate buffers for input rows */ - for (row = 0; row < neighbors; row++) - bufs[row] = Rast_allocate_d_input_buf(); - - cur_row = -100; + bufs = G_malloc(threads * sizeof *bufs); + for (t = 0; t < threads; t++) + for (row = 0; row < neighbors; row++) + bufs[t][row] = Rast_allocate_d_input_buf(); /* open old map */ - infile = Rast_open_old(rastin->answer, ""); + infile = G_malloc(threads * sizeof *infile); + for (t = 0; t < threads; t++) + infile[t] = Rast_open_old(rastin->answer, ""); /* reset window to current region */ Rast_set_output_window(&dst_w); - outbuf = Rast_allocate_d_output_buf(); + outbuf = G_calloc((size_t)bufrows * dst_w.cols, sizeof(DCELL)); /* open new map */ outfile = Rast_open_new(rastout->answer, DCELL_TYPE); - switch (neighbors) { - case 1: /* nearest */ - for (row = 0; row < dst_w.rows; row++) { - double north = Rast_row_to_northing(row + 0.5, &dst_w); - double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5; - int maprow0 = (int)floor(maprow_f + 0.5); + int computed = 0; + int written = 0; + int cur_row = -100; - G_percent(row, dst_w.rows, 2); + while (written < dst_w.rows) { + int range = bufrows; - read_rows(infile, maprow0); + if (range > dst_w.rows - written) + range = dst_w.rows - written; + int start = written; + int end = written + range; - for (col = 0; col < dst_w.cols; col++) { - double east = Rast_col_to_easting(col + 0.5, &dst_w); - double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5; - int mapcol0 = (int)floor(mapcol_f + 0.5); +#pragma omp parallel private(row, col) firstprivate(cur_row) if(threads > 1) + { + int t_id = FIRST_THREAD; - double c = bufs[0][mapcol0]; +#if defined(_OPENMP) + t_id = omp_get_thread_num(); +#endif - if (Rast_is_d_null_value(&c)) { - Rast_set_d_null_value(&outbuf[col], 1); - } - else { - outbuf[col] = c; - } - } + switch (neighbors) { + case 1: /* nearest */ +#pragma omp for schedule(static) + for (row = start; row < end; row++) { + double north = Rast_row_to_northing(row + 0.5, &dst_w); + double maprow_f = + Rast_northing_to_row(north, &src_w) - 0.5; + int maprow0 = (int)floor(maprow_f + 0.5); - Rast_put_d_row(outfile, outbuf); - } - break; - - case 2: /* bilinear */ - for (row = 0; row < dst_w.rows; row++) { - double north = Rast_row_to_northing(row + 0.5, &dst_w); - double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5; - int maprow0 = (int)floor(maprow_f); - double v = maprow_f - maprow0; - - G_percent(row, dst_w.rows, 2); - - read_rows(infile, maprow0); - - for (col = 0; col < dst_w.cols; col++) { - double east = Rast_col_to_easting(col + 0.5, &dst_w); - double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5; - int mapcol0 = (int)floor(mapcol_f); - int mapcol1 = mapcol0 + 1; - double u = mapcol_f - mapcol0; - - double c00 = bufs[0][mapcol0]; - double c01 = bufs[0][mapcol1]; - double c10 = bufs[1][mapcol0]; - double c11 = bufs[1][mapcol1]; - - if (Rast_is_d_null_value(&c00) || - Rast_is_d_null_value(&c01) || - Rast_is_d_null_value(&c10) || - Rast_is_d_null_value(&c11)) { - Rast_set_d_null_value(&outbuf[col], 1); - } - else { - outbuf[col] = - Rast_interp_bilinear(u, v, c00, c01, c10, c11); + G_percent(computed, dst_w.rows, 2); + + read_rows(infile[t_id], maprow0, t_id, &cur_row); + + for (col = 0; col < dst_w.cols; col++) { + double east = Rast_col_to_easting(col + 0.5, &dst_w); + double mapcol_f = + Rast_easting_to_col(east, &src_w) - 0.5; + int mapcol0 = (int)floor(mapcol_f + 0.5); + + double c = bufs[t_id][0][mapcol0]; + + if (Rast_is_d_null_value(&c)) { + Rast_set_d_null_value(&outbuf + [(row - + start) * dst_w.cols + + col], 1); + } + else { + outbuf[(row - start) * dst_w.cols + col] = c; + } + } + +#pragma omp atomic update + computed++; } - } + break; + + case 2: /* bilinear */ +#pragma omp for schedule(static) + for (row = start; row < end; row++) { + double north = Rast_row_to_northing(row + 0.5, &dst_w); + double maprow_f = + Rast_northing_to_row(north, &src_w) - 0.5; + int maprow0 = (int)floor(maprow_f); + double v = maprow_f - maprow0; + + G_percent(computed, dst_w.rows, 2); + + read_rows(infile[t_id], maprow0, t_id, &cur_row); + + for (col = 0; col < dst_w.cols; col++) { + double east = Rast_col_to_easting(col + 0.5, &dst_w); + double mapcol_f = + Rast_easting_to_col(east, &src_w) - 0.5; + int mapcol0 = (int)floor(mapcol_f); + int mapcol1 = mapcol0 + 1; + double u = mapcol_f - mapcol0; + + double c00 = bufs[t_id][0][mapcol0]; + double c01 = bufs[t_id][0][mapcol1]; + double c10 = bufs[t_id][1][mapcol0]; + double c11 = bufs[t_id][1][mapcol1]; + + if (Rast_is_d_null_value(&c00) || + Rast_is_d_null_value(&c01) || + Rast_is_d_null_value(&c10) || + Rast_is_d_null_value(&c11)) { + Rast_set_d_null_value(&outbuf + [(row - + start) * dst_w.cols + + col], 1); + } + else { + outbuf[(row - start) * dst_w.cols + col] = + Rast_interp_bilinear(u, v, c00, c01, c10, + c11); + } + } - Rast_put_d_row(outfile, outbuf); - } - break; - - case 4: /* bicubic */ - for (row = 0; row < dst_w.rows; row++) { - double north = Rast_row_to_northing(row + 0.5, &dst_w); - double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5; - int maprow1 = (int)floor(maprow_f); - int maprow0 = maprow1 - 1; - double v = maprow_f - maprow1; - - G_percent(row, dst_w.rows, 2); - - read_rows(infile, maprow0); - - for (col = 0; col < dst_w.cols; col++) { - double east = Rast_col_to_easting(col + 0.5, &dst_w); - double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5; - int mapcol1 = (int)floor(mapcol_f); - int mapcol0 = mapcol1 - 1; - int mapcol2 = mapcol1 + 1; - int mapcol3 = mapcol1 + 2; - double u = mapcol_f - mapcol1; - - double c00 = bufs[0][mapcol0]; - double c01 = bufs[0][mapcol1]; - double c02 = bufs[0][mapcol2]; - double c03 = bufs[0][mapcol3]; - - double c10 = bufs[1][mapcol0]; - double c11 = bufs[1][mapcol1]; - double c12 = bufs[1][mapcol2]; - double c13 = bufs[1][mapcol3]; - - double c20 = bufs[2][mapcol0]; - double c21 = bufs[2][mapcol1]; - double c22 = bufs[2][mapcol2]; - double c23 = bufs[2][mapcol3]; - - double c30 = bufs[3][mapcol0]; - double c31 = bufs[3][mapcol1]; - double c32 = bufs[3][mapcol2]; - double c33 = bufs[3][mapcol3]; - - if (Rast_is_d_null_value(&c00) || - Rast_is_d_null_value(&c01) || - Rast_is_d_null_value(&c02) || - Rast_is_d_null_value(&c03) || - Rast_is_d_null_value(&c10) || - Rast_is_d_null_value(&c11) || - Rast_is_d_null_value(&c12) || - Rast_is_d_null_value(&c13) || - Rast_is_d_null_value(&c20) || - Rast_is_d_null_value(&c21) || - Rast_is_d_null_value(&c22) || - Rast_is_d_null_value(&c23) || - Rast_is_d_null_value(&c30) || - Rast_is_d_null_value(&c31) || - Rast_is_d_null_value(&c32) || - Rast_is_d_null_value(&c33)) { - Rast_set_d_null_value(&outbuf[col], 1); +#pragma omp atomic update + computed++; } - else { - outbuf[col] = Rast_interp_bicubic(u, v, - c00, c01, c02, c03, - c10, c11, c12, c13, - c20, c21, c22, c23, - c30, c31, c32, c33); + break; + + case 4: /* bicubic */ +#pragma omp for schedule(static) + for (row = start; row < end; row++) { + double north = Rast_row_to_northing(row + 0.5, &dst_w); + double maprow_f = + Rast_northing_to_row(north, &src_w) - 0.5; + int maprow1 = (int)floor(maprow_f); + int maprow0 = maprow1 - 1; + double v = maprow_f - maprow1; + + G_percent(computed, dst_w.rows, 2); + + read_rows(infile[t_id], maprow0, t_id, &cur_row); + + for (col = 0; col < dst_w.cols; col++) { + double east = Rast_col_to_easting(col + 0.5, &dst_w); + double mapcol_f = + Rast_easting_to_col(east, &src_w) - 0.5; + int mapcol1 = (int)floor(mapcol_f); + int mapcol0 = mapcol1 - 1; + int mapcol2 = mapcol1 + 1; + int mapcol3 = mapcol1 + 2; + double u = mapcol_f - mapcol1; + + double c00 = bufs[t_id][0][mapcol0]; + double c01 = bufs[t_id][0][mapcol1]; + double c02 = bufs[t_id][0][mapcol2]; + double c03 = bufs[t_id][0][mapcol3]; + + double c10 = bufs[t_id][1][mapcol0]; + double c11 = bufs[t_id][1][mapcol1]; + double c12 = bufs[t_id][1][mapcol2]; + double c13 = bufs[t_id][1][mapcol3]; + + double c20 = bufs[t_id][2][mapcol0]; + double c21 = bufs[t_id][2][mapcol1]; + double c22 = bufs[t_id][2][mapcol2]; + double c23 = bufs[t_id][2][mapcol3]; + + double c30 = bufs[t_id][3][mapcol0]; + double c31 = bufs[t_id][3][mapcol1]; + double c32 = bufs[t_id][3][mapcol2]; + double c33 = bufs[t_id][3][mapcol3]; + + if (Rast_is_d_null_value(&c00) || + Rast_is_d_null_value(&c01) || + Rast_is_d_null_value(&c02) || + Rast_is_d_null_value(&c03) || + Rast_is_d_null_value(&c10) || + Rast_is_d_null_value(&c11) || + Rast_is_d_null_value(&c12) || + Rast_is_d_null_value(&c13) || + Rast_is_d_null_value(&c20) || + Rast_is_d_null_value(&c21) || + Rast_is_d_null_value(&c22) || + Rast_is_d_null_value(&c23) || + Rast_is_d_null_value(&c30) || + Rast_is_d_null_value(&c31) || + Rast_is_d_null_value(&c32) || + Rast_is_d_null_value(&c33)) { + Rast_set_d_null_value(&outbuf + [(row - + start) * dst_w.cols + + col], 1); + } + else { + outbuf[(row - start) * dst_w.cols + col] = + Rast_interp_bicubic(u, v, c00, c01, c02, c03, + c10, c11, c12, c13, c20, + c21, c22, c23, c30, c31, + c32, c33); + } + } + +#pragma omp atomic update + computed++; } - } + break; + + case 5: /* lanczos */ +#pragma omp for schedule(static) + for (row = start; row < end; row++) { + double north = Rast_row_to_northing(row + 0.5, &dst_w); + double maprow_f = + Rast_northing_to_row(north, &src_w) - 0.5; + int maprow1 = (int)floor(maprow_f + 0.5); + int maprow0 = maprow1 - 2; + double v = maprow_f - maprow1; + + G_percent(computed, dst_w.rows, 2); + + read_rows(infile[t_id], maprow0, t_id, &cur_row); + + for (col = 0; col < dst_w.cols; col++) { + double east = Rast_col_to_easting(col + 0.5, &dst_w); + double mapcol_f = + Rast_easting_to_col(east, &src_w) - 0.5; + int mapcol2 = (int)floor(mapcol_f + 0.5); + int mapcol0 = mapcol2 - 2; + int mapcol4 = mapcol2 + 2; + double u = mapcol_f - mapcol2; + double c[25]; + int ci = 0, i, j, do_lanczos = 1; + + for (i = 0; i < 5; i++) { + for (j = mapcol0; j <= mapcol4; j++) { + c[ci] = bufs[t_id][i][j]; + if (Rast_is_d_null_value(&(c[ci]))) { + Rast_set_d_null_value(&outbuf + [(row - + start) * + dst_w.cols + col], + 1); + do_lanczos = 0; + break; + } + ci++; + } + if (!do_lanczos) + break; + } - Rast_put_d_row(outfile, outbuf); - } - break; - - case 5: /* lanczos */ - for (row = 0; row < dst_w.rows; row++) { - double north = Rast_row_to_northing(row + 0.5, &dst_w); - double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5; - int maprow1 = (int)floor(maprow_f + 0.5); - int maprow0 = maprow1 - 2; - double v = maprow_f - maprow1; - - G_percent(row, dst_w.rows, 2); - - read_rows(infile, maprow0); - - for (col = 0; col < dst_w.cols; col++) { - double east = Rast_col_to_easting(col + 0.5, &dst_w); - double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5; - int mapcol2 = (int)floor(mapcol_f + 0.5); - int mapcol0 = mapcol2 - 2; - int mapcol4 = mapcol2 + 2; - double u = mapcol_f - mapcol2; - double c[25]; - int ci = 0, i, j, do_lanczos = 1; - - for (i = 0; i < 5; i++) { - for (j = mapcol0; j <= mapcol4; j++) { - c[ci] = bufs[i][j]; - if (Rast_is_d_null_value(&(c[ci]))) { - Rast_set_d_null_value(&outbuf[col], 1); - do_lanczos = 0; - break; + if (do_lanczos) { + outbuf[(row - start) * dst_w.cols + col] = + Rast_interp_lanczos(u, v, c); } - ci++; } - if (!do_lanczos) - break; - } - if (do_lanczos) { - outbuf[col] = Rast_interp_lanczos(u, v, c); +#pragma omp atomic update + computed++; } + break; } + } - Rast_put_d_row(outfile, outbuf); + /* write to output map */ + for (row = start; row < end; row++) { + Rast_put_d_row(outfile, &outbuf[(row - start) * dst_w.cols]); } - break; + written = end; + } G_percent(dst_w.rows, dst_w.rows, 2); - Rast_close(infile); + for (t = 0; t < threads; t++) + Rast_close(infile[t]); Rast_close(outfile); diff --git a/raster/r.resamp.interp/r.resamp.interp.html b/raster/r.resamp.interp/r.resamp.interp.html index bb7d6fc7449..b8cb239891d 100644 --- a/raster/r.resamp.interp/r.resamp.interp.html +++ b/raster/r.resamp.interp/r.resamp.interp.html @@ -36,6 +36,20 @@

NOTES

attempt to implement the latter would violate the integrity of the interpolation method. +

PERFORMANCE

+

By specifying the number of parallel processes with nprocs option, +r.resamp.interp can run significantly faster, see benchmarks below. + +

+ benchmark for number of cells +
+ Figure: Benchmark shows execution time for different + number of cells. See benchmark scripts in source code. +
+

To reduce the memory requirements to minimum, set option memory to zero. +To take advantage of the parallelization, GRASS GIS +needs to compiled with OpenMP enabled. +

EXAMPLE

Resample elevation raster map to a higher resolution (from 500m to 250m; diff --git a/raster/r.resamp.interp/r_resamp_interp_benchmark_size.png b/raster/r.resamp.interp/r_resamp_interp_benchmark_size.png new file mode 100644 index 00000000000..b21ef6126db Binary files /dev/null and b/raster/r.resamp.interp/r_resamp_interp_benchmark_size.png differ diff --git a/raster/r.resamp.interp/testsuite/test_r_resamp_interp.py b/raster/r.resamp.interp/testsuite/test_r_resamp_interp.py new file mode 100644 index 00000000000..eb1f14ac849 --- /dev/null +++ b/raster/r.resamp.interp/testsuite/test_r_resamp_interp.py @@ -0,0 +1,137 @@ +from grass.gunittest.case import TestCase +from grass.gunittest.main import test + + +class TestResampInterp(TestCase): + """ + + Used dataset: nc_spm_08_grass7 + + Test cases: + test_finite_options: Test output with finite kernels only. + test_infinite_box_options: Test output with infinite kernels + in conjunction with the box kernel. + """ + + test_options = { + "nearest": { + "n": 5625000, + "null_cells": 0, + "cells": 5625000, + "min": 55.5787925720215, + "max": 156.329864501953, + "range": 100.751071929932, + "mean": 110.375436614522, + "mean_of_abs": 110.375436614522, + "stddev": 20.3153283525752, + "variance": 412.712566072944, + "coeff_var": 18.4056606938053, + "sum": 620861830.956688, + }, + "bilinear": { + "n": 5615504, + "null_cells": 9496, + "cells": 5625000, + "min": 55.9232139587402, + "max": 156.316369018555, + "range": 100.393155059814, + "mean": 110.382539288015, + "mean_of_abs": 110.382539288015, + "stddev": 20.304736015932, + "variance": 412.282304676688, + "coeff_var": 18.3948803378694, + "sum": 619853590.902006, + }, + "bicubic": { + "n": 5601275, + "null_cells": 23725, + "cells": 5625000, + "min": 55.8029828833008, + "max": 156.356721084961, + "range": 100.55373820166, + "mean": 110.404233965113, + "mean_of_abs": 110.404233965113, + "stddev": 20.2966878636654, + "variance": 411.955538235063, + "coeff_var": 18.3839759896156, + "sum": 618404475.60294, + }, + "lanczos": { + "n": 5596536, + "null_cells": 28464, + "cells": 5625000, + "min": 55.7825176472632, + "max": 156.36299154775, + "range": 100.580473900487, + "mean": 110.396555800561, + "mean_of_abs": 110.396555800561, + "stddev": 20.2935834200418, + "variance": 411.829528026196, + "coeff_var": 18.3824425253842, + "sum": 617838298.813847, + }, + } + + # TODO: replace by unified handing of maps + to_remove = [] + + def check_univar(self, test_case): + """Checks multiple output against unviariate reference from dict""" + for method in self.test_options: + self.assertRasterFitsUnivar( + raster="{}_raster_{}".format(test_case, method), + reference=self.test_options[method], + precision=1e-5, + ) + self.assertRasterFitsUnivar( + raster="{}_raster_threaded_{}".format(test_case, method), + reference=self.test_options[method], + precision=1e-5, + ) + + @classmethod + def setUpClass(cls): + cls.use_temp_region() + # run with resolution of 6 + cls.runModule("g.region", raster="elevation", res=6) + + @classmethod + def tearDownClass(cls): + cls.del_temp_region() + if cls.to_remove: + cls.runModule( + "g.remove", flags="f", type="raster", name=",".join(cls.to_remove) + ) + + def test_resamp_interp(self): + """Test for all methods of interpolation.""" + test_case = "test_resamp_interp" + outputs = [ + "{}_raster_{}".format(test_case, filter) for filter in self.test_options + ] + outputs_threaded = [ + "{}_raster_threaded_{}".format(test_case, filter) + for filter in self.test_options + ] + self.to_remove.extend(outputs) + self.to_remove.extend(outputs_threaded) + + for method in self.test_options: + self.assertModule( + "r.resamp.interp", + input="elevation", + output="{}_raster_{}".format(test_case, method), + method=method, + ) + self.assertModule( + "r.resamp.interp", + input="elevation", + output="{}_raster_threaded_{}".format(test_case, method), + method=method, + nprocs=8, + ) + self.check_univar(test_case) + + +if __name__ == "__main__": + test()