GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-6596d49e63
secpar2d.c
Go to the documentation of this file.
1 /*!
2  * \file secpar2d.c
3  *
4  * \author H. Mitasova, L. Mitas, I. Kosinovsky, D. Gerdes Fall 1994 (original
5  * authors) \author modified by McCauley in August 1995 \author modified by
6  * Mitasova in August 1995 \author H. Mitasova (University of Illinois) \author
7  * L. Mitas (University of Illinois) \author I. Kosinovsky, (USA-CERL) \author
8  * D.Gerdes (USA-CERL)
9  *
10  * \copyright
11  * (C) 1994-1995 by Helena Mitasova and the GRASS Development Team
12  *
13  * \copyright
14  * This program is free software under the
15  * GNU General Public License (>=v2).
16  * Read the file COPYING that comes with GRASS
17  * for details.
18  *
19  */
20 
21 #include <stdio.h>
22 #include <math.h>
23 #include <unistd.h>
24 #include <grass/gis.h>
25 #include <grass/bitmap.h>
26 #include <grass/interpf.h>
27 
28 /*!
29  * Compute slope aspect and curvatures
30  *
31  * Computes slope, aspect and curvatures (depending on cond1, cond2) for
32  * derivative arrays adx,...,adxy between columns ngstc and nszc.
33  */
35  struct interp_params *params, int ngstc, /*!< starting column */
36  int nszc, /*!< ending column */
37  int k, /*!< current row */
38  struct BM *bitmask, double *gmin, double *gmax, double *c1min,
39  double *c1max, double *c2min, double *c2max, /*!< min,max interp. values */
40  int cond1,
41  int cond2 /*!< determine if particular values need to be computed */
42 )
43 {
44  double dnorm1, ro, /* rad to deg conv */
45  dx2 = 0, dy2 = 0, grad2 = 0, /* gradient squared */
46  slp = 0, grad, /* gradient */
47  oor = 0, /* aspect (orientation) */
48  curn = 0, /* profile curvature */
49  curh = 0, /* tangential curvature */
50  curm = 0, /* mean curvature */
51  temp, /* temp variable */
52  dxy2; /* temp variable square of part diriv. */
53 
54  double gradmin;
55  int i, got, bmask = 1;
56  static int first_time_g = 1;
57 
58  ro = M_R2D;
59  gradmin = 0.001;
60 
61  for (i = ngstc; i <= nszc; i++) {
62  if (bitmask != NULL) {
63  bmask = BM_get(bitmask, i, k);
64  }
65  got = 0;
66  if (bmask == 1) {
67  while ((got == 0) && (cond1)) {
68  dx2 = (double)(params->adx[i] * params->adx[i]);
69  dy2 = (double)(params->ady[i] * params->ady[i]);
70  grad2 = dx2 + dy2;
71  grad = sqrt(grad2);
72  /* slope in % slp = 100. * grad; */
73  /* slope in degrees */
74  slp = ro * atan(grad);
75  if (grad <= gradmin) {
76  oor = 0.;
77  got = 3;
78  if (cond2) {
79  curn = 0.;
80  curh = 0.;
81  got = 3;
82  break;
83  }
84  }
85  if (got == 3)
86  break;
87 
88  /***********aspect from r.slope.aspect, with adx, ady computed
89  from interpol. function RST
90  **************************/
91 
92  if (params->adx[i] == 0.) {
93  if (params->ady[i] > 0.)
94  oor = 90;
95  else
96  oor = 270;
97  }
98  else {
99  oor = ro * atan2(params->ady[i], params->adx[i]);
100  if (oor <= 0.)
101  oor = 360. + oor;
102  }
103 
104  got = 1;
105  } /* while */
106  if ((got != 3) && (cond2)) {
107 
108  dnorm1 = sqrt(grad2 + 1.);
109  dxy2 = 2. * (double)(params->adxy[i] * params->adx[i] *
110  params->ady[i]);
111 
112  curn = (double)(params->adxx[i] * dx2 + dxy2 +
113  params->adyy[i] * dy2) /
114  (grad2 * dnorm1 * dnorm1 * dnorm1);
115 
116  curh = (double)(params->adxx[i] * dy2 - dxy2 +
117  params->adyy[i] * dx2) /
118  (grad2 * dnorm1);
119 
120  temp = grad2 + 1.;
121  curm = .5 *
122  ((1. + dy2) * params->adxx[i] - dxy2 +
123  (1. + dx2) * params->adyy[i]) /
124  (temp * dnorm1);
125  }
126  if (first_time_g) {
127  first_time_g = 0;
128  *gmin = *gmax = slp;
129  *c1min = *c1max = curn;
130  *c2min = *c2max = curh;
131  }
132  *gmin = amin1(*gmin, slp);
133  *gmax = amax1(*gmax, slp);
134  *c1min = amin1(*c1min, curn);
135  *c1max = amax1(*c1max, curn);
136  *c2min = amin1(*c2min, curh);
137  *c2max = amax1(*c2max, curh);
138  if (cond1) {
139  params->adx[i] = (FCELL)slp;
140  params->ady[i] = (FCELL)oor;
141  if (cond2) {
142  params->adxx[i] = (FCELL)curn;
143  params->adyy[i] = (FCELL)curh;
144  params->adxy[i] = (FCELL)curm;
145  }
146  }
147  } /* bmask == 1 */
148  }
149  return 1;
150 }
#define NULL
Definition: ccmath.h:32
int BM_get(struct BM *, int, int)
Gets 'val' from the bitmap.
Definition: bitmap.c:217
float FCELL
Definition: gis.h:629
#define M_R2D
Definition: gis.h:167
double amin1(double, double)
Definition: minmax.c:65
double amax1(double, double)
Definition: minmax.c:52
int IL_secpar_loop_2d(struct interp_params *params, int ngstc, int nszc, int k, struct BM *bitmask, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, int cond1, int cond2)
Definition: secpar2d.c:34
Definition: bitmap.h:17
DCELL * adxy
Definition: interpf.h:89
DCELL * adyy
Definition: interpf.h:89
DCELL * adx
Definition: interpf.h:89
DCELL * ady
Definition: interpf.h:89
DCELL * adxx
Definition: interpf.h:89