GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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)) {
104  N_put_array_2d_value_null(data, x, y);
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)) {
120  N_put_array_2d_value_null(data, x, y);
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)) {
136  N_put_array_2d_value_null(data, x, y);
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 =
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) {
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  * */
385 void 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)
416  else if (type == FCELL_TYPE)
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) {
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 }
#define CELL_TYPE
Definition: raster.h:11
#define RASTER3D_USE_CACHE_DEFAULT
Definition: raster3d.h:20
int Rast3d_mask_file_exists(void)
Returns 1 if the 3d mask file exists.
Definition: mask.c:64
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 void void void G_fatal_error(const char *,...) __attribute__((format(printf
void Rast_put_d_row(int, const DCELL *)
Writes the next row for dcell file (DCELL version)
void * Rast3d_open_cell_old(const char *, const char *, RASTER3D_Region *, int, int)
Opens existing g3d-file name in mapset. Tiles are stored in memory with type which must be any of FCE...
Definition: raster3d/open.c:77
2D/3D raster map header (used also for region)
Definition: gis.h:423
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
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:414
int Rast3d_put_float(RASTER3D_Map *, int, int, int, float)
Is equivalent to Rast3d_put_value (map, x, y, z, &value, FCELL_TYPE).
Definition: putvalue.c:21
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
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
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
double DCELL
Definition: gis.h:614
int depths
Definition: N_pde.h:168
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
#define RASTER3D_DEFAULT_WINDOW
Definition: raster3d.h:29
int count
int Rast3d_close(RASTER3D_Map *)
Close 3D raster map files.
#define G_incr_void_ptr(ptr, size)
Definition: defs/gis.h:100
const char * G_find_raster3d(const char *, const char *)
Search for a 3D raster map in current search path or in a specified mapset.
Definition: find_rast3d.c:28
#define Rast_is_f_null_value(fcellVal)
Definition: defs/raster.h:412
#define NULL
Definition: ccmath.h:32
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
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
#define x
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_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
void Rast3d_fatal_error(const char *,...) __attribute__((format(printf
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
void G_message(const char *,...) __attribute__((format(printf
#define DCELL_TYPE
Definition: raster.h:13
void Rast3d_mask_on(RASTER3D_Map *)
Turns on the mask for map. Do not invoke this function after the first tile has been read since the r...
Definition: mask.c:332
int Rast_open_new(const char *, RASTER_MAP_TYPE)
Opens a new raster map.
Definition: raster/open.c:997
size_t Rast_cell_size(RASTER_MAP_TYPE)
Returns size of a raster cell in bytes.
Definition: alloc_cell.c:39
void * Rast3d_open_new_opt_tile_size(const char *, int, RASTER3D_Region *, int, int)
Opens new g3d-file with name in the current mapset. This method tries to compute optimal tile size ba...
Definition: open2.c:98
int cols
Definition: N_pde.h:168
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
void Rast_put_f_row(int, const FCELL *)
Writes the next row for fcell file (FCELL version)
void Rast_get_row(int, void *, int, RASTER_MAP_TYPE)
Get raster row.
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
int Rast3d_put_double(RASTER3D_Map *, int, int, int, double)
Is equivalent to Rast3d_put_value (map, x, y, z, &value, DCELL_TYPE).
Definition: putvalue.c:57
void G_get_set_window(struct Cell_head *)
Get the current working window (region)
#define RASTER3D_USE_CACHE_XY
Definition: raster3d.h:24
int Rast_open_old(const char *, const char *)
Open an existing integer raster map (cell)
Definition: raster/open.c:112
int Rast3d_mask_is_on(RASTER3D_Map *)
Returns 1 if the mask for map is turned on. Returns 0 otherwise.
Definition: mask.c:365
int depths
number of depths for 3D data
Definition: gis.h:446
float FCELL
Definition: gis.h:615
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:62
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
int cols
Number of columns for 2D data.
Definition: gis.h:442
int type
Definition: N_pde.h:133
void Rast_put_c_row(int, const CELL *)
Writes the next row for cell file (CELL version)
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
#define RASTER3D_TILE_SAME_AS_FILE
Definition: raster3d.h:12
void * Rast_allocate_buf(RASTER_MAP_TYPE)
Allocate memory for a raster map of given type.
Definition: alloc_cell.c:55
int CELL
Definition: gis.h:613
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_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
RASTER_MAP_TYPE Rast_get_map_type(int)
Determine raster type from descriptor.
Definition: raster/open.c:918
#define _(str)
Definition: glocale.h:10
#define FCELL_TYPE
Definition: raster.h:12
int Rast3d_flush_all_tiles(RASTER3D_Map *)
Definition: cache.c:280
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
const char * name
Definition: named_colr.c:7
void Rast3d_get_value(RASTER3D_Map *, int, int, int, void *, int)
Returns in *value the resampled cell-value of the cell with window-coordinate (x, y...
Definition: getvalue.c:22
int Rast3d_tile_type_map(RASTER3D_Map *)
Returns the type in which tiles of map are stored in memory.
Definition: headerinfo.c:161
int rows
Number of rows for 2D data.
Definition: gis.h:438
int type
Definition: N_pde.h:167
#define Rast_is_c_null_value(cellVal)
Definition: defs/raster.h:410
int rows
Definition: N_pde.h:168
void Rast_close(int)
Close a raster map.
Definition: raster/close.c:99
void Rast3d_get_window(RASTER3D_Region *)
Stores the current default window in window.
void Rast3d_mask_off(RASTER3D_Map *)
Turns off the mask for map. This is the default. Do not invoke this function after the first tile has...
Definition: mask.c:349
int Rast3d_mask_is_off(RASTER3D_Map *)
Returns 1 if the mask for map is turned off. Returns 0 otherwise.
Definition: mask.c:380