GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
tin.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/tin.c
3 
4  \brief Vector library - TIN
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  (C) 2001-2009 by the GRASS Development Team
9 
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12 
13  \author Radim Blazek
14  */
15 
16 #include <grass/vector.h>
17 
18 /*!
19  \brief Calculates z coordinate for point from TIN
20 
21  \param Map pointer to vector map
22  \param tx,ty point coordinates
23  \param[out] tz z-coordinate of point
24  \param[out] angle angle
25  \param[out] slope slope
26 
27  \return 1 on success,
28  \return 0 point is not in area,
29  \return -1 area has not 4 points or has island
30  */
31 int
33  double tx, double ty, double *tz, double *angle, double *slope)
34 {
35  int i, area, n_points;
36  struct Plus_head *Plus;
37  struct P_area *Area;
38  static struct line_pnts *Points;
39  static int first_time = 1;
40  double *x, *y, *z;
41  double vx1, vx2, vy1, vy2, vz1, vz2;
42  double a, b, c, d;
43 
44  /* TODO angle, slope */
45 
46  Plus = &(Map->plus);
47  if (first_time == 1) {
48  Points = Vect_new_line_struct();
49  first_time = 0;
50  }
51 
52  area = Vect_find_area(Map, tx, ty);
53  G_debug(3, "TIN: area = %d", area);
54  if (area == 0)
55  return 0;
56 
57  Area = Plus->Area[area];
58  if (Area->n_isles > 0)
59  return -1;
60 
61  Vect_get_area_points(Map, area, Points);
62  n_points = Points->n_points;
63  if (n_points != 4)
64  return -1;
65 
66  x = Points->x;
67  y = Points->y;
68  z = Points->z;
69  for (i = 0; i < 3; i++) {
70  G_debug(3, "TIN: %d %f %f %f", i, x[i], y[i], z[i]);
71  }
72 
73  vx1 = x[1] - x[0];
74  vy1 = y[1] - y[0];
75  vz1 = z[1] - z[0];
76  vx2 = x[2] - x[0];
77  vy2 = y[2] - y[0];
78  vz2 = z[2] - z[0];
79 
80  a = vy1 * vz2 - vy2 * vz1;
81  b = vz1 * vx2 - vz2 * vx1;
82  c = vx1 * vy2 - vx2 * vy1;
83  d = -a * x[0] - b * y[0] - c * z[0];
84 
85  /* OK ? */
86  *tz = -(d + a * tx + b * ty) / c;
87  G_debug(3, "TIN: z = %f", *tz);
88 
89  return 1;
90 }
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:891
int n_points
Number of points.
Definition: dig_structs.h:1692
#define x
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
Basic topology-related info.
Definition: dig_structs.h:784
double b
Definition: r_raster.c:39
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1286
plus_t n_isles
Number of islands inside.
Definition: dig_structs.h:1632
int Vect_find_area(struct Map_info *, double, double)
Find the nearest area.
Vector map info.
Definition: dig_structs.h:1259
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
Area (topology) info.
Definition: dig_structs.h:1605
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
int Vect_get_area_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_tin_get_z(struct Map_info *Map, double tx, double ty, double *tz, double *angle, double *slope)
Calculates z coordinate for point from TIN.
Definition: tin.c:32