← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/gpde/n_arrays_io.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gpde/n_arrays_io.c
 
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
#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
 * */
 
47
N_array_2d *N_read_rast_to_array_2d(char *name, N_array_2d * array)
 
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)
 
84
            G_fatal_error
 
85
                ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
 
86
        if (data->rows != rows)
 
87
            G_fatal_error
 
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
 * */
 
173
void N_write_array_2d_to_rast(N_array_2d * array, char *name)
 
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
 * */
 
253
N_array_3d *N_read_rast3d_to_array_3d(char *name, N_array_3d * array,
 
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)
 
298
            G_fatal_error
 
299
                ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
 
300
        if (data->rows != rows)
 
301
            G_fatal_error
 
302
                ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
 
303
        if (data->depths != depths)
 
304
            G_fatal_error
 
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
 * */
 
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)
 
404
        G_fatal_error
 
405
            ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
 
406
    if (data->rows != rows)
 
407
        G_fatal_error
 
408
            ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
 
409
    if (data->depths != depths)
 
410
        G_fatal_error
 
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
}