GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * Written by H. Mitasova, L. Mitas, I. Kosinovsky, D. Gerdes Fall 1994 00004 * University of Illinois 00005 * US Army Construction Engineering Research Lab 00006 * Copyright 1994, H. Mitasova (University of Illinois), 00007 * L. Mitas (University of Illinois), 00008 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL) 00009 * 00010 * modified by McCauley in August 1995 00011 * modified by Mitasova in August 1995 00012 * 00013 */ 00014 00015 #include <stdio.h> 00016 #include <math.h> 00017 #include <unistd.h> 00018 #include <grass/gis.h> 00019 #include <grass/bitmap.h> 00020 #include <grass/interpf.h> 00021 00022 int IL_secpar_loop_2d(struct interp_params *params, int ngstc, /* starting column */ 00023 int nszc, /* ending column */ 00024 int k, /* current row */ 00025 struct BM *bitmask, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, /* min,max interp. 00026 * values */ 00027 int cond1, int cond2 /* determine if particular values need to 00028 * be computed */ 00029 ) 00030 00031 /* 00032 * Computes slope, aspect and curvatures (depending on cond1, cond2) for 00033 * derivative arrays adx,...,adxy between columns ngstc and nszc. 00034 */ 00035 { 00036 double dnorm1, ro, /* rad to deg conv */ 00037 dx2 = 0, dy2 = 0, grad2 = 0, /* gradient squared */ 00038 slp = 0, grad, /* gradient */ 00039 oor = 0, /* aspect (orientation) */ 00040 curn = 0, /* profile curvature */ 00041 curh = 0, /* tangential curvature */ 00042 curm = 0, /* mean curvature */ 00043 temp, /* temp variable */ 00044 dxy2; /* temp variable square of part diriv. */ 00045 00046 double gradmin; 00047 int i, got, bmask = 1; 00048 static int first_time_g = 1; 00049 00050 ro = 57.295779; 00051 gradmin = 0.001; 00052 00053 00054 for (i = ngstc; i <= nszc; i++) { 00055 if (bitmask != NULL) { 00056 bmask = BM_get(bitmask, i, k); 00057 } 00058 got = 0; 00059 if (bmask == 1) { 00060 while ((got == 0) && (cond1)) { 00061 dx2 = (double)(params->adx[i] * params->adx[i]); 00062 dy2 = (double)(params->ady[i] * params->ady[i]); 00063 grad2 = dx2 + dy2; 00064 grad = sqrt(grad2); 00065 /* slope in % slp = 100. * grad; */ 00066 /* slope in degrees */ 00067 slp = ro * atan(grad); 00068 if (grad <= gradmin) { 00069 oor = 0.; 00070 got = 3; 00071 if (cond2) { 00072 curn = 0.; 00073 curh = 0.; 00074 got = 3; 00075 break; 00076 } 00077 } 00078 if (got == 3) 00079 break; 00080 00081 /***********aspect from r.slope.aspect, with adx, ady computed 00082 from interpol. function RST **************************/ 00083 00084 if (params->adx[i] == 0.) { 00085 if (params->ady[i] > 0.) 00086 oor = 90; 00087 else 00088 oor = 270; 00089 } 00090 else { 00091 oor = ro * atan2(params->ady[i], params->adx[i]); 00092 if (oor <= 0.) 00093 oor = 360. + oor; 00094 } 00095 00096 got = 1; 00097 } /* while */ 00098 if ((got != 3) && (cond2)) { 00099 00100 dnorm1 = sqrt(grad2 + 1.); 00101 dxy2 = 00102 2. * (double)(params->adxy[i] * params->adx[i] * 00103 params->ady[i]); 00104 00105 00106 curn = 00107 (double)(params->adxx[i] * dx2 + dxy2 + 00108 params->adyy[i] * dy2) / (grad2 * dnorm1 * 00109 dnorm1 * dnorm1); 00110 00111 curh = 00112 (double)(params->adxx[i] * dy2 - dxy2 + 00113 params->adyy[i] * dx2) / (grad2 * dnorm1); 00114 00115 temp = grad2 + 1.; 00116 curm = 00117 .5 * ((1. + dy2) * params->adxx[i] - dxy2 + 00118 (1. + dx2) * params->adyy[i]) / (temp * dnorm1); 00119 } 00120 if (first_time_g) { 00121 first_time_g = 0; 00122 *gmin = *gmax = slp; 00123 *c1min = *c1max = curn; 00124 *c2min = *c2max = curh; 00125 } 00126 *gmin = amin1(*gmin, slp); 00127 *gmax = amax1(*gmax, slp); 00128 *c1min = amin1(*c1min, curn); 00129 *c1max = amax1(*c1max, curn); 00130 *c2min = amin1(*c2min, curh); 00131 *c2max = amax1(*c2max, curh); 00132 if (cond1) { 00133 params->adx[i] = (FCELL) slp; 00134 params->ady[i] = (FCELL) oor; 00135 if (cond2) { 00136 params->adxx[i] = (FCELL) curn; 00137 params->adyy[i] = (FCELL) curh; 00138 params->adxy[i] = (FCELL) curm; 00139 } 00140 } 00141 } /* bmask == 1 */ 00142 } 00143 return 1; 00144 }