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