GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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 * */
73N_geom_data *N_init_geom_data_3d(RASTER3D_Region *region3d,
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 * */
117N_geom_data *N_init_geom_data_2d(struct Cell_head *region, N_geom_data *geodata)
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}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
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:66
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:29
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_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:117
void N_free_geom_data(N_geom_data *geom)
Release memory of a pde geometry data structure.
Definition n_geom.c:49
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition proj3.c:146
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