GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-2bdcc05b3a
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)) {
102  N_put_array_2d_value_null(data, x, y);
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)) {
118  N_put_array_2d_value_null(data, x, y);
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)) {
134  N_put_array_2d_value_null(data, x, y);
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(
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) {
284  }
285  if (type == 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) {
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  * */
377 void 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)
408  &region, DCELL_TYPE, 32);
409  else if (type == FCELL_TYPE)
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) {
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 }
#define NULL
Definition: ccmath.h:32
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:61
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_get_set_window(struct Cell_head *)
Get the current working window (region)
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 G_incr_void_ptr(ptr, size)
Definition: defs/gis.h:81
void G_message(const char *,...) __attribute__((format(printf
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_flush_all_tiles(RASTER3D_Map *)
Definition: cache.c:283
int Rast3d_mask_is_off(RASTER3D_Map *)
Returns 1 if the mask for map is turned off. Returns 0 otherwise.
Definition: mask.c:374
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:18
int Rast3d_tile_type_map(RASTER3D_Map *)
Returns the type in which tiles of map are stored in memory.
Definition: headerinfo.c:156
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:80
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:345
int Rast3d_mask_file_exists(void)
Returns 1 if the 3d mask file exists.
Definition: mask.c:64
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:329
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:94
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:53
void Rast3d_get_window(RASTER3D_Region *)
Stores the current default window in window.
void Rast3d_fatal_error(const char *,...) __attribute__((format(printf
int Rast3d_close(RASTER3D_Map *)
Close 3D raster map files.
int Rast3d_mask_is_on(RASTER3D_Map *)
Returns 1 if the mask for map is turned on. Returns 0 otherwise.
Definition: mask.c:360
void * Rast_allocate_buf(RASTER_MAP_TYPE)
Allocate memory for a raster map of given type.
Definition: alloc_cell.c:54
#define Rast_is_f_null_value(fcellVal)
Definition: defs/raster.h:403
void Rast_close(int)
Close a raster map.
Definition: raster/close.c:99
int Rast_open_old(const char *, const char *)
Open an existing integer raster map (cell)
Definition: raster/open.c:112
void Rast_put_f_row(int, const FCELL *)
Writes the next row for fcell file (FCELL version)
int Rast_open_new(const char *, RASTER_MAP_TYPE)
Opens a new raster map.
Definition: raster/open.c:1013
size_t Rast_cell_size(RASTER_MAP_TYPE)
Returns size of a raster cell in bytes.
Definition: alloc_cell.c:38
void Rast_put_d_row(int, const DCELL *)
Writes the next row for dcell file (DCELL version)
void Rast_put_c_row(int, const CELL *)
Writes the next row for cell file (CELL version)
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:405
#define Rast_is_c_null_value(cellVal)
Definition: defs/raster.h:401
RASTER_MAP_TYPE Rast_get_map_type(int)
Determine raster type from descriptor.
Definition: raster/open.c:932
void Rast_get_row(int, void *, int, RASTER_MAP_TYPE)
Get raster row.
float FCELL
Definition: gis.h:629
double DCELL
Definition: gis.h:628
int CELL
Definition: gis.h:627
#define _(str)
Definition: glocale.h:10
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
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
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
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
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
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.
Definition: n_arrays_io.c:377
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:249
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
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:171
const char * name
Definition: named_colr.c:6
#define RASTER3D_DEFAULT_WINDOW
Definition: raster3d.h:29
#define RASTER3D_USE_CACHE_DEFAULT
Definition: raster3d.h:19
#define RASTER3D_TILE_SAME_AS_FILE
Definition: raster3d.h:11
#define RASTER3D_USE_CACHE_XY
Definition: raster3d.h:23
#define FCELL_TYPE
Definition: raster.h:12
#define DCELL_TYPE
Definition: raster.h:13
#define CELL_TYPE
Definition: raster.h:11
2D/3D raster map header (used also for region)
Definition: gis.h:439
int depths
number of depths for 3D data
Definition: gis.h:462
int rows
Number of rows for 2D data.
Definition: gis.h:454
int cols
Number of columns for 2D data.
Definition: gis.h:458
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