[GRASS-SVN] r51734 - in grass-addons/grass7/raster: . r.massvov
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 25 07:07:09 EDT 2012
Author: maxi
Date: 2012-05-25 04:07:09 -0700 (Fri, 25 May 2012)
New Revision: 51734
Added:
grass-addons/grass7/raster/r.massvov/
grass-addons/grass7/raster/r.massvov/Makefile
grass-addons/grass7/raster/r.massvov/filters.c
grass-addons/grass7/raster/r.massvov/filters.h
grass-addons/grass7/raster/r.massvov/general_f.c
grass-addons/grass7/raster/r.massvov/general_f.h
grass-addons/grass7/raster/r.massvov/main.c
grass-addons/grass7/raster/r.massvov/r.massmov.html
Log:
new addon: landslide runout simulation model
Added: grass-addons/grass7/raster/r.massvov/Makefile
===================================================================
--- grass-addons/grass7/raster/r.massvov/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.massvov/Makefile 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.massmov
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass-addons/grass7/raster/r.massvov/filters.c
===================================================================
--- grass-addons/grass7/raster/r.massvov/filters.c (rev 0)
+++ grass-addons/grass7/raster/r.massvov/filters.c 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,353 @@
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+
+#define nullo -999.9f
+
+float gradx3(float **matrix, int row, int col, float dx, int abs) {
+ float v;
+
+ if (matrix[row][col] != nullo && matrix[row][col + 1] != nullo
+ && matrix[row - 1][col + 1] != nullo
+ && matrix[row + 1][col + 1] != nullo
+ && matrix[row][col - 1] != nullo
+ && matrix[row - 1][col - 1] != nullo
+ && matrix[row + 1][col - 1] != nullo) {
+
+ if (abs == 1){
+ v = (fabs((matrix[row - 1][col + 1]) + 2 * fabs(matrix[row][col + 1])
+ + fabs(matrix[row + 1][col + 1]))
+ - (fabs(matrix[row - 1][col - 1]) + 2 * fabs(matrix[row][col - 1])
+ + fabs(matrix[row + 1][col - 1]))) / (8 * dx);
+ return v;
+ }
+ else{
+ v = ((matrix[row - 1][col + 1] + 2 * matrix[row][col + 1]
+ + matrix[row + 1][col + 1])
+ - (matrix[row - 1][col - 1] + 2 * matrix[row][col - 1]
+ + matrix[row + 1][col - 1])) / (8 * dx);
+ return v;
+ }
+ } else {
+ return nullo;
+ }
+}
+
+float grady3(float **matrix, int row, int col, float dy, int abs) {
+ float v;
+
+ if (matrix[row][col] != nullo && matrix[row - 1][col - 1] != nullo
+ && matrix[row - 1][col] != nullo
+ && matrix[row - 1][col + 1] != nullo
+ && matrix[row + 1][col - 1] != nullo
+ && matrix[row + 1][col] != nullo
+ && matrix[row + 1][col + 1] != nullo) {
+
+ if (abs == 1){
+ v = ((fabs(matrix[row - 1][col - 1]) + 2 * fabs(matrix[row - 1][col])
+ + fabs(matrix[row - 1][col + 1]))
+ - fabs((matrix[row + 1][col - 1]) + 2 * fabs(matrix[row + 1][col])
+ + fabs(matrix[row + 1][col + 1]))) / (8 * dy);
+ return v;
+ }
+ else {
+ v = ((matrix[row - 1][col - 1] + 2 * matrix[row - 1][col]
+ + matrix[row - 1][col + 1])
+ - (matrix[row + 1][col - 1] + 2 * matrix[row + 1][col]
+ + matrix[row + 1][col + 1])) / (8 * dy);
+ return v;
+ }
+
+ } else {
+ return nullo;
+ }
+}
+
+float gradx2(float **matrix, int row, int col, float dx, int abs) {
+ float v;
+
+ if (matrix[row][col] != nullo && matrix[row][col + 1] != nullo
+ && matrix[row][col - 1] != nullo) {
+
+ if (abs == 1){
+ v = (fabs(matrix[row][col + 1]) - fabs(matrix[row][col - 1])) / (2 * dx);
+ return v;
+ }
+ else {
+ v = (matrix[row][col + 1] - matrix[row][col - 1]) / (2 * dx);
+ return v;
+ }
+ } else {
+ return nullo;
+ }
+}
+
+float gradPx2(float **matrix1, float **matrix2, float **matrix3, int row,
+ int col, float dx) {
+ float v;
+
+ if (matrix1[row][col] != nullo && matrix2[row][col] != nullo
+ && (cos(atan(matrix3[row][col]))) != nullo
+ && matrix1[row][col + 1] != nullo && matrix2[row][col + 1] != nullo
+ && (cos(atan(matrix3[row][col + 1]))) != nullo
+ && matrix1[row][col - 1] != nullo && matrix2[row][col - 1] != nullo
+ && (cos(atan(matrix3[row][col - 1]))) != nullo) {
+ v = ((9.8 * (matrix1[row][col + 1] + matrix2[row][col + 1])
+ * (cos(atan(matrix3[row][col + 1]))))
+ - (9.8 * (matrix1[row][col - 1] + matrix2[row][col - 1])
+ * (cos(atan(matrix3[row][col - 1]))))) / (2 * dx);
+ return v;
+
+ } else
+ return nullo;
+}
+
+float grady2(float **matrix, int row, int col, float dy, int abs) {
+ float v;
+
+ if (matrix[row][col] != nullo && matrix[row + 1][col] != nullo
+ && matrix[row - 1][col] != nullo) {
+
+ if (abs == 1){
+ v = (fabs(matrix[row - 1][col]) - fabs(matrix[row + 1][col])) / (2 * dy);
+ return v;
+ }
+
+ else {
+ v = (matrix[row - 1][col] - matrix[row + 1][col]) / (2 * dy);
+ return v;
+ }
+ } else {
+ return nullo;
+ }
+}
+
+/* calcolo del gradiente combinato della somma di 2 matrici (usato per P)
+ * gradPy2 (pendenza y, matrice 1, matrice 2, riga, col, res y)
+ *
+ * */
+float gradPy2(float **matrix1, float **matrix2, float **matrix3, int row,
+ int col, float dy) {
+ float v;
+
+ if (matrix1[row][col] != nullo && matrix2[row][col] != nullo
+ && (cos(atan(matrix3[row][col]))) != nullo
+ && matrix1[row + 1][col] != nullo && matrix2[row + 1][col] != nullo
+ && (cos(atan(matrix3[row + 1][col]))) != nullo
+ && matrix1[row - 1][col] != nullo && matrix2[row - 1][col] != nullo
+ && (cos(atan(matrix3[row - 1][col]))) != nullo) {
+ v = ((9.8 * (matrix1[row - 1][col] + matrix2[row - 1][col])
+ * (cos(atan(matrix3[row - 1][col]))))
+ - (9.8 * (matrix1[row + 1][col] + matrix2[row + 1][col])
+ * (cos(atan(matrix3[row + 1][col]))))) / (2 * dy);
+ return v;
+
+ } else
+ return nullo;
+}
+
+float lax(float **matrix, int row, int col, float laxfactor) {
+
+ float gg = 0.0;
+ float hh = 0.0;
+ float v;
+
+ if (matrix[row][col] != nullo) {
+
+ if (matrix[row - 1][col - 1] != nullo) {
+ gg = gg + 2 * matrix[row - 1][col - 1];
+ hh = hh + 2;
+ }
+
+
+ if (matrix[row - 1][col] != nullo) {
+ gg = gg + 3 * matrix[row - 1][col];
+ hh = hh + 3;
+ }
+
+ if (matrix[row - 1][col + 1] != nullo) {
+ gg = gg + 2 * matrix[row - 1][col + 1];
+ hh = hh + 2;
+ }
+
+ if (matrix[row][col - 1] != nullo) {
+ gg = gg + 3 * matrix[row][col - 1];
+ hh = hh + 3;
+ }
+
+ if (matrix[row][col + 1] != nullo) {
+ gg = gg + 3 * matrix[row][col + 1];
+ hh = hh + 3;
+ }
+
+ if (matrix[row + 1][col - 1] != nullo) {
+ gg = gg + 2 * matrix[row + 1][col - 1];
+ hh = hh + 2;
+ }
+
+ if (matrix[row + 1][col] != nullo) {
+ gg = gg + 3 * matrix[row + 1][col];
+ hh = hh + 3;
+ }
+
+ if (matrix[row + 1][col + 1] != nullo) {
+ gg = gg + 2 * matrix[row + 1][col + 1];
+ hh = hh + 2;
+ }
+
+ if (/*gg != 0.0 &&*/ hh != 0.0)
+ v = ((1 - laxfactor) * matrix[row][col] + laxfactor * (gg / hh));
+ else
+ v = matrix[row][col];
+
+ return v;
+ }
+
+ else {
+ return nullo;
+ }
+
+}
+
+float filter_lax(float **matrix, int row, int col, float laxfactor, float **filter_matrix,float threshold,float val) {
+
+
+ float gg = 0.0;
+ float hh = 0.0;
+ float v;
+
+
+ if (matrix[row][col] != nullo && (filter_matrix[row][col] > threshold)) {
+
+ if ((matrix[row - 1][col - 1] != nullo)&&(filter_matrix[row - 1][col - 1] > threshold)) {
+ gg = gg + 2 * matrix[row - 1][col - 1];
+ hh = hh + 2;
+ }
+
+
+ if ((matrix[row - 1][col] != nullo)&&(filter_matrix[row - 1][col] > threshold)) {
+ gg = gg + 3 * matrix[row - 1][col];
+ hh = hh + 3;
+ }
+
+ if ((matrix[row - 1][col + 1] != nullo)&&(filter_matrix[row - 1][col + 1] > threshold)) {
+ gg = gg + 2 * matrix[row - 1][col + 1];
+ hh = hh + 2;
+ }
+
+ if ((matrix[row][col - 1] != nullo)&&(filter_matrix[row][col - 1] > threshold)) {
+ gg = gg + 3 * matrix[row][col - 1];
+ hh = hh + 3;
+ }
+
+ if ((matrix[row][col + 1] != nullo)&&(filter_matrix[row][col + 1] > threshold)) {
+ gg = gg + 3 * matrix[row][col + 1];
+ hh = hh + 3;
+ }
+
+ if ((matrix[row + 1][col - 1] != nullo)&&(filter_matrix[row + 1][col - 1] > threshold)) {
+ gg = gg + 2 * matrix[row + 1][col - 1];
+ hh = hh + 2;
+ }
+
+ if ((matrix[row + 1][col] != nullo)&&(filter_matrix[row + 1][col] > threshold)) {
+ gg = gg + 3 * matrix[row + 1][col];
+ hh = hh + 3;
+ }
+
+ if ((matrix[row + 1][col + 1] != nullo)&&(filter_matrix[row + 1][col + 1] > threshold)) {
+ gg = gg + 2 * matrix[row + 1][col + 1];
+ hh = hh + 2;
+ }
+
+ if (/*gg != 0.0 &&*/ hh != 0.0)
+ v = ((1 - laxfactor) * matrix[row][col] + laxfactor * (gg / hh));
+ else
+ v = matrix[row][col];
+
+ return v;
+ }
+
+ else {
+ v=val;
+ return v;
+ }
+
+}
+
+
+float filter_lax_print(float **matrix, int row, int col, float laxfactor, float **filter_matrix,float threshold,float val) {
+
+
+ float gg = 0.0;
+ float hh = 0.0;
+ float v;
+
+ if (matrix[row][col] != nullo && (filter_matrix[row][col] > threshold)) {
+ if ((matrix[row - 1][col - 1] != nullo)&&(filter_matrix[row - 1][col - 1] > threshold)) {
+ gg = gg + 2 * matrix[row - 1][col - 1];
+ hh = hh + 2;
+ while (getchar() != 'y') {}
+ }
+
+
+ if ((matrix[row - 1][col] != nullo)&&(filter_matrix[row - 1][col] > threshold)) {
+ gg = gg + 3 * matrix[row - 1][col];
+ hh = hh + 3;
+ while (getchar() != 'y') {}
+ }
+
+ if ((matrix[row - 1][col + 1] != nullo)&&(filter_matrix[row - 1][col + 1] > threshold)) {
+ gg = gg + 2 * matrix[row - 1][col + 1];
+ hh = hh + 2;
+ while (getchar() != 'y') {}
+ }
+
+ if ((matrix[row][col - 1] != nullo)&&(filter_matrix[row][col - 1] > threshold)) {
+ gg = gg + 3 * matrix[row][col - 1];
+ hh = hh + 3;
+ while (getchar() != 'y') {}
+ }
+
+ if ((matrix[row][col + 1] != nullo)&&(filter_matrix[row][col + 1] > threshold)) {
+ gg = gg + 3 * matrix[row][col + 1];
+ hh = hh + 3;
+ while (getchar() != 'y') {}
+ }
+
+ if ((matrix[row + 1][col - 1] != nullo)&&(filter_matrix[row + 1][col - 1] > threshold)) {
+ gg = gg + 2 * matrix[row + 1][col - 1];
+ hh = hh + 2;
+ while (getchar() != 'y') {}
+ }
+
+ if ((matrix[row + 1][col] != nullo)&&(filter_matrix[row + 1][col] > threshold)) {
+ gg = gg + 3 * matrix[row + 1][col];
+ hh = hh + 3;
+ while (getchar() != 'y') {}
+ }
+
+ if ((matrix[row + 1][col + 1] != nullo)&&(filter_matrix[row + 1][col + 1] > threshold)) {
+ gg = gg + 2 * matrix[row + 1][col + 1];
+ hh = hh + 2;
+ while (getchar() != 'y') {}
+ }
+
+ if (/*gg != 0.0 &&*/ hh != 0.0)
+ v = ((1 - laxfactor) * matrix[row][col] + laxfactor * (gg / hh));
+ else
+ v = matrix[row][col];
+
+ return v;
+ }
+
+ else {
+ v=val;
+ return v;
+ }
+}
Added: grass-addons/grass7/raster/r.massvov/filters.h
===================================================================
--- grass-addons/grass7/raster/r.massvov/filters.h (rev 0)
+++ grass-addons/grass7/raster/r.massvov/filters.h 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,12 @@
+float gradx3(float **matrix, int row, int col, float dx, int abs);
+float grady3(float **matrix, int row, int col, float dy, int abs);
+float gradx2(float **matrix, int row, int col, float dx, int abs);
+float grady2(float **matrix, int row, int col, float dy, int abs);
+float lax(float **matrix, int row, int col, float laxfactor);
+float filter_lax(float **matrix, int row, int col, float laxfactor, float **filter_matrix,float threshold,float val);
+float gradPx2(float **matrix1, float **matrix2, float **matrix3, int row,
+ int col, float dx);
+float gradPy2(float **matrix1, float **matrix2, float **matrix3, int row,
+ int col, float dy);
+float isocell(float **matrix,int row,int col,float val);
+float filter_lax_print(float **matrix, int row, int col, float laxfactor, float **filter_matrix,float threshold,float val);
Added: grass-addons/grass7/raster/r.massvov/general_f.c
===================================================================
--- grass-addons/grass7/raster/r.massvov/general_f.c (rev 0)
+++ grass-addons/grass7/raster/r.massvov/general_f.c 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,261 @@
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include <grass/linkm.h>
+
+#define HBExp 1
+#define verysmall 0.01
+#define min(A,B) ((A) < (B) ? (A):(B))
+#define max(A,B) ((A) > (B) ? (A):(B))
+
+float **G_alloc_fmatrix(int rows, int cols) {
+ float **m;
+ int i;
+ m = (float **) G_calloc(rows, sizeof(float *));
+ m[0] = (float *) G_calloc(rows * cols, sizeof(float));
+ for (i = 1; i < rows; i++)
+ m[i] = m[i - 1] + cols;
+ return m;
+}
+
+int **G_alloc_imatrix(int rows, int cols) {
+ int **mmm;
+ int i;
+ mmm = (int **) G_calloc(rows, sizeof(int *));
+ mmm[0] = (int *) G_calloc(rows * cols, sizeof(int));
+ for (i = 1; i < rows; i++)
+ mmm[i] = mmm[i - 1] + cols;
+ return mmm;
+}
+
+void G_free_fmatrix(float **m) {
+ G_free(m[0]);
+ G_free(m);
+ m = NULL;
+ return;
+}
+
+void G_free_imatrix(int **mmm) {
+ G_free(mmm[0]);
+ G_free(mmm);
+ mmm = NULL;
+ return;
+}
+
+int check_rheol_par(int rheol_type, float chezy, float visco, float rho) {
+ if (rheol_type == 2) {
+ if (chezy > 0)
+ return 1;
+ else
+ return -2;
+ }
+
+ if (rheol_type == 3) {
+ if (visco > 0 && rho > 0)
+ return 1;
+ else
+ return -3;
+ }
+}
+
+float t_frict(float **h, int row, int col, float b_frict) {
+ float t;
+
+ if (h[row][col] > verysmall) {
+ t = tan((M_PI * b_frict) / 180.0);
+ } else {
+ t = 5 * tan((M_PI * b_frict) / 180.0);
+ }
+ return t;
+}
+
+float t_voellmy(float v, float **h, int row, int col, float b_frict,
+ float chezy) {
+ float t;
+ if (h[row][col] > verysmall) {
+ t = tan((M_PI * b_frict) / 180.0) + pow(v, 2) / (chezy * h[row][col]);
+ } else {
+ t = 5 * tan((M_PI * b_frict) / 180.0);
+ }
+ return t;
+}
+
+float t_visco(float v, float **h, int row, int col, float b_frict, float rho,
+ float visco, float ystress) {
+ float t;
+
+ if (h[row][col] > verysmall) {
+ if (ystress > 0) {
+ t =
+ tan((M_PI * b_frict) / 180.0)
+ + (1.5 * ystress
+ + (3 * visco * pow(v, HBExp) / h[row][col]))
+ / (rho * h[row][col]);
+ } else {
+ t = tan((M_PI * b_frict) / 180.0)
+ + (3 * visco * pow(v, HBExp) / h[row][col] / h[row][col]);
+ }
+ } else {
+ if (ystress > 0) {
+ t = 5 * ystress;
+ } else {
+ t = 5 * tan((M_PI * b_frict) / 180.0); ///// CHECK settare a 0 ???
+ }
+ }
+ return t;
+}
+
+
+float veldt(float ua, float t, float g_x, float p_x, float i_x, float t_x) {
+ float v;
+
+ if (ua>0)
+ v = max(0, ua + t * (g_x + p_x + i_x - t_x));
+
+ else if (ua<0)
+ v = min(ua + t * (g_x + p_x - i_x + t_x),0);
+
+ else {
+ if ((g_x+p_x)>0)
+ v = max(0, t * (g_x + p_x + i_x - t_x));
+
+ else if ((g_x+p_x)<0){
+ v = min(t * (g_x + p_x - i_x + t_x),0);
+ }
+ else
+ v=0;
+ }
+
+ return v;
+}
+
+float shift0(float **m, int r, int c, int maxR, int maxC, int minR, int minC, int n, int w) {
+ float v;
+
+ if ((r+n<minR) || (r+n>maxR))
+ v=0;
+ else if ((c+w<minC)|| (c+w>maxC))
+ v=0;
+ else
+ v = m[r+n][c+w];
+
+ return v;
+}
+
+
+void out_print (float **matrix,char *name, int nr, int nc){
+ int row,col;
+ float *outrast;
+ int outfd;
+
+ outrast = Rast_allocate_f_buf();
+ outfd = Rast_open_fp_new(name);
+ for (row = 0; row < nr; row++) {
+ for (col = 0; col < nc; col++) {
+ if (matrix[row][col]==0)
+ Rast_set_f_null_value(&outrast[col], 1);
+ else
+ ((FCELL *) outrast)[col] = matrix[row][col];
+ }
+ Rast_put_f_row(outfd, outrast);
+ }
+ G_free(outrast);
+ Rast_close(outfd);
+}
+
+
+void out_sum_print (float **matrix1, float **matrix2, float **matrix3, float **matrix4,char *name, int nr, int nc, int mode, float threshold){
+ int row,col;
+ float *outrast;
+ int outfd;
+ /*struct History history; * holds meta-data */
+
+ outrast = Rast_allocate_f_buf();
+ outfd = Rast_open_fp_new(name);
+ for (row = 0; row < nr; row++) {
+ for (col = 0; col < nc; col++) {
+ if (mode==1){
+ if (matrix1[row][col]+matrix2[row][col]>0)
+ ((FCELL *) outrast)[col] = matrix1[row][col]+matrix2[row][col];
+ else
+ Rast_set_f_null_value(&outrast[col], 1);
+ }
+ if (mode==2){
+ if (matrix1[row][col]+matrix2[row][col]>0)
+ if (matrix2[row][col]>threshold)
+ ((FCELL *) outrast)[col] = sqrt(pow(matrix3[row][col],2)+pow(matrix4[row][col],2));
+ else
+ ((FCELL *) outrast)[col] = 0.0;
+ else
+ Rast_set_f_null_value(&outrast[col], 1);
+ }
+ }
+ Rast_put_f_row(outfd, outrast);
+ /**
+ Rast_short_history(name, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(name, &history);
+ **/
+ }
+ G_free(outrast);
+ Rast_close(outfd);
+}
+
+float nash_sutcliffe(float **m_t1,float**m_t2, int nr, int nc){
+ float sum_den=0;
+ float sum_num=0;
+ float sum_ave=0;
+ int row,col,c=0;
+ float ns_res,ave;
+
+ for (row = 0; row < nr; row++) {
+ for (col = 0; col < nc; col++) {
+ sum_ave +=m_t1[row][col];
+ c+=1;
+ }
+ }
+
+ ave = sum_ave/c;
+
+ for (row = 0; row < nr; row++) {
+ for (col = 0; col < nc; col++) {
+ sum_num += pow(m_t1[row][col]-m_t2[row][col],2);
+ sum_den += pow(m_t1[row][col]-ave,2);
+ }
+ }
+ if (sum_den==0){
+ ns_res = -9999999;
+ G_message("WARNING: -inf value obtained for NS coefficient; the new value -9999999 has been set");
+ } else {
+ ns_res = 1 - (sum_num/sum_den);
+ }
+ return ns_res;
+}
+
+
+void report_input(float ifrict,float rho,float ystress,float visco,float chezy,float bfrict,float fluid,float NS_thres,int t,int delta){
+ fprintf(stdout, "-----------Input data:-----------\n");
+ fprintf(stdout, "Internal friction angle = %.2f\n",ifrict);
+ if (rho!=-1)
+ fprintf(stdout, "Density = %.2f\n",rho);
+ if (ystress!=-1)
+ fprintf(stdout, "Yield stress = %.2f\n",ystress);
+ if (visco!=-1)
+ fprintf(stdout, "Viscosity = %.2f\n",visco);
+ if (chezy!=-1)
+ fprintf(stdout, "Chézy coefficient = %.2f\n",chezy);
+ if (bfrict!=-1)
+ fprintf(stdout, "Basal friction angle = %.2f\n",bfrict);
+ if (fluid!=-1)
+ fprintf(stdout, "Fluidization rate = %.2f\n",fluid);
+ if (NS_thres!=-1)
+ fprintf(stdout, "Nash-Sutcliffe threshold = %.4f\n",NS_thres);
+ fprintf(stdout, "Maximum timesteps number = %i\n",t);
+ if (delta!=-1)
+ fprintf(stdout, "Reporting time frequency = %i\n",delta);
+ fprintf(stdout, "---------------------------------\n");
+}
Added: grass-addons/grass7/raster/r.massvov/general_f.h
===================================================================
--- grass-addons/grass7/raster/r.massvov/general_f.h (rev 0)
+++ grass-addons/grass7/raster/r.massvov/general_f.h 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,16 @@
+float **G_alloc_fmatrix(int rows, int cols);
+int **G_alloc_imatrix(int rows, int cols);
+void G_free_fmatrix(float **m);
+void G_free_imatrix(int **mmm);
+int check_rheol_par(int rheol_type, float chezy, float visco, float rho);
+float t_frict(float **h, int row, int col, float b_frict);
+float t_voellmy(float vel, float **h, int row, int col, float b_frict,
+ float chezy);
+float t_visco(float v, float **h, int row, int col, float b_frict, float rho,
+ float visco, float ystress);
+float veldt(float ua, float t, float g_x, float p_x, float i_x, float t_x);
+float shift0(float **m, int r, int c, int maxR, int maxC, int minR, int minC, int n, int w);
+void out_print (float **matrix,char *name, int nr, int nc);
+void out_sum_print (float **matrix1,float **matrix2,float **matrix3,float **matrix4, char *name, int nr, int nc, int mode, float threshold);
+float nash_sutcliffe(float **m_t1,float**m_t2, int nr, int nc);
+void report_input(float ifrict,float rho,float ystress,float visco,float chezy,float bfrict,float fluid,float NS_thres,int t,int delta);
Added: grass-addons/grass7/raster/r.massvov/main.c
===================================================================
--- grass-addons/grass7/raster/r.massvov/main.c (rev 0)
+++ grass-addons/grass7/raster/r.massvov/main.c 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,1122 @@
+/*****************************************************************************
+ *
+ * MODULE: r.massmov
+ *
+ * AUTHOR(S): Original author:
+ * Santiago Begueria (EEAD-CSIC)
+ * Monia Molinari <monia.molinari supsi.ch>
+ * Massimiliano Cannata <massimiliano.cannata supsi.ch>
+ *
+ *
+ * PURPOSE: This module is intended to provide the capability of simulating
+ * mass movement (fast landslide) over complex topography. This module
+ * applies the Shallow Water Equation (SWE) with different types of
+ * rheologies (frictional,Voellmy,viscoplastic)
+ *
+ *
+ * COPYRIGHT: (C) 2012 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * Licence (>=2). Read the file COPYING that cames with GRASS
+ * for details.
+ *
+ ****************************************************************************/
+
+/* libraries*/
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "filters.h"
+#include "general_f.h"
+
+/* numerical stability control */
+#define CFLlimsup 0.5 /* Higher value of the Courant-Friedrichs-Levy */
+#define CFLliminf 0.3 /* Lower value of the Courant-Friedrichs-Levy */
+#define MinNLoops 1 /* Minimum number of internal loops */
+#define MaxNLoops 124 /* Maximum number of internal loops */
+#define InitialLoops 1 /* Initial number of internal loops */
+#define laxfactor 0.5 /* Low pass filtering for the lax function (0 = no filtering) */
+#define dTLeap 0.3 /* Fraction of dT, used for time extrapolation of fluxes */
+#define verysmall 0.01 /* threshold to avoid div by very small values */
+#define nullo -999.9f
+
+/* timestep control */
+#define dT 1.0 /* Timeslice */
+
+/* other constants */
+#define grav 9.8 /* Gravity acceleration */
+
+/* functions definition */
+#define min(A,B) ((A) < (B) ? (A):(B))
+#define max(A,B) ((A) > (B) ? (A):(B))
+
+/*
+ * global functions declaration
+ */
+
+/*
+ * main function
+ *
+ *
+ */
+
+int main(int argc, char *argv[]) {
+ struct Cell_head cellhd; /* it stores region information and header information of rasters */
+ struct Cell_head window;
+ struct History history; /* holds meta-data */
+
+ /* mapset name locator */
+ char *mapset_ELEV;
+ char *mapset_HINI;
+ char *mapset_DIST;
+
+ /* input buffer */
+ FCELL *inrast_ELEV;
+ FCELL *inrast_HINI;
+ FCELL *inrast_DIST;
+
+ /* output buffer */
+ float *outrast_H;
+ float *outrast_V;
+
+ /* input file descriptor */
+ int infd_ELEV;
+ int infd_HINI;
+ int infd_DIST;
+
+ /* output file descriptor */
+ int outfd_H;
+ int outfd_V;
+
+ /* cell counters */
+ int row, col;
+ int nrows, ncols;
+ float res_ns, res_ew, ca;
+
+ /* input & output arguments */
+ int TIMESTEPS, DELTAT, RHEOL_TYPE;
+ float RHO, YSTRESS, CHEZY, VISCO, BFRICT, IFRICT, FLUID,NS_THRES;
+ char RHEOL;
+ char name1[20];
+ char *name_ELEV;
+ char *name_HINI;
+ char *name_DIST;
+ char *result_H;
+ char *result_V;
+ char *result_HMAX;
+ char *result_VMAX;
+
+ /* memory matrix */
+ float **m_ELEV, **m_HINI, **m_DIST;
+ int **m_OUTLET;
+ float **m_gx, **m_gy, **m_slope;
+ float **m_H, **m_Hloop_dt, **m_Hloop, **m_V, **m_Vloop, **m_Vloop_dt, **m_U,
+ **m_Uloop, **m_Uloop_dt, **m_HUloop, **m_HVloop, **m_Hold, **m_Hmax, **m_Vmax;
+ float **m_K, **m_Kloop, **m_Hland, **TEST;
+
+ /*int variables */
+
+ float elev_fcell_val, dist_fcell_val, hini_fcell_val;
+ int outlet_cell_val;
+
+ float gx, gy, slope;
+ float k_act, k_pas;
+ float G_x, G_y, I_x, I_y, P_x, P_y, T, T_x, T_y, T_x_b, T_y_b, T_b;
+ float ddt, dt;
+ float Uloop_a, Vloop_a, Uloop_b, Vloop_b, vel, vel_b, Uloop_dt, Vloop_dt, Hloop_a, Hloop_dt, dH_dT, H_fluid;
+ float CFL, CFL_u, CFL_v;
+ float Vol_gone_loop, Vol_exp_loop, mbe;
+ float NS,NSold;
+ int testPar;
+
+ /***********************************************************************************************************************/
+ /* GRASS structure */
+ struct GModule *module;
+ struct Option *input_ELEV, *input_HINI,*input_DIST,
+ *input_RHEOL, *input_RHO, *input_YSTRESS, *input_CHEZY,
+ *input_VISCO, *input_BFRICT, *input_IFRICT, *input_FLUID,
+ *input_TIMESTEPS, *input_DELTAT, *output_H, *output_VEL, *input_NS_THRES, *output_HMAX, *output_VMAX;
+ struct Flag *flag_i;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]);
+
+ /* initialize module */
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("landslide"));
+ G_add_keyword(_("model"));
+ module->description =_("Estimation of run-out and deposition of landslide phenomena over a complex topography");
+
+ /* define the different options */
+
+ input_ELEV = G_define_option();
+ input_ELEV->key = "elev";
+ input_ELEV->type = TYPE_STRING;
+ input_ELEV->required = YES;
+ input_ELEV->gisprompt = "old,cell,raster";
+ input_ELEV->description = _("Name of elevation raster map");
+ input_ELEV->guisection = _("Input options");
+
+ input_HINI = G_define_option();
+ input_HINI->key = "h_ini";
+ input_HINI->type = TYPE_STRING;
+ input_HINI->required = YES;
+ input_HINI->gisprompt = "old,cell,raster";
+ input_HINI->description = _("Name of landslide initial body thickness raster map");
+ input_HINI->guisection = _("Input options");
+
+ input_DIST = G_define_option();
+ input_DIST->key = "fluiddist";
+ input_DIST->type = TYPE_STRING;
+ input_DIST->required = YES;
+ input_DIST->gisprompt = "old,cell,raster";
+ input_DIST->description = _("Name of distance from the landlide toe raster map");
+ input_DIST->guisection = _("Input options");
+
+ input_RHEOL = G_define_option();
+ input_RHEOL->key = "rheology";
+ input_RHEOL->type = TYPE_STRING;
+ input_RHEOL->required = YES;
+ input_RHEOL->options = "frictional,Voellmy,viscoplastic";
+ input_RHEOL->description = _("Name of rheological law");
+ input_RHEOL->guisection = _("Input options");
+
+ input_RHO = G_define_option();
+ input_RHO->key = "rho";
+ input_RHO->type = TYPE_DOUBLE;
+ input_RHO->required = NO;
+ input_RHO->multiple = NO;
+ input_RHO->description =_("Density of the flow [Kg/m3]. Required only for viscous rheologies.");
+ input_RHO->guisection = _("Input options");
+
+ input_YSTRESS = G_define_option();
+ input_YSTRESS->key = "ystress";
+ input_YSTRESS->type = TYPE_DOUBLE;
+ input_YSTRESS->required = NO;
+ input_YSTRESS->multiple = NO;
+ input_YSTRESS->description =_("Apparent yield stress [Pa]. Used only for viscous rheologies (optional).");
+ input_YSTRESS->guisection = _("Input options");
+
+ input_VISCO = G_define_option();
+ input_VISCO->key = "visco";
+ input_VISCO->type = TYPE_DOUBLE;
+ input_VISCO->required = NO;
+ input_VISCO->multiple = NO;
+ input_VISCO->description = _("Dynamic viscosity [Pa*s]. Required only for viscous rheologies");
+ input_VISCO->guisection = _("Input options");
+
+ input_CHEZY = G_define_option();
+ input_CHEZY->key = "chezy";
+ input_CHEZY->type = TYPE_DOUBLE;
+ input_CHEZY->required = NO;
+ input_CHEZY->multiple = NO;
+ input_CHEZY->description =_("Chezy roughness coefficient [m/s2]. Required only for Voellmy rheology");
+ input_CHEZY->guisection = _("Input options");
+
+ input_BFRICT = G_define_option();
+ input_BFRICT->key = "bfrict";
+ input_BFRICT->type = TYPE_DOUBLE;
+ input_BFRICT->required = NO;
+ input_BFRICT->multiple = NO;
+ input_BFRICT->description = _("Angle of basal friction [deg]");
+ input_BFRICT->guisection = _("Input options");
+
+ input_IFRICT = G_define_option();
+ input_IFRICT->key = "ifrict";
+ input_IFRICT->type = TYPE_DOUBLE;
+ input_IFRICT->required = YES;
+ input_IFRICT->multiple = NO;
+ input_IFRICT->description = _("Angle of internal friction [deg]");
+ input_IFRICT->guisection = _("Input options");
+
+ input_FLUID = G_define_option();
+ input_FLUID->key = "fluid";
+ input_FLUID->type = TYPE_DOUBLE;
+ input_FLUID->required = YES;
+ input_FLUID->multiple = NO;
+ input_FLUID->description = _("Upward velocity of transition from solid to fluid of the landsliding mass [m/s]");
+ input_FLUID->guisection = _("Input options");
+
+ input_TIMESTEPS = G_define_option();
+ input_TIMESTEPS->key = "timesteps";
+ input_TIMESTEPS->type = TYPE_INTEGER;
+ input_TIMESTEPS->required = YES;
+ input_TIMESTEPS->multiple = NO;
+ input_TIMESTEPS->description = _(
+ "Maximum number of time steps of the simulation [s]");
+ input_TIMESTEPS->guisection = _("Input options");
+
+ input_DELTAT = G_define_option();
+ input_DELTAT->key = "deltatime";
+ input_DELTAT->type = TYPE_INTEGER;
+ input_DELTAT->required = NO;
+ input_DELTAT->multiple = NO;
+ input_DELTAT->description = _("Reporting time frequency [s]");
+ input_DELTAT->guisection = _("Input options");
+
+ input_NS_THRES = G_define_option();
+ input_NS_THRES->key = "ns_thres";
+ input_NS_THRES->type = TYPE_DOUBLE;
+ input_NS_THRES->required = NO;
+ input_NS_THRES->multiple = NO;
+ input_NS_THRES->description = _("Threshold for Nash-Sutcliffe convergence [-]");
+ input_NS_THRES->guisection = _("Input options");
+
+ output_H = G_define_option();
+ output_H->key = "h";
+ output_H->type = TYPE_STRING;
+ output_H->required = NO;
+ output_H->gisprompt = "new,cell,raster";
+ output_H->description = _("Prefix for flow thickness output raster maps");
+ output_H->guisection = _("Output options");
+
+ output_HMAX = G_define_option();
+ output_HMAX->key = "h_max";
+ output_HMAX->type = TYPE_STRING;
+ output_HMAX->required = NO;
+ output_HMAX->gisprompt = "new,cell,raster";
+ output_HMAX->description = _("Prefix for maximum flow thickness output raster maps");
+ output_HMAX->guisection = _("Output options");
+
+ output_VEL = G_define_option();
+ output_VEL->key = "v";
+ output_VEL->type = TYPE_STRING;
+ output_VEL->required = NO;
+ output_VEL->gisprompt = "new,cell,raster";
+ output_VEL->description = _("Prefix for flow velocity output raster maps");
+ output_VEL->guisection = _("Output options");
+
+ output_VMAX = G_define_option();
+ output_VMAX->key = "v_max";
+ output_VMAX->type = TYPE_STRING;
+ output_VMAX->required = NO;
+ output_VMAX->gisprompt = "new,cell,raster";
+ output_VMAX->description = _("Prefix for maximum flow velocity output raster maps");
+ output_VMAX->guisection = _("Output options");
+
+ flag_i = G_define_flag();
+ flag_i->key = 'i';
+ flag_i->description = _("Print input data");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /***********************************************************************************************************************/
+ /* get entered parameters */
+
+ name_ELEV = input_ELEV->answer;
+ name_DIST = input_DIST->answer;
+ name_HINI = input_HINI->answer;
+
+ result_H = output_H->answer;
+ result_V = output_VEL->answer;
+ result_HMAX = output_HMAX->answer;
+ result_VMAX = output_VMAX->answer;
+
+ G_debug(1,"Getting numeric input data...");
+
+ if (strcmp(input_RHEOL->answer, "frictional") == 0)
+ RHEOL_TYPE = 1;
+ if (strcmp(input_RHEOL->answer, "Voellmy") == 0)
+ RHEOL_TYPE = 2;
+ if (strcmp(input_RHEOL->answer, "viscoplastic") == 0)
+ RHEOL_TYPE = 3;
+
+
+ if (input_RHO->answer)
+ sscanf(input_RHO->answer, "%f", &RHO);
+ else
+ RHO = -1;
+
+ if (input_YSTRESS->answer)
+ sscanf(input_YSTRESS->answer, "%f", &YSTRESS);
+ else
+ YSTRESS = -1;
+
+ if (input_CHEZY->answer)
+ sscanf(input_CHEZY->answer, "%f", &CHEZY);
+ else
+ CHEZY = -1;
+
+ if (input_VISCO->answer)
+ sscanf(input_VISCO->answer, "%f", &VISCO);
+ else
+ VISCO = -1;
+
+ if (input_BFRICT->answer)
+ sscanf(input_BFRICT->answer, "%f", &BFRICT);
+ else
+ BFRICT = -1;
+
+ if (input_NS_THRES->answer)
+ sscanf(input_NS_THRES->answer, "%f", &NS_THRES);
+ else
+ NS_THRES = -1;
+
+ sscanf(input_IFRICT->answer, "%f", &IFRICT);
+
+ sscanf(input_FLUID->answer, "%f", &FLUID);
+
+ if (input_DELTAT->answer)
+ DELTAT = atoi(input_DELTAT->answer);
+ else
+ DELTAT = -1;
+
+ TIMESTEPS = atoi(input_TIMESTEPS->answer);
+
+ /* check outputs required */
+ if (!((result_H) || (result_V) || (result_HMAX) || (result_VMAX))){
+ G_fatal_error(_("No output specified"));
+ }
+
+ /* check rheology parameters */
+ testPar = check_rheol_par(RHEOL_TYPE, CHEZY, VISCO, RHO);
+
+ if (testPar == -2)
+ G_fatal_error(_("For the selected rheology Chezy parameter is required"));
+
+ if (testPar == -3)
+ G_fatal_error(_("For the selected rheology viscosity and density parameters are required"));
+
+ /* report Input Data */
+ if (flag_i->answer){
+ report_input(IFRICT,RHO, YSTRESS, VISCO, CHEZY, BFRICT,FLUID, NS_THRES, TIMESTEPS, DELTAT);
+ }
+
+ G_debug(2,"Verifiyng raster maps input data...");
+
+ /* find maps in mapset */
+ mapset_ELEV = (char *) G_find_raster2(name_ELEV, "");
+ if (mapset_ELEV == NULL)
+ G_fatal_error(_("cell file [%s] not found"), name_ELEV);
+
+ mapset_DIST = (char *) G_find_raster2(name_DIST, "");
+ if (mapset_DIST == NULL)
+ G_fatal_error(_("cell file [%s] not found"), name_DIST);
+
+ mapset_HINI = (char *) G_find_raster2(name_HINI, "");
+ if (mapset_HINI == NULL)
+ G_fatal_error(_("cell file [%s] not found"), name_HINI);
+
+ /* rast_open_old */
+ if ((infd_ELEV = Rast_open_old(name_ELEV, mapset_ELEV)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), name_ELEV);
+
+
+ if ((infd_DIST = Rast_open_old(name_DIST, mapset_DIST)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), name_DIST);
+
+ if ((infd_HINI = Rast_open_old(name_HINI, mapset_HINI)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), name_HINI);
+
+
+ /* controlling if we can open input raster */
+ Rast_get_cellhd(name_ELEV, mapset_ELEV, &cellhd);
+ G_debug(3, "number of rows %d", cellhd.rows);
+
+ Rast_get_cellhd(name_DIST, mapset_DIST, &cellhd);
+ G_debug(3, "number of rows %d", cellhd.rows);
+
+ Rast_get_cellhd(name_HINI, mapset_HINI, &cellhd);
+ G_debug(3, "number of rows %d", cellhd.rows);
+
+
+ /* allocate input buffer */
+ inrast_ELEV = Rast_allocate_f_buf();
+ inrast_DIST = Rast_allocate_f_buf();
+ inrast_HINI = Rast_allocate_f_buf();
+
+ G_debug(1,"Getting region extension...");
+
+ /* get windows rows & cols */
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ G_get_window(&window);
+ res_ew = window.ew_res;
+ res_ns = window.ns_res;
+ ca = res_ns * res_ew;
+
+
+ G_debug(2,"Allocating memory matrix...");
+
+ /* allocate memory matrix */
+ m_ELEV = G_alloc_fmatrix(nrows, ncols);
+ m_OUTLET = G_alloc_imatrix(nrows, ncols);
+ m_DIST = G_alloc_fmatrix(nrows, ncols);
+ m_HINI = G_alloc_fmatrix(nrows, ncols);
+ m_gx = G_alloc_fmatrix(nrows, ncols);
+ m_gy = G_alloc_fmatrix(nrows, ncols);
+ m_slope = G_alloc_fmatrix(nrows, ncols);
+ m_U = G_alloc_fmatrix(nrows, ncols);
+ m_Uloop = G_alloc_fmatrix(nrows, ncols);
+ m_Uloop_dt = G_alloc_fmatrix(nrows, ncols);
+ m_V = G_alloc_fmatrix(nrows, ncols);
+ m_Vloop = G_alloc_fmatrix(nrows, ncols);
+ m_Vloop_dt = G_alloc_fmatrix(nrows, ncols);
+ m_H = G_alloc_fmatrix(nrows, ncols);
+ m_Hloop = G_alloc_fmatrix(nrows, ncols);
+ m_Hloop_dt = G_alloc_fmatrix(nrows, ncols);
+ m_Hland = G_alloc_fmatrix(nrows, ncols);
+ m_K = G_alloc_fmatrix(nrows, ncols);
+ m_Kloop = G_alloc_fmatrix(nrows, ncols);
+ m_HUloop = G_alloc_fmatrix(nrows, ncols);
+ m_HVloop = G_alloc_fmatrix(nrows, ncols);
+
+ if (NS_THRES!=-1){
+ m_Hold= G_alloc_fmatrix(nrows, ncols);
+ }
+
+ if (result_HMAX){
+ m_Hmax = G_alloc_fmatrix(nrows, ncols);
+ }
+
+ if (result_VMAX){
+ m_Vmax = G_alloc_fmatrix(nrows, ncols);
+ }
+
+ /* read rows */
+ G_debug(2,"Reading input maps...");
+ for (row = 0; row < nrows; row++) {
+
+ /* read a line input maps into buffers*/
+ Rast_get_f_row(infd_ELEV, inrast_ELEV, row);
+ Rast_get_f_row(infd_DIST, inrast_DIST, row);
+ Rast_get_f_row(infd_HINI, inrast_HINI, row);
+
+ for (col = 0; col < ncols; col++) {
+ elev_fcell_val = ((FCELL *) inrast_ELEV)[col];
+ dist_fcell_val = ((FCELL *) inrast_DIST)[col];
+ hini_fcell_val = ((FCELL *) inrast_HINI)[col];
+
+ /* elevation map */
+ if (Rast_is_f_null_value(&elev_fcell_val) == 1) {
+ m_ELEV[row][col] = nullo;
+ } else {
+ m_ELEV[row][col] = ((FCELL *) inrast_ELEV)[col];
+ }
+
+ /* distance map */
+ if (Rast_is_f_null_value(&dist_fcell_val) == 1) {
+ m_DIST[row][col] = nullo;
+ } else {
+ m_DIST[row][col] = ((FCELL *) inrast_DIST)[col];
+ }
+
+ /* hini map */
+ if (Rast_is_f_null_value(&hini_fcell_val) == 1) {
+ m_HINI[row][col] = nullo;
+ } else {
+ m_HINI[row][col] = ((FCELL *) inrast_HINI)[col];
+ }
+
+ }
+ }
+
+ /* memory cleanup */
+ G_free(inrast_ELEV);
+ G_free(inrast_DIST);
+ G_free(inrast_HINI);
+
+ /* closing raster map */
+ Rast_close(infd_ELEV);
+ Rast_close(infd_DIST);
+ Rast_close(infd_HINI);
+
+
+ /* earth pressure coefficient */
+ k_act = (1 - sin((M_PI * IFRICT) / 180.0))
+ / (1 + sin((M_PI * IFRICT) / 180.0));
+
+ k_pas = (1 + sin((M_PI * IFRICT) / 180.0))
+ / (1 - sin((M_PI * IFRICT) / 180.0));
+
+ G_debug(1,"K_pas value = %f",k_pas);
+ G_debug(1,"K_act value = %f",k_act);
+
+
+ /* inizializzazione matrici in memoria: U,V,H,gx,gz,slope */
+
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+
+ m_H[row][col] = 0;
+ m_Hloop[row][col] = 0;
+ m_Hloop_dt[row][col] = 0;
+ m_Hland[row][col] = 0;
+ m_U[row][col] = 0;
+ m_Uloop[row][col] = 0;
+ m_Uloop_dt[row][col] = 0;
+ m_V[row][col] = 0;
+ m_Vloop[row][col] = 0;
+ m_Vloop_dt[row][col] = 0;
+ m_K[row][col] = nullo;
+ m_gx[row][col] = nullo;
+ m_gy[row][col] = nullo;
+ m_slope[row][col] = nullo;
+ m_HUloop[row][col] = 0;
+ m_HVloop[row][col] = 0;
+ }
+ }
+
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ if (row==0 || row==1 || col==0 || col==1 || row==nrows-1 || col==ncols-1 || row==nrows-2 || col==ncols-2)
+ m_OUTLET[row][col]=1;
+ else
+ m_OUTLET[row][col]=0;
+ }
+ }
+
+
+ if (NS_THRES!=-1){
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ m_Hold[row][col] = 0;
+ }
+ }
+ }
+
+
+ if (result_HMAX){
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ m_Hmax[row][col] = 0;
+ }
+ }
+ }
+
+ if (result_VMAX){
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ m_Vmax[row][col] = 0;
+ }
+ }
+ }
+
+
+
+ /*---------------------------------------------------------- *\
+ /* calcolo matrici gx, gy, slope */
+
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+
+ gx = gradx2(m_ELEV, row, col, res_ew, 0);
+ gy = grady2(m_ELEV, row, col, res_ns, 0);
+
+ /* slope calcolata con Horn 3x3 */
+ slope = sqrt(pow(gradx3(m_ELEV, row, col, res_ew, 0), 2) + pow(grady3(m_ELEV, row, col, res_ns, 0), 2));
+ m_gx[row][col] = gx;
+ m_gy[row][col] = gy;
+ m_slope[row][col] = slope;
+
+ }
+ }
+
+
+ /* Starting loops*/
+ int t=1;
+ int i;
+ int loop, n_loops = 1;
+ float incr_volexp = 0;
+ float Vol_gone=0;
+ float Vol_exp=0;
+ float Vol_actual=0;
+ int NS_count = 0;
+
+
+ /* while: for each timestep or until NS_count<3 */
+ while (NS_count<3 && t<=TIMESTEPS) {
+ i = t - 1;
+
+ float CFL_max = 0;
+ int stable = 0;
+
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+
+ /* Ricalcolo H in funzione di H_fluid e calcolo incremento di Volume aspettato */
+ if (m_DIST[row][col] <= i * dT * FLUID
+ && m_DIST[row][col] > (i - 1) * dT * FLUID) {
+ H_fluid = m_HINI[row][col];
+ m_H[row][col] += H_fluid;
+ } else {
+ H_fluid = 0.0;
+ }
+ incr_volexp += H_fluid/(cos(atan(m_slope[row][col])))*ca;
+
+
+ /* Calcolo Hland */
+ if (m_DIST[row][col] > i * dT * FLUID) {
+ m_Hland[row][col] = m_HINI[row][col];
+ } else {
+ m_Hland[row][col] = 0.0;
+ }
+ }
+ }
+
+ G_debug(1,"\nTIMESTEP %i",t);
+ G_debug(1,"Volume expected increment = %f",incr_volexp);
+
+ /* ciclo while */
+ while (stable==0) {
+ int exit = 0;
+
+
+ G_debug(2,"Updating maps and volumes");
+
+ /* Aggiornamento Hloop,Vloop e Uloop */
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+ m_Hloop[row][col] = m_H[row][col];
+ m_Uloop[row][col] = m_U[row][col];
+ m_Vloop[row][col] = m_V[row][col];
+ }
+ }
+
+ /* Aggiornamento Volumi */
+ Vol_gone_loop = Vol_gone;
+ Vol_exp_loop = Vol_exp;
+ Vol_actual = 0;
+
+
+ /* Per ogni iterazione */
+ G_debug(1,"\n---NLOOPS=%i---",n_loops);
+ for (loop = 1; loop <= n_loops && exit == 0; loop++) {
+ G_debug(1,"\n-LOOP=%i",loop);
+
+
+ G_debug(2,"Updating K and Kloop matrix");
+
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+ if ((gradx2(m_Uloop, row, col, res_ew, 0)
+ + grady2(m_Vloop, row, col, res_ns, 0)) >= 0)
+ m_K[row][col] = k_act;
+ else
+ m_K[row][col] = k_pas;
+ }
+ }
+
+
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+ m_Kloop[row][col] = lax(m_K, row, col, laxfactor);
+ }
+ }
+
+
+ /* ciclo row-col per calcolare G, P, I, T, Tb */
+ G_debug(2,"Calculating G, P, I, T, Tb");
+
+ for (row = 1; row < nrows - 1 && exit == 0; row++) {
+ for (col = 1; col < ncols - 1 && exit == 0; col++) {
+
+ /* Calcolo G_x e G_y */
+
+ if (m_Hloop[row][col] > verysmall) {
+ G_x = -grav * (sin(atan(m_gx[row][col])));
+ } else {
+ G_x = 0.0;
+ }
+
+
+ if (m_Hloop[row][col] > verysmall) {
+ G_y = -grav * (sin(atan(m_gy[row][col])));
+ } else {
+ G_y = 0.0;
+ }
+
+
+ /* Calcolo I_x e I_y */
+
+ I_x = -((m_Uloop[row][col] * (cos(atan(m_gx[row][col])))
+ * gradx2(m_Uloop, row, col, res_ew, 1))
+ + (m_Vloop[row][col] * (cos(atan(m_gy[row][col])))
+ * grady2(m_Uloop, row, col, res_ew, 1)));
+
+ I_y = -((m_Uloop[row][col] * (cos(atan(m_gx[row][col])))
+ * gradx2(m_Vloop, row, col, res_ew, 1))
+ + (m_Vloop[row][col] * (cos(atan(m_gy[row][col])))
+ * grady2(m_Vloop, row, col, res_ns, 1)));
+
+
+ /* Calcolo P_x e P_y */
+
+ if (m_Hloop[row][col] > verysmall) {
+ P_x = -m_Kloop[row][col] * cos(atan(m_gx[row][col]))
+ * gradPx2(m_Hloop, m_Hland, m_gx, row, col,
+ res_ew);
+ } else {
+ P_x = 0.0;
+ }
+
+ if (m_Hloop[row][col] > verysmall) {
+ P_y = -m_Kloop[row][col] * cos(atan(m_gy[row][col]))
+ * gradPy2(m_Hloop, m_Hland, m_gy, row, col,
+ res_ns);
+ } else {
+ P_y = 0.0;
+ }
+
+
+ /* Calcolo T_x e T_y */
+
+ vel = sqrt(
+ pow(m_Uloop[row][col], 2)
+ + pow(m_Vloop[row][col], 2));
+
+
+ if (RHEOL_TYPE == 1) {
+ T = t_frict(m_Hloop, row, col, BFRICT);
+ }
+
+ if (RHEOL_TYPE == 2) {
+ T = t_voellmy(vel, m_Hloop, row, col, BFRICT, CHEZY);
+ }
+
+ if (RHEOL_TYPE == 3) {
+ T = t_visco(vel, m_Hloop, row, col, BFRICT, RHO, VISCO,
+ YSTRESS);
+ }
+
+ if (m_Hloop[row][col] > verysmall && vel > verysmall) {
+ T_x = fabs(m_Uloop[row][col]) / vel * grav
+ * (cos(atan(m_gx[row][col]))) * T;
+ T_y = fabs(m_Vloop[row][col]) / vel * grav
+ * (cos(atan(m_gy[row][col]))) * T;
+ }
+
+ else {
+ T_x = grav * (cos(atan(m_gx[row][col]))) * T;
+ T_y = grav * (cos(atan(m_gy[row][col]))) * T;
+ }
+
+
+ /* Stima del flusso a t + ddt */
+ ddt = dTLeap * dT / n_loops;
+
+ /* Lax su Uloop e Vloop se is_df altrimenti 0 */
+
+ Uloop_a = filter_lax(m_Uloop, row, col, laxfactor, m_Hloop, 0.01, 0.0);
+ Vloop_a = filter_lax(m_Vloop, row, col, laxfactor, m_Hloop, 0.01, 0.0);
+
+ /* Calcolo U_loop_b e V_loop_b */
+
+ Uloop_b = veldt(Uloop_a, ddt, G_x, P_x, I_x, T_x);
+ Vloop_b = veldt(Vloop_a, ddt, G_y, P_y, I_y, T_y);
+
+
+ /* Calcolo T_x_b, T_y_b in funzione delle nuove Uloop_b, Vloop_b */
+
+ vel_b = sqrt(pow(Uloop_b, 2) + pow(Vloop_b, 2));
+
+
+ if (RHEOL_TYPE == 1) {
+ T_b = T;
+ }
+
+ if (RHEOL_TYPE == 2) {
+ T_b = t_voellmy(vel_b, m_Hloop, row, col, BFRICT,
+ CHEZY);
+ }
+
+ if (RHEOL_TYPE == 3) {
+ T_b = t_visco(vel_b, m_Hloop, row, col, BFRICT, RHO,
+ VISCO, YSTRESS);
+ }
+
+
+ if (m_Hloop[row][col] > verysmall && vel_b > verysmall) {
+ T_x_b = fabs(Uloop_b) / vel_b * grav
+ * (cos(atan(m_gx[row][col]))) * T_b;
+ T_y_b = fabs(Vloop_b) / vel_b * grav
+ * (cos(atan(m_gy[row][col]))) * T_b;
+ } else {
+ T_x_b = grav * (cos(atan(m_gx[row][col]))) * T_b;
+ T_y_b = grav * (cos(atan(m_gy[row][col]))) * T_b;
+ }
+
+
+ /* Stima del flusso a t + dt*/
+
+ dt = dT / n_loops;
+ Uloop_dt = veldt(Uloop_a, dt, G_x, P_x, I_x, T_x_b);
+ Vloop_dt = veldt(Vloop_a, dt, G_y, P_y, I_y, T_y_b);
+
+ /* Calcolo del valore Courant-Friendrichs-Levy */
+
+ CFL_u = dt * sqrt(2)
+ * fabs((cos(atan(m_gx[row][col]))) * Uloop_dt)
+ / res_ew;
+
+ CFL_v = dt * sqrt(2)
+ * fabs((cos(atan(m_gy[row][col]))) * Vloop_dt)
+ / res_ns;
+
+ CFL = max(CFL_u,CFL_v);
+
+
+ if (CFL > CFL_max)
+ CFL_max = CFL;
+
+ /* Verifica della condizione di stabilità */
+
+ if (CFL_max > CFLlimsup && n_loops < MaxNLoops) {
+ exit = 1;
+ n_loops += 1;
+ CFL_max=0;
+ } else {
+ m_Uloop_dt[row][col] = Uloop_dt;
+ m_Vloop_dt[row][col] = Vloop_dt;
+
+
+ if (m_Hloop[row][col] > verysmall) {
+ m_HUloop[row][col] = m_Hloop[row][col] * Uloop_dt;
+ m_HVloop[row][col] = m_Hloop[row][col] * Vloop_dt;
+ } else {
+ m_HUloop[row][col] = 0.0;
+ m_HVloop[row][col] = 0.0;
+ }
+ }
+ } /*chiusura FOR cols*/
+ } /*chiusura FOR rows*/
+
+ G_debug(2,"CFL_max=%f",CFL_max);
+
+ if (exit==0) {
+ G_debug(2,"Calculating Hloop_dt without mbe for loop");
+
+ for (row = 1; row < nrows - 1 && exit == 0; row++) {
+ for (col = 1; col < ncols - 1 && exit == 0; col++) {
+
+ /* calcolo dH/dT */
+ dH_dT = -cos(atan(m_gx[row][col]))* (shift0(m_HUloop,row,col,nrows-1,ncols-1,1,1,0,-1)-shift0(m_HUloop,row,col,nrows-1,ncols-1,1,1,0,1))/(2*res_ew/cos(atan(m_gx[row][col])))
+ -cos(atan(m_gy[row][col]))* (shift0(m_HVloop,row,col,nrows-1,ncols-1,1,1,1,0)-shift0(m_HVloop,row,col,nrows-1,ncols-1,1,1,-1,0))/(2*res_ns/cos(atan(m_gy[row][col])));
+
+ /* Lax su Hloop e calcolo Hloop_dt (senza mbe) */
+
+ if (dH_dT==0) {
+ Hloop_a=m_Hloop[row][col];
+ } else {
+ Hloop_a = filter_lax(m_Hloop, row, col, laxfactor, m_Hloop, 0.01, 0.0);
+ }
+
+ Hloop_dt = Hloop_a - dT/n_loops * dH_dT;
+
+
+ /* Vol_gone_loop */
+ Vol_actual += Hloop_dt/(cos(atan(m_slope[row][col])))*ca;
+ if (m_OUTLET[row][col]==1){
+ /*Vol_gone_no_mbe = Hloop_dt/(cos(atan(m_slope[row][col])))*ca;*/
+ Vol_gone_loop += Hloop_dt/(cos(atan(m_slope[row][col])))*ca;
+ Hloop_dt = 0.0;
+ }
+
+ /* matrice H_loop_dtemp */
+ m_Hloop_dt[row][col] = Hloop_dt;
+
+ }/*chiusura FOR col*/
+ }/*chiusura FOR row*/
+
+
+ /* m.b.e */
+ G_debug(2,"Calculating mass balance error");
+
+ Vol_exp_loop += incr_volexp/n_loops;
+
+ if (Vol_exp_loop > 0.01) {
+ mbe = Vol_actual/Vol_exp_loop;
+ } else {
+ mbe = 1;
+ }
+
+ G_debug(2,"mbe=%f",mbe);
+
+ /*Vol_gone_loop += Vol_gone_no_mbe/mbe;*/
+ Vol_actual = 0;
+
+ /* Hloop_dt con mbe */
+ G_debug(2,"Calculating Hloop_dt with mbe for loop");
+
+ for (row = 1; row < nrows - 1 && exit == 0; row++) {
+ for (col = 1; col < ncols - 1 && exit == 0; col++) {
+ if (mbe > 0.01){
+ m_Hloop_dt[row][col]=m_Hloop_dt[row][col]/mbe;
+ }
+ }
+ }
+
+ /* loop<n_loops */
+ if (loop<n_loops) {
+ G_debug(2,"loop<n_loops");
+
+ for (row = 1; row < nrows - 1 && exit == 0; row++) {
+ for (col = 1; col < ncols - 1 && exit == 0; col++) {
+ m_Hloop[row][col]=m_Hloop_dt[row][col];
+ m_Uloop[row][col]=m_Uloop_dt[row][col];
+ m_Vloop[row][col]=m_Vloop_dt[row][col];
+ }
+ }
+ }
+
+ }/* chiusura IF exit*/
+
+ } /*chiusura FOR loops*/
+
+ if (exit==0) {
+ G_debug(2,"FINISH N_LOOPS");
+
+
+ /*------Nash-Sutcliffe-------*/
+ if (NS_THRES!=-1){
+ G_debug(1,"Calculating Nash-Sutcliffe value");
+ if (t==1){
+ NSold = nash_sutcliffe(m_Hold,m_Hloop_dt,nrows, ncols);
+ } else {
+ NS = nash_sutcliffe(m_Hold,m_Hloop_dt, nrows, ncols);
+ G_debug(1,"NS old value=%f",NSold);
+ G_debug(1,"NS value=%f",NS);
+
+ if(fabs(NS-NSold) <= NS_THRES){
+ NS_count += 1;
+ } else {
+ NS_count = 0;
+ }
+ NSold = NS;
+ }
+ G_debug(1,"NS count=%i",NS_count);
+ }
+
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+ m_H[row][col]=m_Hloop_dt[row][col];
+ m_U[row][col]=m_Uloop_dt[row][col];
+ m_V[row][col]=m_Vloop_dt[row][col];
+ if (NS_THRES!=-1){
+ m_Hold[row][col]=m_Hloop_dt[row][col];
+ }
+ }
+ }
+
+ if (result_HMAX){
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+ if(m_H[row][col]+m_Hland[row][col] > m_Hmax[row][col])
+ m_Hmax[row][col] = m_H[row][col]+m_Hland[row][col];
+ }
+ }
+ }
+
+ if (result_VMAX){
+ for (row = 1; row < nrows - 1; row++) {
+ for (col = 1; col < ncols - 1; col++) {
+ if (m_H[row][col]>verysmall){
+ if ((sqrt(pow(m_U[row][col],2)+pow(m_V[row][col],2)))>m_Vmax[row][col]){
+ m_Vmax[row][col]=(sqrt(pow(m_U[row][col],2)+pow(m_V[row][col],2)));
+ }
+ }
+ }
+ }
+ }
+
+ stable=1;
+ G_debug(1,"Vol_gone_loop=%f\n",Vol_gone_loop);
+
+
+ Vol_exp=Vol_exp_loop;
+ Vol_gone=Vol_gone_loop;
+
+ incr_volexp=0.0;
+ if (CFL_max<=CFLliminf) {
+ n_loops=max(n_loops-1,MinNLoops);
+ }
+
+ /*------- DEFAULT OUTPUT SERIES--------*/
+
+ if((DELTAT!=-1 && t%DELTAT==0) || (t==TIMESTEPS) || (NS_count==3)){
+
+ if(result_H){
+ sprintf(name1,"%s_t%d",result_H,t);
+ out_sum_print(m_Hland,m_H,m_U,m_V,name1,nrows,ncols,1,0.01);
+ Rast_short_history(name1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(name1, &history);
+ }
+ if(result_V){
+ sprintf(name1,"%s_t%d",result_V,t);
+ out_sum_print(m_Hland,m_H,m_U,m_V,name1,nrows,ncols,2,0.01);
+ Rast_short_history(name1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(name1, &history);
+ }
+ }
+
+
+ if((result_HMAX) && ((t==TIMESTEPS) || (NS_count==3))){
+ sprintf(name1,"%s",result_HMAX);
+ out_print(m_Hmax,name1,nrows,ncols);
+ Rast_short_history(name1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(name1, &history);
+ }
+
+ if((result_VMAX) && ((t==TIMESTEPS) || (NS_count==3))){
+ sprintf(name1,"%s",result_VMAX);
+ out_print(m_Vmax,name1,nrows,ncols);
+ Rast_short_history(name1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(name1, &history);
+ }
+
+
+ } /*chiusura IF exit*/
+ } /*chiusura WHILE*/
+ t++;
+ } /* chiusura WHILE NS_count*/
+
+ if (Vol_gone>0){
+ G_warning("The flow is out of the current computational domain. You may need to change the region extension.");
+ G_warning("Volume gone: %f",Vol_gone);
+ }
+
+
+ /* deallocate memory matrix */
+ G_free_fmatrix(m_ELEV);
+ G_free_imatrix(m_OUTLET);
+ G_free_fmatrix(m_DIST);
+ G_free_fmatrix(m_HINI);
+ G_free_fmatrix(m_gx);
+ G_free_fmatrix(m_gy);
+ G_free_fmatrix(m_slope);
+ G_free_fmatrix(m_H);
+ G_free_fmatrix(m_Hloop);
+ G_free_fmatrix(m_Hloop_dt);
+ G_free_fmatrix(m_Hland);
+ G_free_fmatrix(m_U);
+ G_free_fmatrix(m_Uloop);
+ G_free_fmatrix(m_Uloop_dt);
+ G_free_fmatrix(m_V);
+ G_free_fmatrix(m_Vloop);
+ G_free_fmatrix(m_Vloop_dt);
+ G_free_fmatrix(m_K);
+ G_free_fmatrix(m_Kloop);
+ G_free_fmatrix(m_HUloop);
+ G_free_fmatrix(m_HVloop);
+
+ if (NS_THRES!=-1)
+ G_free_fmatrix(m_Hold);
+
+ if (result_VMAX)
+ G_free_fmatrix(m_Vmax);
+
+ if (result_HMAX)
+ G_free_fmatrix(m_Hmax);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/raster/r.massvov/r.massmov.html
===================================================================
--- grass-addons/grass7/raster/r.massvov/r.massmov.html (rev 0)
+++ grass-addons/grass7/raster/r.massvov/r.massmov.html 2012-05-25 11:07:09 UTC (rev 51734)
@@ -0,0 +1,78 @@
+<h2>DESCRIPTION</h2>
+
+<p>MassMov2D is a numerical model that allows users to simulate the expansion (runout)
+and deposition of mass movements over a complex topography by approximating the
+heterogeneous sliding mass to a homogeneous one-phase fluid (following the approach
+proposed by Savage and Hutter (1989) and Iverson and Denlinger (2001).
+The model describes the mass movements as a two-dimensional flux taking advantage
+of the shallow water equations. This formula is derived from the general Navier-
+Stokes equations under the hypothesis that the vertical components of velocity and
+pressure are negligible with respect to the horizontal components, and that the vertical
+pressure profile can be considered as almost hydrostatic (Kinnmark 1985).</p>
+
+<p>The required inputs can be classified in three categories based on the information type:
+<ul>
+<li> raster maps of the topography, in particular the sliding surface topography (digital terrain
+model without the sliding body), the initial sliding mass thickness and the 'distance map',
+representing the cells distance from the collapsing body lower limit;</li>
+<li> numerical parameters for the characterization of the mass material, density [kg/m3], apparent
+yield stress [Pa], Chezy roughness coefficient [m/s2], dynamic viscosity [Pa*s], basal friction
+angle [deg], internal friction angle of the sliding mass during the expansion [deg] and the fluid rate [m/s].This last parameter provides information on the transaction velocity of the
+sliding mass when passes from a solid state to a fluid state; together with the 'distance map'
+it allows to define the amount of mass mobilized as a function of time. It is worth noting that
+depending on the selected rheological law different sets of parameters are mandatory;</li>
+<li>control parameters that are used to stop the simulation, like maximum time step number and/or Nash-Sutcliff tresholds value.</li>
+</ul>
+
+<p>The model outputs a series of flux velocity map and deposit depth raster map at different time step according to the setted <i>deltatime</i> parameter; Additionally the module outputs two raster maps representing the maximum thickness and velocity registered during the simulation.</p>
+
+<h2>NOTES</h2>
+
+<p>The generation of the model input maps, in case the simulation refer to en existing collapse and pre and post event DTM is available, can be performed taking advantage of the GRASS modules; in
+particular:</p>
+<ul>
+<li>the sliding surface can be calculated by subtracting the collapsing body from the pre-event
+DTM (r.mapcalc) </li>
+<li>the collapsing body thickness can evaluated by considering the negative differences between
+the post and pre-event DTM multiplied for the cosine of the slope (r.mapcalc and
+r.slope.aspect)</li>
+<li>the distance map from the landslide toe can be obtained by applying the r.grow.distance module to the rasterized limits of the landslide</li>
+</ul>
+
+<h2>DIAGNOSTICS</h2>
+
+<p>The module has been tested in several cases (see references), but up to now most of the simulations was done using a Voellmy rheology thus other rheology laws should be better investigated.</p>
+
+
+<h2>REFERENCES</h2>
+
+<p>Savage S B and Hutter K 1989 The motion of a finite mass of granular material down a rough
+incline. Journal of Fluid Mechanics 199:177-215</p>
+
+<p>Iverson R M and Denlinger R P 2001 Flow of variably fluidized granular masses across
+threedimensional terrain: 1, Coulomb mixture theory. Journal of Geophysical Research 106:537-52</p>
+
+<p>Kinnmark I P E 1985 The shallow water equations: Formulation, analysis and application. In
+Brebia C A and Orszag S A (eds) Lecture Notes in Engineering 15. Berlin, Springer-Verlag:1-187</p>
+
+<p>Begueria S, Van Asch T W J, Malet J P and Grondahl S 2009 A GIS based numerical model for simulating the kinematics
+of mud and debris flows over complex terrain. Nat Hazards Earth Syst Sci, 9, 1897-1909.</p>
+
+<p>Molinari M, Cannata M, Begueria S and Ambrosi C 2012 GIS-based Calibration of MassMov2D.
+Transactions in GIS, 2012, 16(2):215-231</p>
+
+<h2>SEE ALSO</h2>
+
+<a href="r.grow.distance.html">r.grow.distance</a><br>
+<a href="r.slope.aspect.html">r.slope.aspect</a><br>
+<a href="r.mapcalc.html">r.mapcalc</a><br>
+
+<h2>AUTHORS</h2>
+
+<p><i>Original version of program:</i>
+<br>Santiago Begueria
+
+<p><i>The current version of the program (ported to GRASS7.0)</i>:
+<br>Monia Molinari, massimiliano Cannata, Santiago Begueria.<br>
+
+<p><i>Last changed: $Date: 2011-11-08 13:24:20 -0800 (Tue, 08 Nov 2011) $</i>
More information about the grass-commit
mailing list