GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
n_geom.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: part of the gpde library
9 * allocation, destroying and initializing the geometric struct
10 *
11 * COPYRIGHT: (C) 2000 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18 
19 
20 #include <grass/N_pde.h>
21 
22 /* *************************************************************** *
23  * *********** Konstruktor *************************************** *
24  * *************************************************************** */
25 /*!
26  * \brief Allocate the pde geometry data structure and return a pointer to the new allocated structure
27  *
28  * \return N_geom_data *
29  * */
31 {
32  N_geom_data *geom = (N_geom_data *) G_calloc(1, sizeof(N_geom_data));
33 
34  geom->area = NULL;
35  geom->planimetric = 1;
36  geom->dim = 0;
37 
38  return geom;
39 }
40 
41 /* *************************************************************** *
42  * *********** Destruktor **************************************** *
43  * *************************************************************** */
44 /*!
45  * \brief Release memory of a pde geometry data structure
46  *
47  * \param geom N_geom_data *
48  * \return void
49  * */
51 {
52  if (geom->area != NULL)
53  G_free(geom->area);
54 
55  G_free(geom);
56  return;
57 }
58 
59 /* *************************************************************** *
60  * *************************************************************** *
61  * *************************************************************** */
62 /*!
63  * \brief Initiate a pde geometry data structure with a 3d region
64  *
65  * If the projection is not planimetric, a double array will be created based on the
66  * number of rows of the provided region
67  *
68  * \param region3d RASTER3D_Region *
69  * \param geodata N_geom_data * - if a NULL pointer is given, a new structure will be allocatet and returned
70  *
71  * \return N_geom_data *
72  * */
74 {
75  N_geom_data *geom = geodata;
76  struct Cell_head region2d;
77 
78 #pragma omp critical
79  {
80 
81  G_debug(2,
82  "N_init_geom_data_3d: initializing the geometry structure");
83 
84  if (geom == NULL)
85  geom = N_alloc_geom_data();
86 
87  geom->dz = region3d->tb_res * G_database_units_to_meters_factor(); /*this function is not thread safe */
88  geom->depths = region3d->depths;
89  geom->dim = 3;
90 
91  /*convert the 3d into a 2d region and begin the area calculation */
92  G_get_set_window(&region2d); /*this function is not thread safe */
93  Rast3d_region_to_cell_head(region3d, &region2d);
94  }
95 
96  return N_init_geom_data_2d(&region2d, geom);
97 }
98 
99 
100 /* *************************************************************** *
101  * *************************************************************** *
102  * *************************************************************** */
103 /*!
104  * \brief Initiate a pde geometry data structure with a 2d region
105  *
106  * If the projection is not planimetric, a double array will be created based on the
107  * number of rows of the provided region storing all computed areas for each row
108  *
109  * \param region sruct Cell_head *
110  * \param geodata N_geom_data * - if a NULL pointer is given, a new structure will be allocatet and returned
111  *
112  * \return N_geom_data *
113  * */
115  N_geom_data * geodata)
116 {
117  N_geom_data *geom = geodata;
118  struct Cell_head backup;
119  double meters;
120  short ll = 0;
121  int i;
122 
123 
124  /*create an openmp lock to assure that only one thread at a time will access this function */
125 #pragma omp critical
126  {
127  G_debug(2,
128  "N_init_geom_data_2d: initializing the geometry structure");
129 
130  /*make a backup from this region */
131  G_get_set_window(&backup); /*this function is not thread safe */
132  /*set the current region */
133  Rast_set_window(region); /*this function is not thread safe */
134 
135  if (geom == NULL)
136  geom = N_alloc_geom_data();
137 
138  meters = G_database_units_to_meters_factor(); /*this function is not thread safe */
139 
140  /*set the dim to 2d if it was not initiated with 3, that's a bit ugly :( */
141  if (geom->dim != 3)
142  geom->dim = 2;
143 
144  geom->planimetric = 1;
145  geom->rows = region->rows;
146  geom->cols = region->cols;
147  geom->dx = region->ew_res * meters;
148  geom->dy = region->ns_res * meters;
149  geom->Az = geom->dy * geom->dx; /*square meters in planimetric proj */
150  /*depths and dz are initialized with a 3d region */
151 
152  /*Begin the area calculation */
153  ll = G_begin_cell_area_calculations(); /*this function is not thread safe */
154 
155  /*if the projection is not planimetric, calc the area for each row */
156  if (ll == 2) {
157  G_debug(2,
158  "N_init_geom_data_2d: calculating the areas for non parametric projection");
159  geom->planimetric = 0;
160 
161  if (geom->area != NULL)
162  G_free(geom->area);
163  else
164  geom->area = G_calloc(geom->rows, sizeof(double));
165 
166  /*fill the area vector */
167  for (i = 0; i < geom->rows; i++) {
168  geom->area[i] = G_area_of_cell_at_row(i); /*square meters */
169  }
170  }
171 
172  /*restore the old region */
173  Rast_set_window(&backup); /*this function is not thread safe */
174  }
175 
176  return geom;
177 }
178 
179 /* *************************************************************** *
180  * *************************************************************** *
181  * *************************************************************** */
182 /*!
183  * \brief Get the areay size in square meter of one cell (x*y) at row
184  *
185  * This function works for two and three dimensions
186  *
187  * \param geom N_geom_data *
188  * \param row int
189  * \return area double
190  *
191  * */
193 {
194  if (geom->planimetric) {
195  G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->Az);
196  return geom->Az;
197  }
198  else {
199  G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->area[row]);
200  return geom->area[row];
201  }
202 
203  return 0.0;
204 }
double tb_res
Definition: raster3d.h:57
2D/3D raster map header (used also for region)
Definition: gis.h:412
int dim
Definition: N_pde.h:107
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
#define NULL
Definition: ccmath.h:32
int depths
Definition: N_pde.h:115
#define G_calloc(m, n)
Definition: defs/gis.h:113
double dy
Definition: N_pde.h:110
void Rast_set_window(struct Cell_head *)
Establishes &#39;window&#39; as the current working window.
N_geom_data * N_init_geom_data_2d(struct Cell_head *region, N_geom_data *geodata)
Initiate a pde geometry data structure with a 2d region.
Definition: n_geom.c:114
Geometric information about the structured grid.
Definition: N_pde.h:103
double dz
Definition: N_pde.h:111
int rows
Definition: N_pde.h:116
double Az
Definition: N_pde.h:113
void G_get_set_window(struct Cell_head *)
Get the current working window (region)
int cols
Definition: N_pde.h:117
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
Get the areay size in square meter of one cell (x*y) at row.
Definition: n_geom.c:192
int G_begin_cell_area_calculations(void)
Begin cell area calculations.
Definition: gis/area.c:46
int planimetric
Definition: N_pde.h:105
N_geom_data * N_alloc_geom_data(void)
Allocate the pde geometry data structure and return a pointer to the new allocated structure...
Definition: n_geom.c:30
int cols
Number of columns for 2D data.
Definition: gis.h:431
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:452
void N_free_geom_data(N_geom_data *geom)
Release memory of a pde geometry data structure.
Definition: n_geom.c:50
double * area
Definition: N_pde.h:106
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition: proj3.c:138
void Rast3d_region_to_cell_head(RASTER3D_Region *, struct Cell_head *)
Returns in region2d the 2d portion of region3d.
Definition: region.c:48
N_geom_data * N_init_geom_data_3d(RASTER3D_Region *region3d, N_geom_data *geodata)
Initiate a pde geometry data structure with a 3d region.
Definition: n_geom.c:73
double dx
Definition: N_pde.h:109
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:448
double G_area_of_cell_at_row(int)
Cell area in specified row.
Definition: gis/area.c:87
int rows
Number of rows for 2D data.
Definition: gis.h:427
int G_debug(int, const char *,...) __attribute__((format(printf