GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /***************************************************************************** 00003 * 00004 * MODULE: Grass PDE Numerical Library 00005 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006 00006 * soerengebbert <at> gmx <dot> de 00007 * 00008 * PURPOSE: gradient management functions 00009 * part of the gpde library 00010 * 00011 * COPYRIGHT: (C) 2000 by the GRASS Development Team 00012 * 00013 * This program is free software under the GNU General Public 00014 * License (>=v2). Read the file COPYING that comes with GRASS 00015 * for details. 00016 * 00017 *****************************************************************************/ 00018 00019 #include <grass/N_pde.h> 00020 00029 void N_calc_gradient_field_2d_stats(N_gradient_field_2d * field) 00030 { 00031 double minx, miny; 00032 double maxx, maxy; 00033 double sumx, sumy; 00034 int nonullx, nonully; 00035 00036 G_debug(3, 00037 "N_calc_gradient_field_2d_stats: compute gradient field stats"); 00038 00039 N_calc_array_2d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0); 00040 N_calc_array_2d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0); 00041 00042 if (minx < miny) 00043 field->min = minx; 00044 else 00045 field->min = miny; 00046 00047 if (maxx > maxy) 00048 field->max = maxx; 00049 else 00050 field->max = maxy; 00051 00052 field->sum = sumx + sumy; 00053 field->nonull = nonullx + nonully; 00054 field->mean = field->sum / (double)field->nonull; 00055 00056 return; 00057 } 00058 00105 N_gradient_field_2d *N_compute_gradient_field_2d(N_array_2d * pot, 00106 N_array_2d * weight_x, 00107 N_array_2d * weight_y, 00108 N_geom_data * geom, 00109 N_gradient_field_2d * 00110 gradfield) 00111 { 00112 int i, j; 00113 int rows, cols; 00114 double dx, dy, p1, p2, r1, r2, mean, grad, res; 00115 N_gradient_field_2d *field = gradfield; 00116 00117 00118 if (pot->cols != weight_x->cols || pot->cols != weight_y->cols) 00119 G_fatal_error 00120 ("N_compute_gradient_field_2d: the arrays are not of equal size"); 00121 00122 if (pot->rows != weight_x->rows || pot->rows != weight_y->rows) 00123 G_fatal_error 00124 ("N_compute_gradient_field_2d: the arrays are not of equal size"); 00125 00126 if (pot->cols != geom->cols || pot->rows != geom->rows) 00127 G_fatal_error 00128 ("N_compute_gradient_field_2d: array sizes and geometry data are different"); 00129 00130 00131 G_debug(3, "N_compute_gradient_field_2d: compute gradient field"); 00132 00133 rows = pot->rows; 00134 cols = pot->cols; 00135 dx = geom->dx; 00136 dy = geom->dy; 00137 00138 if (field == NULL) { 00139 field = N_alloc_gradient_field_2d(cols, rows); 00140 } 00141 else { 00142 if (field->cols != geom->cols || field->rows != geom->rows) 00143 G_fatal_error 00144 ("N_compute_gradient_field_2d: gradient field sizes and geometry data are different"); 00145 } 00146 00147 00148 for (j = 0; j < rows; j++) 00149 for (i = 0; i < cols - 1; i++) { 00150 grad = 0; 00151 mean = 0; 00152 00153 /* Only compute if the arrays are not null */ 00154 if (!N_is_array_2d_value_null(pot, i, j) && 00155 !N_is_array_2d_value_null(pot, i + 1, j)) { 00156 p1 = N_get_array_2d_d_value(pot, i, j); 00157 p2 = N_get_array_2d_d_value(pot, i + 1, j); 00158 grad = (p1 - p2) / dx; /* gradient */ 00159 } 00160 if (!N_is_array_2d_value_null(weight_x, i, j) && 00161 !N_is_array_2d_value_null(weight_x, i + 1, j)) { 00162 r1 = N_get_array_2d_d_value(weight_x, i, j); 00163 r2 = N_get_array_2d_d_value(weight_x, i + 1, j); 00164 mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */ 00165 } 00166 00167 res = mean * grad; 00168 00169 N_put_array_2d_d_value(field->x_array, i + 1, j, res); 00170 00171 } 00172 00173 for (j = 0; j < rows - 1; j++) 00174 for (i = 0; i < cols; i++) { 00175 grad = 0; 00176 mean = 0; 00177 00178 /* Only compute if the arrays are not null */ 00179 if (!N_is_array_2d_value_null(pot, i, j) && 00180 !N_is_array_2d_value_null(pot, i, j + 1)) { 00181 p1 = N_get_array_2d_d_value(pot, i, j); 00182 p2 = N_get_array_2d_d_value(pot, i, j + 1); 00183 grad = (p1 - p2) / dy; /* gradient */ 00184 } 00185 if (!N_is_array_2d_value_null(weight_y, i, j) && 00186 !N_is_array_2d_value_null(weight_y, i, j + 1)) { 00187 r1 = N_get_array_2d_d_value(weight_y, i, j); 00188 r2 = N_get_array_2d_d_value(weight_y, i, j + 1); 00189 mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */ 00190 } 00191 00192 res = -1 * mean * grad; 00193 00194 N_put_array_2d_d_value(field->y_array, i, j + 1, res); 00195 00196 } 00197 00198 /*Compute gradient field statistics */ 00199 N_calc_gradient_field_2d_stats(field); 00200 00201 return field; 00202 } 00203 00242 void 00243 N_compute_gradient_field_components_2d(N_gradient_field_2d * field, 00244 N_array_2d * x_comp, 00245 N_array_2d * y_comp) 00246 { 00247 int i, j; 00248 int rows, cols; 00249 double vx, vy; 00250 N_array_2d *x = x_comp; 00251 N_array_2d *y = y_comp; 00252 N_gradient_2d grad; 00253 00254 00255 if (!x) 00256 G_fatal_error("N_compute_gradient_components_2d: x array is empty"); 00257 if (!y) 00258 G_fatal_error("N_compute_gradient_components_2d: y array is empty"); 00259 00260 cols = field->x_array->cols; 00261 rows = field->x_array->rows; 00262 00263 /*Check the array sizes */ 00264 if (x->cols != cols || x->rows != rows) 00265 G_fatal_error 00266 ("N_compute_gradient_components_2d: the size of the x array doesn't fit the gradient field size"); 00267 if (y->cols != cols || y->rows != rows) 00268 G_fatal_error 00269 ("N_compute_gradient_components_2d: the size of the y array doesn't fit the gradient field size"); 00270 00271 for (j = 0; j < rows; j++) 00272 for (i = 0; i < cols; i++) { 00273 N_get_gradient_2d(field, &grad, i, j); 00274 00275 /* in case a gradient is zero, we expect a no flow boundary */ 00276 if (grad.WC == 0.0 || grad.EC == 0.0) 00277 vx = (grad.WC + grad.EC); 00278 else 00279 vx = (grad.WC + grad.EC) / 2; 00280 if (grad.NC == 0.0 || grad.SC == 0.0) 00281 vy = (grad.NC + grad.SC); 00282 else 00283 vy = (grad.NC + grad.SC) / 2; 00284 00285 N_put_array_2d_d_value(x, i, j, vx); 00286 N_put_array_2d_d_value(y, i, j, vy); 00287 } 00288 00289 return; 00290 } 00291 00300 void N_calc_gradient_field_3d_stats(N_gradient_field_3d * field) 00301 { 00302 double minx, miny, minz; 00303 double maxx, maxy, maxz; 00304 double sumx, sumy, sumz; 00305 int nonullx, nonully, nonullz; 00306 00307 G_debug(3, 00308 "N_calc_gradient_field_3d_stats: compute gradient field stats"); 00309 00310 N_calc_array_3d_stats(field->x_array, &minx, &maxx, &sumx, &nonullx, 0); 00311 N_calc_array_3d_stats(field->y_array, &miny, &maxy, &sumy, &nonully, 0); 00312 N_calc_array_3d_stats(field->z_array, &minz, &maxz, &sumz, &nonullz, 0); 00313 00314 if (minx <= minz && minx <= miny) 00315 field->min = minx; 00316 if (miny <= minz && miny <= minx) 00317 field->min = miny; 00318 if (minz <= minx && minz <= miny) 00319 field->min = minz; 00320 00321 if (maxx >= maxz && maxx >= maxy) 00322 field->max = maxx; 00323 if (maxy >= maxz && maxy >= maxx) 00324 field->max = maxy; 00325 if (maxz >= maxx && maxz >= maxy) 00326 field->max = maxz; 00327 00328 field->sum = sumx + sumy + sumz; 00329 field->nonull = nonullx + nonully + nonullz; 00330 field->mean = field->sum / (double)field->nonull; 00331 00332 return; 00333 } 00334 00335 00386 N_gradient_field_3d *N_compute_gradient_field_3d(N_array_3d * pot, 00387 N_array_3d * weight_x, 00388 N_array_3d * weight_y, 00389 N_array_3d * weight_z, 00390 N_geom_data * geom, 00391 N_gradient_field_3d * 00392 gradfield) 00393 { 00394 int i, j, k; 00395 int cols, rows, depths; 00396 double dx, dy, dz, p1, p2, r1, r2, mean, grad, res; 00397 N_gradient_field_3d *field = gradfield; 00398 00399 00400 if (pot->cols != weight_x->cols || pot->cols != weight_y->cols || 00401 pot->cols != weight_z->cols) 00402 G_fatal_error 00403 ("N_compute_gradient_field_3d: the arrays are not of equal size"); 00404 00405 if (pot->rows != weight_x->rows || pot->rows != weight_y->rows || 00406 pot->rows != weight_z->rows) 00407 G_fatal_error 00408 ("N_compute_gradient_field_3d: the arrays are not of equal size"); 00409 00410 if (pot->depths != weight_x->depths || pot->depths != weight_y->depths || 00411 pot->depths != weight_z->depths) 00412 G_fatal_error 00413 ("N_compute_gradient_field_3d: the arrays are not of equal size"); 00414 00415 if (pot->cols != geom->cols || pot->rows != geom->rows || 00416 pot->depths != geom->depths) 00417 G_fatal_error 00418 ("N_compute_gradient_field_3d: array sizes and geometry data are different"); 00419 00420 G_debug(3, "N_compute_gradient_field_3d: compute gradient field"); 00421 00422 cols = geom->cols; 00423 rows = geom->rows; 00424 depths = geom->depths; 00425 dx = geom->dx; 00426 dy = geom->dy; 00427 dz = geom->dz; 00428 00429 if (gradfield == NULL) { 00430 field = N_alloc_gradient_field_3d(cols, rows, depths); 00431 } 00432 else { 00433 if (field->cols != geom->cols || field->rows != geom->rows || 00434 field->depths != geom->depths) 00435 G_fatal_error 00436 ("N_compute_gradient_field_3d: gradient field sizes and geometry data are different"); 00437 } 00438 00439 for (k = 0; k < depths; k++) 00440 for (j = 0; j < rows; j++) 00441 for (i = 0; i < cols - 1; i++) { 00442 grad = 0; 00443 mean = 0; 00444 00445 /*Only compute if the arrays are not null */ 00446 if (!N_is_array_3d_value_null(pot, i, j, k) && 00447 !N_is_array_3d_value_null(pot, i + 1, j, k)) { 00448 p1 = N_get_array_3d_d_value(pot, i, j, k); 00449 p2 = N_get_array_3d_d_value(pot, i + 1, j, k); 00450 grad = (p1 - p2) / dx; /* gradient */ 00451 } 00452 if (!N_is_array_3d_value_null(weight_x, i, j, k) && 00453 !N_is_array_3d_value_null(weight_x, i + 1, j, k)) { 00454 r1 = N_get_array_3d_d_value(weight_x, i, j, k); 00455 r2 = N_get_array_3d_d_value(weight_x, i + 1, j, k); 00456 mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */ 00457 } 00458 00459 res = mean * grad; 00460 00461 G_debug(6, 00462 "N_compute_gradient_field_3d: X-direction insert value %6.5g at %i %i %i ", 00463 res, k, j, i + 1); 00464 00465 N_put_array_3d_d_value(field->x_array, i + 1, j, k, res); 00466 00467 } 00468 00469 for (k = 0; k < depths; k++) 00470 for (j = 0; j < rows - 1; j++) 00471 for (i = 0; i < cols; i++) { 00472 grad = 0; 00473 mean = 0; 00474 00475 /* Only compute if the arrays are not null */ 00476 if (!N_is_array_3d_value_null(pot, i, j, k) && 00477 !N_is_array_3d_value_null(pot, i, j + 1, k)) { 00478 p1 = N_get_array_3d_d_value(pot, i, j, k); 00479 p2 = N_get_array_3d_d_value(pot, i, j + 1, k); 00480 grad = (p1 - p2) / dy; /* gradient */ 00481 } 00482 if (!N_is_array_3d_value_null(weight_y, i, j, k) && 00483 !N_is_array_3d_value_null(weight_y, i, j + 1, k)) { 00484 r1 = N_get_array_3d_d_value(weight_y, i, j, k); 00485 r2 = N_get_array_3d_d_value(weight_y, i, j + 1, k); 00486 mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */ 00487 } 00488 00489 res = -1 * mean * grad; /*invert the direction, because we count from north to south, 00490 * but the gradient is defined in y direction */ 00491 00492 G_debug(6, 00493 "N_compute_gradient_field_3d: Y-direction insert value %6.5g at %i %i %i ", 00494 res, k, j + 1, i); 00495 00496 N_put_array_3d_d_value(field->y_array, i, j + 1, k, res); 00497 00498 } 00499 00500 for (k = 0; k < depths - 1; k++) 00501 for (j = 0; j < rows; j++) 00502 for (i = 0; i < cols; i++) { 00503 grad = 0; 00504 mean = 0; 00505 00506 /* Only compute if the arrays are not null */ 00507 if (!N_is_array_3d_value_null(pot, i, j, k) && 00508 !N_is_array_3d_value_null(pot, i, j, k + 1)) { 00509 p1 = N_get_array_3d_d_value(pot, i, j, k); 00510 p2 = N_get_array_3d_d_value(pot, i, j, k + 1); 00511 grad = (p1 - p2) / dz; /* gradient */ 00512 } 00513 if (!N_is_array_3d_value_null(weight_z, i, j, k) && 00514 !N_is_array_3d_value_null(weight_z, i, j, k + 1)) { 00515 r1 = N_get_array_3d_d_value(weight_z, i, j, k); 00516 r2 = N_get_array_3d_d_value(weight_z, i, j, k + 1); 00517 mean = N_calc_harmonic_mean(r1, r2); /*harmonical mean */ 00518 } 00519 00520 res = mean * grad; 00521 00522 G_debug(6, 00523 "N_compute_gradient_field_3d: Z-direction insert value %6.5g at %i %i %i ", 00524 res, k + 1, j, i); 00525 00526 N_put_array_3d_d_value(field->z_array, i, j, k + 1, res); 00527 00528 } 00529 00530 /*Compute gradient field statistics */ 00531 N_calc_gradient_field_3d_stats(field); 00532 00533 return field; 00534 } 00535 00578 void 00579 N_compute_gradient_field_components_3d(N_gradient_field_3d * field, 00580 N_array_3d * x_comp, 00581 N_array_3d * y_comp, 00582 N_array_3d * z_comp) 00583 { 00584 int i, j, k; 00585 int rows, cols, depths; 00586 double vx, vy, vz; 00587 N_array_3d *x = x_comp; 00588 N_array_3d *y = y_comp; 00589 N_array_3d *z = z_comp; 00590 N_gradient_3d grad; 00591 00592 00593 if (!x) 00594 G_fatal_error("N_compute_gradient_components_3d: x array is empty"); 00595 if (!y) 00596 G_fatal_error("N_compute_gradient_components_3d: y array is empty"); 00597 if (!z) 00598 G_fatal_error("N_compute_gradient_components_3d: z array is empty"); 00599 00600 cols = field->x_array->cols; 00601 rows = field->x_array->rows; 00602 depths = field->x_array->depths; 00603 00604 /*Check the array sizes */ 00605 if (x->cols != cols || x->rows != rows || x->depths != depths) 00606 G_fatal_error 00607 ("N_compute_gradient_components_3d: the size of the x array doesn't fit the gradient field size"); 00608 if (y->cols != cols || y->rows != rows || y->depths != depths) 00609 G_fatal_error 00610 ("N_compute_gradient_components_3d: the size of the y array doesn't fit the gradient field size"); 00611 if (z->cols != cols || z->rows != rows || z->depths != depths) 00612 G_fatal_error 00613 ("N_compute_gradient_components_3d: the size of the z array doesn't fit the gradient field size"); 00614 00615 for (k = 0; k < depths; k++) 00616 for (j = 0; j < rows; j++) 00617 for (i = 0; i < cols; i++) { 00618 N_get_gradient_3d(field, &grad, i, j, k); 00619 /* in case a gradient is zero, we expect a no flow boundary */ 00620 if (grad.WC == 0.0 || grad.EC == 0.0) 00621 vx = (grad.WC + grad.EC); 00622 else 00623 vx = (grad.WC + grad.EC) / 2; 00624 if (grad.NC == 0.0 || grad.SC == 0.0) 00625 vy = (grad.NC + grad.SC); 00626 else 00627 vy = (grad.NC + grad.SC) / 2; 00628 if (grad.TC == 0.0 || grad.BC == 0.0) 00629 vz = (grad.TC + grad.BC); 00630 else 00631 vz = (grad.TC + grad.BC) / 2; 00632 00633 N_put_array_3d_d_value(x, i, j, k, vx); 00634 N_put_array_3d_d_value(y, i, j, k, vy); 00635 N_put_array_3d_d_value(z, i, j, k, vz); 00636 } 00637 00638 00639 return; 00640 }