GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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 * */
73N_geom_data *N_init_geom_data_3d(RASTER3D_Region * region3d, N_geom_data * geodata)
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 * */
114N_geom_data *N_init_geom_data_2d(struct Cell_head * region,
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}
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
double G_area_of_cell_at_row(int row)
Cell area in specified row.
Definition: area.c:87
int G_begin_cell_area_calculations(void)
Begin cell area calculations.
Definition: area.c:46
#define NULL
Definition: ccmath.h:32
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
void G_get_set_window(struct Cell_head *window)
Get the current working window (region)
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
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
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_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
void N_free_geom_data(N_geom_data *geom)
Release memory of a pde geometry data structure.
Definition: n_geom.c:50
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition: proj3.c:138
Geometric information about the structured grid.
Definition: N_pde.h:104
double dx
Definition: N_pde.h:109
int depths
Definition: N_pde.h:115
double dy
Definition: N_pde.h:110
double * area
Definition: N_pde.h:106
double dz
Definition: N_pde.h:111
int dim
Definition: N_pde.h:107
int rows
Definition: N_pde.h:116
int planimetric
Definition: N_pde.h:105
int cols
Definition: N_pde.h:117
double Az
Definition: N_pde.h:113