GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
n_arrays_io.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: IO array management functions
9* part of the gpde library
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#include <math.h>
19
20#include <grass/N_pde.h>
21#include <grass/raster.h>
22#include <grass/glocale.h>
23
24
25/* ******************** 2D ARRAY FUNCTIONS *********************** */
26
27/*!
28 * \brief Read a raster map into a N_array_2d structure
29 *
30 * The raster map will be opened in the current region settings.
31 * If no N_array_2d structure is provided (NULL pointer), a new structure will be
32 * allocated with the same data type as the raster map and the size of the current region.
33 * The array offset will be set to 0.
34 * <br><br>
35 * If a N_array_2d structure is provided, the values from the raster map are
36 * casted to the N_array_2d type. The array must have the same size
37 * as the current region.
38 * <br><br>
39 * The new created or the provided array are returned.
40 * If the reading of the raster map fails, G_fatal_error() will
41 * be invoked.
42 *
43 * \param name * char - the name of an existing raster map
44 * \param array * N_array_2d - an existing array or NULL
45 * \return N_array_2d * - the existing or new allocated array
46 * */
48{
49 int map; /*The rastermap */
50 int x, y, cols, rows, type;
51 void *rast;
52 void *ptr;
53 struct Cell_head region;
54 N_array_2d *data = array;
55
56 /* Get the active region */
57 G_get_set_window(&region);
58
59 /*set the rows and cols */
60 rows = region.rows;
61 cols = region.cols;
62
63 /*open the raster map */
64 map = Rast_open_old(name, "");
65
66 type = Rast_get_map_type(map);
67
68 /*if the array is NULL create a new one with the data type of the raster map */
69 /*the offset is 0 by default */
70 if (data == NULL) {
71 if (type == DCELL_TYPE) {
72 data = N_alloc_array_2d(cols, rows, 0, DCELL_TYPE);
73 }
74 if (type == FCELL_TYPE) {
75 data = N_alloc_array_2d(cols, rows, 0, FCELL_TYPE);
76 }
77 if (type == CELL_TYPE) {
78 data = N_alloc_array_2d(cols, rows, 0, CELL_TYPE);
79 }
80 }
81 else {
82 /*Check the array sizes */
83 if (data->cols != cols)
85 ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
86 if (data->rows != rows)
88 ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
89 }
90
91 rast = Rast_allocate_buf(type);
92
93 G_message(_("Reading raster map <%s> into memory"), name);
94
95 for (y = 0; y < rows; y++) {
96 G_percent(y, rows - 1, 10);
97
98 Rast_get_row(map, rast, y, type);
99
100 for (x = 0, ptr = rast; x < cols;
101 x++, ptr = G_incr_void_ptr(ptr, Rast_cell_size(type))) {
102 if (type == CELL_TYPE) {
103 if (Rast_is_c_null_value(ptr)) {
105 }
106 else {
107 if (data->type == CELL_TYPE)
108 N_put_array_2d_c_value(data, x, y,
109 (CELL) * (CELL *) ptr);
110 if (data->type == FCELL_TYPE)
111 N_put_array_2d_f_value(data, x, y,
112 (FCELL) * (CELL *) ptr);
113 if (data->type == DCELL_TYPE)
114 N_put_array_2d_d_value(data, x, y,
115 (DCELL) * (CELL *) ptr);
116 }
117 }
118 if (type == FCELL_TYPE) {
119 if (Rast_is_f_null_value(ptr)) {
121 }
122 else {
123 if (data->type == CELL_TYPE)
124 N_put_array_2d_c_value(data, x, y,
125 (CELL) * (FCELL *) ptr);
126 if (data->type == FCELL_TYPE)
127 N_put_array_2d_f_value(data, x, y,
128 (FCELL) * (FCELL *) ptr);
129 if (data->type == DCELL_TYPE)
130 N_put_array_2d_d_value(data, x, y,
131 (DCELL) * (FCELL *) ptr);
132 }
133 }
134 if (type == DCELL_TYPE) {
135 if (Rast_is_d_null_value(ptr)) {
137 }
138 else {
139 if (data->type == CELL_TYPE)
140 N_put_array_2d_c_value(data, x, y,
141 (CELL) * (DCELL *) ptr);
142 if (data->type == FCELL_TYPE)
143 N_put_array_2d_f_value(data, x, y,
144 (FCELL) * (DCELL *) ptr);
145 if (data->type == DCELL_TYPE)
146 N_put_array_2d_d_value(data, x, y,
147 (DCELL) * (DCELL *) ptr);
148 }
149 }
150 }
151 }
152
153 /* Close file */
154 Rast_close(map);
155
156 return data;
157}
158
159/*!
160 * \brief Write a N_array_2d struct to a raster map
161 *
162 * A new raster map is created with the same type as the N_array_2d.
163 * The current region is used to open the raster map.
164 * The N_array_2d must have the same size as the current region.
165 If the writing of the raster map fails, G_fatal_error() will
166 * be invoked.
167
168 * \param array N_array_2d *
169 * \param name char * - the name of the raster map
170 * \return void
171 *
172 * */
174{
175 int map; /*The rastermap */
176 int x, y, cols, rows, count, type;
177 CELL *rast = NULL;
178 FCELL *frast = NULL;
179 DCELL *drast = NULL;
180 struct Cell_head region;
181
182 if (!array)
183 G_fatal_error(_("N_array_2d * array is empty"));
184
185 /* Get the current region */
186 G_get_set_window(&region);
187
188 rows = region.rows;
189 cols = region.cols;
190 type = array->type;
191
192 /*Open the new map */
193 map = Rast_open_new(name, type);
194
195 if (type == CELL_TYPE)
196 rast = Rast_allocate_buf(type);
197 if (type == FCELL_TYPE)
198 frast = Rast_allocate_buf(type);
199 if (type == DCELL_TYPE)
200 drast = Rast_allocate_buf(type);
201
202 G_message(_("Write 2d array to raster map <%s>"), name);
203
204 count = 0;
205 for (y = 0; y < rows; y++) {
206 G_percent(y, rows - 1, 10);
207 for (x = 0; x < cols; x++) {
208 if (type == CELL_TYPE)
209 rast[x] = N_get_array_2d_c_value(array, x, y);
210 if (type == FCELL_TYPE)
211 frast[x] = N_get_array_2d_f_value(array, x, y);
212 if (type == DCELL_TYPE)
213 drast[x] = N_get_array_2d_d_value(array, x, y);
214 }
215 if (type == CELL_TYPE)
216 Rast_put_c_row(map, rast);
217 if (type == FCELL_TYPE)
218 Rast_put_f_row(map, frast);
219 if (type == DCELL_TYPE)
220 Rast_put_d_row(map, drast);
221 }
222
223 /* Close file */
224 Rast_close(map);
225}
226
227
228/* ******************** 3D ARRAY FUNCTIONS *********************** */
229
230/*!
231 * \brief Read a volume map into a N_array_3d structure
232 *
233 * The volume map is opened in the current region settings.
234 * If no N_array_3d structure is provided (NULL pointer), a new structure will be
235 * allocated with the same data type as the volume map and the size of the current region.
236 * The array offset will be set to 0.
237 * <br><br>
238 *
239 * If a N_array_3d structure is provided, the values from the volume map are
240 * casted to the N_array_3d type. The array must have the same size
241 * as the current region.
242 * <br><br>
243 *
244 * The new created or the provided array is returned.
245 * If the reading of the volume map fails, Rast3d_fatal_error() will
246 * be invoked.
247 *
248 * \param name * char - the name of an existing volume map
249 * \param array * N_array_3d - an existing array or NULL
250 * \param mask int - 0 = false, 1 = ture : if a mask is presenent, use it with the input volume map
251 * \return N_array_3d * - the existing or new allocated array
252 * */
254 int mask)
255{
256 void *map = NULL; /*The 3D Rastermap */
257 int changemask = 0;
258 int x, y, z, cols, rows, depths, type;
259 double d1 = 0, f1 = 0;
260 N_array_3d *data = array;
261 RASTER3D_Region region;
262
263
264 /*get the current region */
265 Rast3d_get_window(&region);
266
267 cols = region.cols;
268 rows = region.rows;
269 depths = region.depths;
270
271
272 if (NULL == G_find_raster3d(name, ""))
273 Rast3d_fatal_error(_("3D raster map <%s> not found"), name);
274
275 /*Open all maps with default region */
276 map =
277 Rast3d_open_cell_old(name, G_find_raster3d(name, ""), RASTER3D_DEFAULT_WINDOW,
278 RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
279
280 if (map == NULL)
281 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"), name);
282
283 type = Rast3d_tile_type_map(map);
284
285 /*if the array is NULL create a new one with the data type of the volume map */
286 /*the offset is 0 by default */
287 if (data == NULL) {
288 if (type == FCELL_TYPE) {
289 data = N_alloc_array_3d(cols, rows, depths, 0, FCELL_TYPE);
290 }
291 if (type == DCELL_TYPE) {
292 data = N_alloc_array_3d(cols, rows, depths, 0, DCELL_TYPE);
293 }
294 }
295 else {
296 /*Check the array sizes */
297 if (data->cols != cols)
299 ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
300 if (data->rows != rows)
302 ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
303 if (data->depths != depths)
305 ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
306 }
307
308
309 G_message(_("Read g3d map <%s> into the memory"), name);
310
311 /*if requested set the Mask on */
312 if (mask) {
313 if (Rast3d_mask_file_exists()) {
314 changemask = 0;
315 if (Rast3d_mask_is_off(map)) {
316 Rast3d_mask_on(map);
317 changemask = 1;
318 }
319 }
320 }
321
322 for (z = 0; z < depths; z++) { /*From the bottom to the top */
323 G_percent(z, depths - 1, 10);
324 for (y = 0; y < rows; y++) {
325 for (x = 0; x < cols; x++) {
326 if (type == FCELL_TYPE) {
327 Rast3d_get_value(map, x, y, z, &f1, type);
328 if (Rast_is_f_null_value((void *)&f1)) {
329 N_put_array_3d_value_null(data, x, y, z);
330 }
331 else {
332 if (data->type == FCELL_TYPE)
333 N_put_array_3d_f_value(data, x, y, z, f1);
334 if (data->type == DCELL_TYPE)
335 N_put_array_3d_d_value(data, x, y, z, (double)f1);
336 }
337 }
338 else {
339 Rast3d_get_value(map, x, y, z, &d1, type);
340 if (Rast_is_d_null_value((void *)&d1)) {
341 N_put_array_3d_value_null(data, x, y, z);
342 }
343 else {
344 if (data->type == FCELL_TYPE)
345 N_put_array_3d_f_value(data, x, y, z, (float)d1);
346 if (data->type == DCELL_TYPE)
347 N_put_array_3d_d_value(data, x, y, z, d1);
348 }
349
350 }
351 }
352 }
353 }
354
355 /*We set the Mask off, if it was off before */
356 if (mask) {
357 if (Rast3d_mask_file_exists())
358 if (Rast3d_mask_is_on(map) && changemask)
359 Rast3d_mask_off(map);
360 }
361
362 /* Close files and exit */
363 if (!Rast3d_close(map))
364 Rast3d_fatal_error(map, NULL, 0, _("Error closing g3d file"));
365
366 return data;
367}
368
369/*!
370 * \brief Write a N_array_3d struct to a volume map
371 *
372 * A new volume map is created with the same type as the N_array_3d.
373 * The current region is used to open the volume map.
374 * The N_array_3d must have the same size as the current region.
375 * If the writing of the volume map fails, Rast3d_fatal_error() will
376 * be invoked.
377 *
378 *
379 * \param array N_array_3d *
380 * \param name char * - the name of the volume map
381 * \param mask int - 1 = use a 3d mask, 0 do not use a 3d mask
382 * \return void
383 *
384 * */
385void N_write_array_3d_to_rast3d(N_array_3d * array, char *name, int mask)
386{
387 void *map = NULL; /*The 3D Rastermap */
388 int changemask = 0;
389 int x, y, z, cols, rows, depths, count, type;
390 double d1 = 0.0, f1 = 0.0;
391 N_array_3d *data = array;
392 RASTER3D_Region region;
393
394 /*get the current region */
395 Rast3d_get_window(&region);
396
397 cols = region.cols;
398 rows = region.rows;
399 depths = region.depths;
400 type = data->type;
401
402 /*Check the array sizes */
403 if (data->cols != cols)
405 ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
406 if (data->rows != rows)
408 ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
409 if (data->depths != depths)
411 ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
412
413 /*Open the new map */
414 if (type == DCELL_TYPE)
415 map = Rast3d_open_new_opt_tile_size(name, RASTER3D_USE_CACHE_XY, &region, DCELL_TYPE, 32);
416 else if (type == FCELL_TYPE)
417 map = Rast3d_open_new_opt_tile_size(name, RASTER3D_USE_CACHE_XY, &region, FCELL_TYPE, 32);
418
419 if (map == NULL)
420 Rast3d_fatal_error(_("Error opening g3d map <%s>"), name);
421
422 G_message(_("Write 3d array to g3d map <%s>"), name);
423
424 /*if requested set the Mask on */
425 if (mask) {
426 if (Rast3d_mask_file_exists()) {
427 changemask = 0;
428 if (Rast3d_mask_is_off(map)) {
429 Rast3d_mask_on(map);
430 changemask = 1;
431 }
432 }
433 }
434
435 count = 0;
436 for (z = 0; z < depths; z++) { /*From the bottom to the top */
437 G_percent(z, depths - 1, 10);
438 for (y = 0; y < rows; y++) {
439 for (x = 0; x < cols; x++) {
440 if (type == FCELL_TYPE) {
441 f1 = N_get_array_3d_f_value(data, x, y, z);
442 Rast3d_put_float(map, x, y, z, f1);
443 }
444 else if (type == DCELL_TYPE) {
445 d1 = N_get_array_3d_d_value(data, x, y, z);
446 Rast3d_put_double(map, x, y, z, d1);
447 }
448 }
449 }
450 }
451
452 /*We set the Mask off, if it was off before */
453 if (mask) {
454 if (Rast3d_mask_file_exists())
455 if (Rast3d_mask_is_on(map) && changemask)
456 Rast3d_mask_off(map);
457 }
458
459 /* Flush all tile */
460 if (!Rast3d_flush_all_tiles(map))
461 Rast3d_fatal_error("Error flushing tiles with Rast3d_flush_all_tiles");
462 /* Close files and exit */
463 if (!Rast3d_close(map))
464 Rast3d_fatal_error(map, NULL, 0, _("Error closing g3d file"));
465
466 return;
467}
void * G_incr_void_ptr(const void *ptr, size_t size)
Advance void pointer.
Definition: alloc.c:186
#define NULL
Definition: ccmath.h:32
const char * G_find_raster3d(const char *name, const char *mapset)
Search for a 3D raster map in current search path or in a specified mapset.
Definition: find_rast3d.c:28
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
void G_get_set_window(struct Cell_head *window)
Get the current working window (region)
int count
FCELL N_get_array_2d_f_value(N_array_2d *data, int col, int row)
Returns the value of type FCELL at position col, row.
Definition: n_arrays.c:346
void N_put_array_3d_value_null(N_array_3d *data, int col, int row, int depth)
This function writes a null value to the N_array_3d data at position col, row, depth.
Definition: n_arrays.c:1076
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition: n_arrays.c:726
CELL N_get_array_2d_c_value(N_array_2d *data, int col, int row)
Returns the value of type CELL at position col, row.
Definition: n_arrays.c:314
void N_put_array_2d_f_value(N_array_2d *data, int col, int row, FCELL value)
Writes a FCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:554
float N_get_array_3d_f_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: n_arrays.c:961
void N_put_array_3d_f_value(N_array_3d *data, int col, int row, int depth, float value)
This function writes a float value to the N_array_3d data at position col, row, depth.
Definition: n_arrays.c:1148
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: n_arrays.c:378
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition: n_arrays.c:1175
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition: n_arrays.c:72
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: n_arrays.c:990
void N_put_array_2d_value_null(N_array_2d *data, int col, int row)
Writes the null value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:458
void N_put_array_2d_c_value(N_array_2d *data, int col, int row, CELL value)
Writes a CELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:524
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:584
void N_write_array_3d_to_rast3d(N_array_3d *array, char *name, int mask)
Write a N_array_3d struct to a volume map.
Definition: n_arrays_io.c:385
N_array_2d * N_read_rast_to_array_2d(char *name, N_array_2d *array)
Read a raster map into a N_array_2d structure.
Definition: n_arrays_io.c:47
N_array_3d * N_read_rast3d_to_array_3d(char *name, N_array_3d *array, int mask)
Read a volume map into a N_array_3d structure.
Definition: n_arrays_io.c:253
void N_write_array_2d_to_rast(N_array_2d *array, char *name)
Write a N_array_2d struct to a raster map.
Definition: n_arrays_io.c:173
const char * name
Definition: named_colr.c:7
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:62
int type
Definition: N_pde.h:133
int type
Definition: N_pde.h:167
int rows
Definition: N_pde.h:168
int depths
Definition: N_pde.h:168
int cols
Definition: N_pde.h:168
#define x