GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 managment 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 
19 #include <math.h>
20 
21 #include "grass/N_pde.h"
22 #include "grass/glocale.h"
23 
24 /* ******************** 2D ARRAY FUNCTIONS *********************** */
25 
47 {
48  int map; /*The rastermap */
49  int x, y, cols, rows, type;
50  void *rast;
51  void *ptr;
52  struct Cell_head region;
53  N_array_2d *data = array;
54 
55  if (NULL == G_find_cell2(name, ""))
56  G_fatal_error(_("Raster map <%s> not found"), name);
57 
58  /* Get the active region */
59  G_get_set_window(&region);
60 
61  /*set the rows and cols */
62  rows = region.rows;
63  cols = region.cols;
64 
65  /*open the raster map */
66  map = G_open_cell_old(name, G_find_cell2(name, ""));
67  if (map < 0)
68  G_fatal_error(_("Unable to open raster map <%s>"), name);
69 
70  type = G_get_raster_map_type(map);
71 
72  /*if the array is NULL create a new one with the data type of the raster map */
73  /*the offset is 0 by default */
74  if (data == NULL) {
75  if (type == DCELL_TYPE) {
76  data = N_alloc_array_2d(cols, rows, 0, DCELL_TYPE);
77  }
78  if (type == FCELL_TYPE) {
79  data = N_alloc_array_2d(cols, rows, 0, FCELL_TYPE);
80  }
81  if (type == CELL_TYPE) {
82  data = N_alloc_array_2d(cols, rows, 0, CELL_TYPE);
83  }
84  }
85  else {
86  /*Check the array sizes */
87  if (data->cols != cols)
89  ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
90  if (data->rows != rows)
92  ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
93  }
94 
95  rast = G_allocate_raster_buf(type);
96 
97  G_message(_("Reading raster map <%s> into memory"), name);
98 
99  for (y = 0; y < rows; y++) {
100  G_percent(y, rows - 1, 10);
101 
102  if (!G_get_raster_row(map, rast, y, type)) {
103  G_close_cell(map);
104  G_fatal_error(_("Could not get raster row"));
105  }
106 
107  for (x = 0, ptr = rast; x < cols;
108  x++, ptr = G_incr_void_ptr(ptr, G_raster_size(type))) {
109  if (type == CELL_TYPE) {
110  if (G_is_c_null_value(ptr)) {
111  N_put_array_2d_value_null(data, x, y);
112  }
113  else {
114  if (data->type == CELL_TYPE)
115  N_put_array_2d_c_value(data, x, y,
116  (CELL) * (CELL *) ptr);
117  if (data->type == FCELL_TYPE)
118  N_put_array_2d_f_value(data, x, y,
119  (FCELL) * (CELL *) ptr);
120  if (data->type == DCELL_TYPE)
121  N_put_array_2d_d_value(data, x, y,
122  (DCELL) * (CELL *) ptr);
123  }
124  }
125  if (type == FCELL_TYPE) {
126  if (G_is_f_null_value(ptr)) {
127  N_put_array_2d_value_null(data, x, y);
128  }
129  else {
130  if (data->type == CELL_TYPE)
131  N_put_array_2d_c_value(data, x, y,
132  (CELL) * (FCELL *) ptr);
133  if (data->type == FCELL_TYPE)
134  N_put_array_2d_f_value(data, x, y,
135  (FCELL) * (FCELL *) ptr);
136  if (data->type == DCELL_TYPE)
137  N_put_array_2d_d_value(data, x, y,
138  (DCELL) * (FCELL *) ptr);
139  }
140  }
141  if (type == DCELL_TYPE) {
142  if (G_is_d_null_value(ptr)) {
143  N_put_array_2d_value_null(data, x, y);
144  }
145  else {
146  if (data->type == CELL_TYPE)
147  N_put_array_2d_c_value(data, x, y,
148  (CELL) * (DCELL *) ptr);
149  if (data->type == FCELL_TYPE)
150  N_put_array_2d_f_value(data, x, y,
151  (FCELL) * (DCELL *) ptr);
152  if (data->type == DCELL_TYPE)
153  N_put_array_2d_d_value(data, x, y,
154  (DCELL) * (DCELL *) ptr);
155  }
156  }
157  }
158  }
159 
160  /* Close file */
161  if (G_close_cell(map) < 0)
162  G_fatal_error(_("Unable to close input map"));
163 
164  return data;
165 }
166 
182 {
183  int map; /*The rastermap */
184  int x, y, cols, rows, count, type;
185  CELL *rast = NULL;
186  FCELL *frast = NULL;
187  DCELL *drast = NULL;
188  struct Cell_head region;
189 
190  if (!array)
191  G_fatal_error(_("N_array_2d * array is empty"));
192 
193  /* Get the current region */
194  G_get_set_window(&region);
195 
196  rows = region.rows;
197  cols = region.cols;
198  type = array->type;
199 
200  /*Open the new map */
201  map = G_open_raster_new(name, type);
202  if (map < 0)
203  G_fatal_error(_("Unable to create raster map <%s>"), name);
204 
205  if (type == CELL_TYPE)
206  rast = G_allocate_raster_buf(type);
207  if (type == FCELL_TYPE)
208  frast = G_allocate_raster_buf(type);
209  if (type == DCELL_TYPE)
210  drast = G_allocate_raster_buf(type);
211 
212  G_message(_("Write 2d array to raster map <%s>"), name);
213 
214  count = 0;
215  for (y = 0; y < rows; y++) {
216  G_percent(y, rows - 1, 10);
217  for (x = 0; x < cols; x++) {
218  if (type == CELL_TYPE)
219  rast[x] = N_get_array_2d_c_value(array, x, y);
220  if (type == FCELL_TYPE)
221  frast[x] = N_get_array_2d_f_value(array, x, y);
222  if (type == DCELL_TYPE)
223  drast[x] = N_get_array_2d_d_value(array, x, y);
224  }
225  if (type == CELL_TYPE)
226  if (!G_put_c_raster_row(map, rast)) {
227  G_unopen_cell(map); /*unopen the new raster map */
228  G_fatal_error(_("Unable to write raster row %i"), y);
229  }
230  if (type == FCELL_TYPE)
231  if (!G_put_f_raster_row(map, frast)) {
232  G_unopen_cell(map); /*unopen the new raster map */
233  G_fatal_error(_("Unable to write raster row %i"), y);
234  }
235  if (type == DCELL_TYPE)
236  if (!G_put_d_raster_row(map, drast)) {
237  G_unopen_cell(map); /*unopen the new raster map */
238  G_fatal_error(_("Unable to write raster row %i"), y);
239  }
240  }
241 
242  /* Close file */
243  if (G_close_cell(map) < 0)
244  G_fatal_error(_("Unable to close input map"));
245 
246  return;
247 
248 }
249 
250 
251 /* ******************** 3D ARRAY FUNCTIONS *********************** */
252 
277  int mask)
278 {
279  void *map = NULL; /*The 3D Rastermap */
280  int changemask = 0;
281  int x, y, z, cols, rows, depths, type;
282  double d1 = 0, f1 = 0;
283  N_array_3d *data = array;
284  G3D_Region region;
285 
286 
287  /*get the current region */
288  G3d_getWindow(&region);
289 
290  cols = region.cols;
291  rows = region.rows;
292  depths = region.depths;
293 
294 
295  if (NULL == G_find_grid3(name, ""))
296  G3d_fatalError(_("3D raster map <%s> not found"), name);
297 
298  /*Open all maps with default region */
299  map =
300  G3d_openCellOld(name, G_find_grid3(name, ""), G3D_DEFAULT_WINDOW,
301  G3D_TILE_SAME_AS_FILE, G3D_USE_CACHE_DEFAULT);
302 
303  if (map == NULL)
304  G3d_fatalError(_("Unable to open 3D raster map <%s>"), name);
305 
306  type = G3d_tileTypeMap(map);
307 
308  /*if the array is NULL create a new one with the data type of the volume map */
309  /*the offset is 0 by default */
310  if (data == NULL) {
311  if (type == FCELL_TYPE) {
312  data = N_alloc_array_3d(cols, rows, depths, 0, FCELL_TYPE);
313  }
314  if (type == DCELL_TYPE) {
315  data = N_alloc_array_3d(cols, rows, depths, 0, DCELL_TYPE);
316  }
317  }
318  else {
319  /*Check the array sizes */
320  if (data->cols != cols)
322  ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
323  if (data->rows != rows)
325  ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
326  if (data->depths != depths)
328  ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
329  }
330 
331 
332  G_message(_("Read g3d map <%s> into the memory"), name);
333 
334  /*if requested set the Mask on */
335  if (mask) {
336  if (G3d_maskFileExists()) {
337  changemask = 0;
338  if (G3d_maskIsOff(map)) {
339  G3d_maskOn(map);
340  changemask = 1;
341  }
342  }
343  }
344 
345  for (z = 0; z < depths; z++) { /*From the bottom to the top */
346  G_percent(z, depths - 1, 10);
347  for (y = 0; y < rows; y++) {
348  for (x = 0; x < cols; x++) {
349  if (type == FCELL_TYPE) {
350  G3d_getValue(map, x, y, z, &f1, type);
351  if (G_is_f_null_value((void *)&f1)) {
352  N_put_array_3d_value_null(data, x, y, z);
353  }
354  else {
355  if (data->type == FCELL_TYPE)
356  N_put_array_3d_f_value(data, x, y, z, f1);
357  if (data->type == DCELL_TYPE)
358  N_put_array_3d_d_value(data, x, y, z, (double)f1);
359  }
360  }
361  else {
362  G3d_getValue(map, x, y, z, &d1, type);
363  if (G_is_d_null_value((void *)&d1)) {
364  N_put_array_3d_value_null(data, x, y, z);
365  }
366  else {
367  if (data->type == FCELL_TYPE)
368  N_put_array_3d_f_value(data, x, y, z, (float)d1);
369  if (data->type == DCELL_TYPE)
370  N_put_array_3d_d_value(data, x, y, z, d1);
371  }
372 
373  }
374  }
375  }
376  }
377 
378  /*We set the Mask off, if it was off before */
379  if (mask) {
380  if (G3d_maskFileExists())
381  if (G3d_maskIsOn(map) && changemask)
382  G3d_maskOff(map);
383  }
384 
385  /* Close files and exit */
386  if (!G3d_closeCell(map))
387  G3d_fatalError(map, NULL, 0, _("Error closing g3d file"));
388 
389  return data;
390 }
391 
408 void N_write_array_3d_to_rast3d(N_array_3d * array, char *name, int mask)
409 {
410  void *map = NULL; /*The 3D Rastermap */
411  int changemask = 0;
412  int x, y, z, cols, rows, depths, count, type;
413  double d1 = 0.0, f1 = 0.0;
414  N_array_3d *data = array;
415  G3D_Region region;
416 
417  /*get the current region */
418  G3d_getWindow(&region);
419 
420 
421  cols = region.cols;
422  rows = region.rows;
423  depths = region.depths;
424  type = data->type;
425 
426  /*Check the array sizes */
427  if (data->cols != cols)
429  ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
430  if (data->rows != rows)
432  ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
433  if (data->depths != depths)
435  ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
436 
437  /*Open the new map */
438  if (type == DCELL_TYPE)
439  map =
440  G3d_openCellNew(name, DCELL_TYPE, G3D_USE_CACHE_DEFAULT, &region);
441  else if (type == FCELL_TYPE)
442  map =
443  G3d_openCellNew(name, FCELL_TYPE, G3D_USE_CACHE_DEFAULT, &region);
444 
445  if (map == NULL)
446  G3d_fatalError(_("Error opening g3d map <%s>"), name);
447 
448  G_message(_("Write 3d array to g3d map <%s>"), name);
449 
450  /*if requested set the Mask on */
451  if (mask) {
452  if (G3d_maskFileExists()) {
453  changemask = 0;
454  if (G3d_maskIsOff(map)) {
455  G3d_maskOn(map);
456  changemask = 1;
457  }
458  }
459  }
460 
461  count = 0;
462  for (z = 0; z < depths; z++) { /*From the bottom to the top */
463  G_percent(z, depths - 1, 10);
464  for (y = 0; y < rows; y++) {
465  for (x = 0; x < cols; x++) {
466  if (type == FCELL_TYPE) {
467  f1 = N_get_array_3d_f_value(data, x, y, z);
468  G3d_putFloat(map, x, y, z, f1);
469  }
470  else if (type == DCELL_TYPE) {
471  d1 = N_get_array_3d_d_value(data, x, y, z);
472  G3d_putDouble(map, x, y, z, d1);
473  }
474  }
475  }
476  }
477 
478  /*We set the Mask off, if it was off before */
479  if (mask) {
480  if (G3d_maskFileExists())
481  if (G3d_maskIsOn(map) && changemask)
482  G3d_maskOff(map);
483  }
484 
485  /* Close files and exit */
486  if (!G3d_closeCell(map))
487  G3d_fatalError(map, NULL, 0, _("Error closing g3d file"));
488 
489  return;
490 }