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