← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to lib/raster/open.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/raster/open.c
 
1
/*!
 
2
 * \file lib/raster/open.c
 
3
 * 
 
4
 * \brief Raster Library - Open raster file
 
5
 *
 
6
 * (C) 1999-2009 by the GRASS Development Team
 
7
 *
 
8
 * This program is free software under the GNU General Public
 
9
 * License (>=v2). Read the file COPYING that comes with GRASS
 
10
 * for details.
 
11
 *
 
12
 * \author USACERL and many others
 
13
 */
 
14
 
 
15
#include <rpc/types.h>
 
16
#include <rpc/xdr.h>
 
17
#include <unistd.h>
 
18
#include <string.h>
 
19
#include <sys/types.h>
 
20
#include <sys/stat.h>
 
21
#include <fcntl.h>
 
22
 
 
23
#include <grass/config.h>
 
24
#include <grass/gis.h>
 
25
#include <grass/raster.h>
 
26
#include <grass/glocale.h>
 
27
 
 
28
#include "R.h"
 
29
#define FORMAT_FILE "f_format"
 
30
#define NULL_FILE   "null"
 
31
 
 
32
static int new_fileinfo(void)
 
33
{
 
34
    int oldsize = R__.fileinfo_count;
 
35
    int newsize = oldsize;
 
36
    int i;
 
37
 
 
38
    for (i = 0; i < oldsize; i++)
 
39
        if (R__.fileinfo[i].open_mode <= 0) {
 
40
            memset(&R__.fileinfo[i], 0, sizeof(struct fileinfo));
 
41
            R__.fileinfo[i].open_mode = -1;
 
42
            return i;
 
43
        }
 
44
 
 
45
    if (newsize < 20)
 
46
        newsize += 20;
 
47
    else
 
48
        newsize *= 2;
 
49
 
 
50
    R__.fileinfo = G_realloc(R__.fileinfo, newsize * sizeof(struct fileinfo));
 
51
 
 
52
    /* Mark all cell files as closed */
 
53
    for (i = oldsize; i < newsize; i++) {
 
54
        memset(&R__.fileinfo[i], 0, sizeof(struct fileinfo));
 
55
        R__.fileinfo[i].open_mode = -1;
 
56
    }
 
57
 
 
58
    R__.fileinfo_count = newsize;
 
59
 
 
60
    return oldsize;
 
61
}
 
62
 
 
63
/*!
 
64
 * \brief Open raster file
 
65
 *
 
66
 * Arrange for the NULL-value bitmap to be read as well as the raster
 
67
 * map. If no NULL-value bitmap exists, arrange for the production of
 
68
 * NULL-values based on zeros in the raster map. If the map is
 
69
 * floating-point, arrange for quantization to integer for
 
70
 * Rast_get_c_row(), et. al., by reading the quantization rules
 
71
 * for the map using Rast_read_quant(). If the programmer wants to read
 
72
 * the floating point map using uing quant rules other than the ones
 
73
 * stored in map's quant file, he/she should call Rast_set_quant_rules()
 
74
 * after the call to Rast_open_old().
 
75
 *
 
76
 * \param name map name
 
77
 * \param open_mode mode
 
78
 * \param map_type map type (CELL, FCELL, DCELL)
 
79
 *
 
80
 * \return open file descriptor ( >= 0) if successful
 
81
 */
 
82
 
 
83
static int open_raster_new(const char *name, int open_mode,
 
84
                           RASTER_MAP_TYPE map_type);
 
85
 
 
86
/*!
 
87
   \brief Open an existing integer raster map (cell)
 
88
 
 
89
   Opens the existing cell file <i>name</i> in the <i>mapset</i> for
 
90
   reading by Rast_get_row() with mapping into the current window.
 
91
 
 
92
   This routine opens the raster map <i>name</i> in <i>mapset</i> for
 
93
   reading. A nonnegative file descriptor is returned if the open is
 
94
   successful. Otherwise a diagnostic message is printed and a negative
 
95
   value is returned. This routine does quite a bit of work. Since
 
96
   GRASS users expect that all raster maps will be resampled into the
 
97
   current region, the resampling index for the raster map is prepared
 
98
   by this routine after the file is opened. The resampling is based on
 
99
   the active module region (see also \ref The_Region}. Preparation
 
100
   required for reading the various raster file formats (see \ref
 
101
   Raster_File_Format for an explanation of the various raster file
 
102
   formats) is also done.
 
103
 
 
104
   Diagnostics: warning message printed if open fails.
 
105
 
 
106
   \param name map name
 
107
   \param mapset mapset name where raster map <i>name</i> lives
 
108
 
 
109
   \return nonnegative file descriptor (int)
 
110
 */
 
111
int Rast_open_old(const char *name, const char *mapset)
 
112
{
 
113
    int fd = Rast__open_old(name, mapset);
 
114
 
 
115
    /* turn on auto masking, if not already on */
 
116
    Rast__check_for_auto_masking();
 
117
    /*
 
118
       if(R__.auto_mask <= 0)
 
119
       R__.mask_buf = Rast_allocate_c_buf();
 
120
       now we don't ever free it!, so no need to allocate it  (Olga)
 
121
     */
 
122
    /* mask_buf is used for reading MASK file when mask is set and
 
123
       for reading map rows when the null file doesn't exist */
 
124
 
 
125
    return fd;
 
126
}
 
127
 
 
128
/*!  \brief Lower level function, open cell files, supercell files,
 
129
   and the MASK file.
 
130
 
 
131
   Actions:
 
132
   - opens the named cell file, following reclass reference if
 
133
   named layer is a reclass layer.
 
134
   - creates the required mapping between the data and the window
 
135
   for use by the get_map_row family of routines.
 
136
 
 
137
   Diagnostics: Errors other than actual open failure will cause a
 
138
   diagnostic to be delivered thru G_warning() open failure messages
 
139
   are left to the calling routine since the masking logic will want to
 
140
   issue a different warning.
 
141
 
 
142
   Note: This routine does NOT open the MASK layer. If it did we would
 
143
   get infinite recursion.  This routine is called to open the mask by
 
144
   Rast__check_for_auto_masking() which is called by Rast_open_old().
 
145
 
 
146
   \param name map name
 
147
   \param mapset mapset of cell file to be opened
 
148
 
 
149
   \return open file descriptor
 
150
 */
 
151
int Rast__open_old(const char *name, const char *mapset)
 
152
{
 
153
    struct fileinfo *fcb;
 
154
    int cell_fd, fd;
 
155
    char *cell_dir;
 
156
    const char *r_name;
 
157
    const char *r_mapset;
 
158
    struct Cell_head cellhd;
 
159
    int CELL_nbytes = 0;        /* bytes per cell in CELL map */
 
160
    int INTERN_SIZE;
 
161
    int reclass_flag;
 
162
    int MAP_NBYTES;
 
163
    RASTER_MAP_TYPE MAP_TYPE;
 
164
    struct Reclass reclass;
 
165
    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
166
    struct GDAL_link *gdal;
 
167
 
 
168
    Rast__init();
 
169
 
 
170
    G_unqualified_name(name, mapset, xname, xmapset);
 
171
    name = xname;
 
172
    mapset = xmapset;
 
173
 
 
174
    if (!G_find_raster2(name, mapset))
 
175
        G_fatal_error(_("Raster map <%s> not found"),
 
176
                      G_fully_qualified_name(name, mapset));
 
177
 
 
178
    /* Check for reclassification */
 
179
    reclass_flag = Rast_get_reclass(name, mapset, &reclass);
 
180
 
 
181
    switch (reclass_flag) {
 
182
    case 0:
 
183
        r_name = name;
 
184
        r_mapset = mapset;
 
185
        break;
 
186
    case 1:
 
187
        r_name = reclass.name;
 
188
        r_mapset = reclass.mapset;
 
189
        if (!G_find_raster2(r_name, r_mapset))
 
190
            G_fatal_error(_("Unable to open raster map <%s@%s> since it is a reclass "
 
191
                            "of raster map <%s@%s> which does not exist"),
 
192
                          name, mapset, r_name, r_mapset);
 
193
        break;
 
194
    default:                    /* Error reading cellhd/reclass file */
 
195
        G_fatal_error(_("Error reading reclass file for raster map <%s>"),
 
196
                      G_fully_qualified_name(name, mapset));
 
197
        break;
 
198
    }
 
199
 
 
200
    /* read the cell header */
 
201
    Rast_get_cellhd(r_name, r_mapset, &cellhd);
 
202
 
 
203
    /* now check the type */
 
204
    MAP_TYPE = Rast_map_type(r_name, r_mapset);
 
205
    if (MAP_TYPE < 0)
 
206
        G_fatal_error(_("Error reading map type for raster map <%s>"),
 
207
                      G_fully_qualified_name(name, mapset));
 
208
 
 
209
    if (MAP_TYPE == CELL_TYPE)
 
210
        /* set the number of bytes for CELL map */
 
211
    {
 
212
        CELL_nbytes = cellhd.format + 1;
 
213
        if (CELL_nbytes < 1)
 
214
            G_fatal_error(_("Raster map <%s@%s>: format field in header file invalid"),
 
215
                          r_name, r_mapset);
 
216
    }
 
217
 
 
218
    if (cellhd.proj != R__.rd_window.proj)
 
219
        G_fatal_error(_("Raster map <%s> is in different projection than current region. "
 
220
                        "Found <%s>, should be <%s>."),
 
221
                      G_fully_qualified_name(name, mapset),
 
222
                      G_projection_name(cellhd.proj),
 
223
                      G_projection_name(R__.rd_window.proj));
 
224
 
 
225
    if (cellhd.zone != R__.rd_window.zone)
 
226
        G_fatal_error(_("Raster map <%s> is in different zone (%d) than current region (%d)"),
 
227
                      G_fully_qualified_name(name, mapset), cellhd.zone, R__.rd_window.zone);
 
228
 
 
229
    /* when map is int warn if too large cell size */
 
230
    if (MAP_TYPE == CELL_TYPE && (unsigned int)CELL_nbytes > sizeof(CELL))
 
231
        G_fatal_error(_("Raster map <%s>: bytes per cell (%d) too large"),
 
232
                      G_fully_qualified_name(name, mapset), CELL_nbytes);
 
233
 
 
234
    /* record number of bytes per cell */
 
235
    if (MAP_TYPE == FCELL_TYPE) {
 
236
        cell_dir = "fcell";
 
237
        INTERN_SIZE = sizeof(FCELL);
 
238
        MAP_NBYTES = XDR_FLOAT_NBYTES;
 
239
    }
 
240
    else if (MAP_TYPE == DCELL_TYPE) {
 
241
        cell_dir = "fcell";
 
242
        INTERN_SIZE = sizeof(DCELL);
 
243
        MAP_NBYTES = XDR_DOUBLE_NBYTES;
 
244
    }
 
245
    else {                      /* integer */
 
246
        cell_dir = "cell";
 
247
        INTERN_SIZE = sizeof(CELL);
 
248
        MAP_NBYTES = CELL_nbytes;
 
249
    }
 
250
 
 
251
    gdal = Rast_get_gdal_link(r_name, r_mapset);
 
252
    if (gdal) {
 
253
#ifdef HAVE_GDAL
 
254
        cell_fd = -1;
 
255
#else
 
256
        G_fatal_error(_("Raster map <%s@%s> is a GDAL link but GRASS is compiled without GDAL support"),
 
257
                      r_name, r_mapset);
 
258
#endif
 
259
    }
 
260
    else {
 
261
        /* now actually open file for reading */
 
262
        cell_fd = G_open_old(cell_dir, r_name, r_mapset);
 
263
        if (cell_fd < 0)
 
264
            G_fatal_error(_("Unable to open %s file for raster map <%s@%s>"),
 
265
                          cell_dir, r_name, r_mapset);
 
266
    }
 
267
 
 
268
    fd = new_fileinfo();
 
269
    fcb = &R__.fileinfo[fd];
 
270
    fcb->data_fd = cell_fd;
 
271
 
 
272
    fcb->map_type = MAP_TYPE;
 
273
 
 
274
    /* Save cell header */
 
275
    fcb->cellhd = cellhd;
 
276
 
 
277
    /* allocate null bitstream buffers for reading null rows */
 
278
    fcb->null_fd = -1;
 
279
    fcb->null_cur_row = -1;
 
280
    fcb->null_bits = Rast__allocate_null_bits(cellhd.cols);
 
281
 
 
282
    /* mark closed */
 
283
    fcb->open_mode = -1;
 
284
 
 
285
    /* save name and mapset */
 
286
    fcb->name = G_store(name);
 
287
    fcb->mapset = G_store(mapset);
 
288
 
 
289
    /* mark no data row in memory  */
 
290
    fcb->cur_row = -1;
 
291
 
 
292
    /* if reclass, copy reclass structure */
 
293
    if ((fcb->reclass_flag = reclass_flag))
 
294
        fcb->reclass = reclass;
 
295
 
 
296
    fcb->gdal = gdal;
 
297
    if (!gdal)
 
298
        /* check for compressed data format, making initial reads if necessary */
 
299
        if (Rast__check_format(fd) < 0) {
 
300
            close(cell_fd);     /* warning issued by check_format() */
 
301
            G_fatal_error(_("Error reading format for <%s@%s>"),
 
302
                          r_name, r_mapset);
 
303
        }
 
304
 
 
305
    /* create the mapping from cell file to window */
 
306
    Rast__create_window_mapping(fd);
 
307
 
 
308
    /*
 
309
     * allocate the data buffer
 
310
     * number of bytes per cell is cellhd.format+1
 
311
     */
 
312
 
 
313
    /* for reading fcb->data is allocated to be fcb->cellhd.cols * fcb->nbytes 
 
314
       (= XDR_FLOAT/DOUBLE_NBYTES) */
 
315
    fcb->data = (unsigned char *)G_calloc(fcb->cellhd.cols, MAP_NBYTES);
 
316
 
 
317
    /* initialize/read in quant rules for float point maps */
 
318
    if (fcb->map_type != CELL_TYPE) {
 
319
        if (fcb->reclass_flag)
 
320
            Rast_read_quant(fcb->reclass.name, fcb->reclass.mapset,
 
321
                            &(fcb->quant));
 
322
        else
 
323
            Rast_read_quant(fcb->name, fcb->mapset, &(fcb->quant));
 
324
    }
 
325
 
 
326
    /* now mark open for read: this must follow create_window_mapping() */
 
327
    fcb->open_mode = OPEN_OLD;
 
328
    fcb->io_error = 0;
 
329
    fcb->map_type = MAP_TYPE;
 
330
    fcb->nbytes = MAP_NBYTES;
 
331
 
 
332
    if (!gdal) {
 
333
        if (!G_find_file2_misc("cell_misc", NULL_FILE, r_name, r_mapset)) {
 
334
            /* G_warning("unable to find [%s]",path); */
 
335
            fcb->null_file_exists = 0;
 
336
        }
 
337
        else {
 
338
            fcb->null_fd = G_open_old_misc("cell_misc", NULL_FILE, r_name, r_mapset);
 
339
            fcb->null_file_exists = fcb->null_fd >= 0;
 
340
        }
 
341
    }
 
342
 
 
343
    return fd;
 
344
}
 
345
 
 
346
/*!
 
347
   \brief Opens a new cell file in a database (compressed)
 
348
 
 
349
   Opens a new cell file <i>name</i> in the current mapset for writing
 
350
   by Rast_put_row().
 
351
 
 
352
   The file is created and filled with no data it is assumed that the
 
353
   new cell file is to conform to the current window.
 
354
 
 
355
   The file must be written sequentially. Use Rast_open_new_random()
 
356
   for non sequential writes.
 
357
 
 
358
   Note: the open actually creates a temporary file Rast_close() will
 
359
   move the temporary file to the cell file and write out the necessary
 
360
   support files (cellhd, cats, hist, etc.).
 
361
 
 
362
   Diagnostics: warning message printed if open fails
 
363
 
 
364
   Warning: calls to Rast_set_window() made after opening a new cell file
 
365
   may create confusion and should be avoided the new cell file will be
 
366
   created to conform to the window at the time of the open.
 
367
 
 
368
   \param name map name
 
369
 
 
370
   \return open file descriptor ( >= 0) if successful
 
371
   \return negative integer if error
 
372
 */
 
373
int Rast_open_c_new(const char *name)
 
374
{
 
375
    return open_raster_new(name, OPEN_NEW_COMPRESSED, CELL_TYPE);
 
376
}
 
377
 
 
378
/*!
 
379
   \brief Opens a new cell file in a database (uncompressed)
 
380
 
 
381
   See also Rast_open_new().
 
382
 
 
383
   \param name map name
 
384
 
 
385
   \return open file descriptor ( >= 0) if successful
 
386
   \return negative integer if error
 
387
 */
 
388
int Rast_open_c_new_uncompressed(const char *name)
 
389
{
 
390
    return open_raster_new(name, OPEN_NEW_UNCOMPRESSED, CELL_TYPE);
 
391
}
 
392
 
 
393
/*!
 
394
   \brief Save histogram for newly create raster map (cell)
 
395
 
 
396
   If newly created cell files should have histograms, set flag=1
 
397
   otherwise set flag=0. Applies to subsequent opens.
 
398
 
 
399
   \param flag flag indicator
 
400
 */
 
401
void Rast_want_histogram(int flag)
 
402
{
 
403
    R__.want_histogram = flag;
 
404
}
 
405
 
 
406
/*!
 
407
   \brief Sets the format for subsequent opens on new integer cell files
 
408
   (uncompressed and random only).
 
409
 
 
410
   Warning: subsequent put_row calls will only write n+1 bytes per
 
411
   cell. If the data requires more, the cell file will be written
 
412
   incorrectly (but with n+1 bytes per cell)
 
413
 
 
414
   When writing float map: format is -1
 
415
 
 
416
   \param n format
 
417
 */
 
418
void Rast_set_cell_format(int n)
 
419
/* sets the format for integer raster map */
 
420
{
 
421
    R__.nbytes = n + 1;
 
422
    if (R__.nbytes <= 0)
 
423
        R__.nbytes = 1;
 
424
    if (R__.nbytes > sizeof(CELL))
 
425
        R__.nbytes = sizeof(CELL);
 
426
}
 
427
 
 
428
/*!
 
429
   \brief Get cell value format
 
430
 
 
431
   \param v cell
 
432
 
 
433
   \return cell format
 
434
 */
 
435
int Rast_get_cell_format(CELL v)
 
436
{
 
437
    unsigned int i;
 
438
 
 
439
    if (v >= 0)
 
440
        for (i = 0; i < sizeof(CELL); i++)
 
441
            if (!(v /= 256))
 
442
                return i;
 
443
    return sizeof(CELL) - 1;
 
444
}
 
445
 
 
446
/*!
 
447
   \brief Opens new fcell file in a database
 
448
 
 
449
   Opens a new floating-point map <i>name</i> in the current mapset for
 
450
   writing. The type of the file (i.e. either double or float) is
 
451
   determined and fixed at this point. The default is FCELL_TYPE. In
 
452
   order to change this default
 
453
 
 
454
   Use Rast_set_fp_type() where type is one of DCELL_TYPE or FCELL_TYPE.
 
455
 
 
456
   See warnings and notes for Rast_open_new().
 
457
 
 
458
   \param name map name
 
459
 
 
460
   \return nonnegative file descriptor (int)
 
461
 */
 
462
int Rast_open_fp_new(const char *name)
 
463
{
 
464
    return open_raster_new(name, OPEN_NEW_COMPRESSED, R__.fp_type);
 
465
}
 
466
 
 
467
/*!
 
468
   \brief Opens new fcell file in a database (uncompressed)
 
469
 
 
470
   See Rast_open_fp_new() for details.
 
471
 
 
472
   \param name map name
 
473
 
 
474
   \return nonnegative file descriptor (int)
 
475
 */
 
476
int Rast_open_fp_new_uncompressed(const char *name)
 
477
{
 
478
    return open_raster_new(name, OPEN_NEW_UNCOMPRESSED, R__.fp_type);
 
479
}
 
480
 
 
481
#ifdef HAVE_GDAL
 
482
static int open_raster_new_gdal(char *map, char *mapset,
 
483
                                RASTER_MAP_TYPE map_type)
 
484
{
 
485
    int fd;
 
486
    struct fileinfo *fcb;
 
487
 
 
488
    fd = new_fileinfo();
 
489
    fcb = &R__.fileinfo[fd];
 
490
    fcb->data_fd = -1;
 
491
 
 
492
    /* mark closed */
 
493
    fcb->map_type = map_type;
 
494
    fcb->open_mode = -1;
 
495
 
 
496
    fcb->gdal = Rast_create_gdal_link(map, map_type);
 
497
    if (!fcb->gdal)
 
498
        G_fatal_error(_("Unable to create GDAL link"));
 
499
 
 
500
    fcb->cellhd = R__.wr_window;
 
501
    fcb->cellhd.compressed = 0;
 
502
    fcb->nbytes = Rast_cell_size(fcb->map_type);
 
503
    /* for writing fcb->data is allocated to be R__.wr_window.cols * 
 
504
       sizeof(CELL or DCELL or FCELL)  */
 
505
    fcb->data = G_calloc(R__.wr_window.cols, fcb->nbytes);
 
506
 
 
507
    fcb->name = map;
 
508
    fcb->mapset = mapset;
 
509
    fcb->cur_row = 0;
 
510
 
 
511
    fcb->row_ptr = NULL;
 
512
    fcb->temp_name = NULL;
 
513
    fcb->null_temp_name = NULL;
 
514
    fcb->null_cur_row = 0;
 
515
    fcb->null_bits = NULL;
 
516
    fcb->null_fd = -1;
 
517
 
 
518
    if (fcb->map_type != CELL_TYPE)
 
519
        Rast_quant_init(&(fcb->quant));
 
520
 
 
521
    /* init cell stats */
 
522
    /* now works only for int maps */
 
523
    if (fcb->map_type == CELL_TYPE)
 
524
        if ((fcb->want_histogram = R__.want_histogram))
 
525
            Rast_init_cell_stats(&fcb->statf);
 
526
 
 
527
    /* init range and if map is double/float init d/f_range */
 
528
    Rast_init_range(&fcb->range);
 
529
 
 
530
    if (fcb->map_type != CELL_TYPE)
 
531
        Rast_init_fp_range(&fcb->fp_range);
 
532
 
 
533
    /* mark file as open for write */
 
534
    fcb->open_mode = OPEN_NEW_UNCOMPRESSED;
 
535
    fcb->io_error = 0;
 
536
 
 
537
    return fd;
 
538
}
 
539
#endif /* HAVE_GDAL */
 
540
 
 
541
static int open_raster_new(const char *name, int open_mode,
 
542
                           RASTER_MAP_TYPE map_type)
 
543
{
 
544
    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
545
    struct fileinfo *fcb;
 
546
    int fd, cell_fd;
 
547
    char *tempname;
 
548
    char *map;
 
549
    char *mapset;
 
550
    const char *cell_dir;
 
551
    int nbytes;
 
552
 
 
553
    Rast__init();
 
554
 
 
555
    switch (map_type) {
 
556
    case CELL_TYPE:
 
557
        cell_dir = "cell";
 
558
        nbytes = R__.nbytes;
 
559
        break;
 
560
    case FCELL_TYPE:
 
561
        nbytes = XDR_FLOAT_NBYTES;
 
562
        cell_dir = "fcell";
 
563
        break;
 
564
    case DCELL_TYPE:
 
565
        nbytes = XDR_DOUBLE_NBYTES;
 
566
        cell_dir = "fcell";
 
567
        break;
 
568
    default:
 
569
        G_fatal_error(_("Invalid map type <%d>"), map_type);
 
570
        break;
 
571
    }
 
572
 
 
573
    if (G_unqualified_name(name, G_mapset(), xname, xmapset) < 0)
 
574
        G_fatal_error(_("Raster map <%s> is not in the current mapset (%s)"),
 
575
                      name, G_mapset());
 
576
    map = G_store(xname);
 
577
    mapset = G_store(xmapset);
 
578
 
 
579
    /* check for legal grass name */
 
580
    if (G_legal_filename(map) < 0)
 
581
        G_fatal_error(_("<%s> is an illegal file name"), map);
 
582
 
 
583
#ifdef HAVE_GDAL
 
584
    if (G_find_file2("", "GDAL", G_mapset()))
 
585
        return open_raster_new_gdal(map, mapset, map_type);
 
586
#endif
 
587
 
 
588
    /* open a tempfile name */
 
589
    tempname = G_tempfile();
 
590
    cell_fd = creat(tempname, 0666);
 
591
    if (cell_fd < 0) {
 
592
        G_free(mapset);
 
593
        G_free(tempname);
 
594
        G_free(map);
 
595
        G_fatal_error(_("No temp files available: %s"), strerror(errno));
 
596
    }
 
597
 
 
598
    fd = new_fileinfo();
 
599
    fcb = &R__.fileinfo[fd];
 
600
    fcb->data_fd = cell_fd;
 
601
 
 
602
    /*
 
603
     * since we are bypassing the normal open logic
 
604
     * must create the cell element 
 
605
     */
 
606
    G_make_mapset_element(cell_dir);
 
607
 
 
608
    /* mark closed */
 
609
    fcb->map_type = map_type;
 
610
    fcb->open_mode = -1;
 
611
    fcb->gdal = NULL;
 
612
 
 
613
    /* for writing fcb->data is allocated to be R__.wr_window.cols * 
 
614
       sizeof(CELL or DCELL or FCELL)  */
 
615
    fcb->data = (unsigned char *)G_calloc(R__.wr_window.cols,
 
616
                                          Rast_cell_size(fcb->map_type));
 
617
 
 
618
    /*
 
619
     * copy current window into cell header
 
620
     * set format to cell/supercell
 
621
     * for compressed writing
 
622
     *   allocate space to hold the row address array
 
623
     */
 
624
    fcb->cellhd = R__.wr_window;
 
625
 
 
626
    if (open_mode == OPEN_NEW_COMPRESSED && fcb->map_type == CELL_TYPE) {
 
627
        fcb->row_ptr = G_calloc(fcb->cellhd.rows + 1, sizeof(off_t));
 
628
        G_zero(fcb->row_ptr, (fcb->cellhd.rows + 1) * sizeof(off_t));
 
629
        Rast__write_row_ptrs(fd);
 
630
        fcb->cellhd.compressed = R__.compression_type;
 
631
 
 
632
        fcb->nbytes = 1;        /* to the minimum */
 
633
    }
 
634
    else {
 
635
        fcb->nbytes = nbytes;
 
636
        if (open_mode == OPEN_NEW_COMPRESSED) {
 
637
            fcb->row_ptr = G_calloc(fcb->cellhd.rows + 1, sizeof(off_t));
 
638
            G_zero(fcb->row_ptr, (fcb->cellhd.rows + 1) * sizeof(off_t));
 
639
            Rast__write_row_ptrs(fd);
 
640
            fcb->cellhd.compressed = R__.compression_type;
 
641
        }
 
642
        else
 
643
            fcb->cellhd.compressed = 0;
 
644
 
 
645
        if (fcb->map_type != CELL_TYPE) {
 
646
            Rast_quant_init(&(fcb->quant));
 
647
        }
 
648
    }
 
649
 
 
650
    /* save name and mapset, and tempfile name */
 
651
    fcb->name = map;
 
652
    fcb->mapset = mapset;
 
653
    fcb->temp_name = tempname;
 
654
 
 
655
    /* next row to be written (in order) is zero */
 
656
    fcb->cur_row = 0;
 
657
 
 
658
    /* open a null tempfile name */
 
659
    tempname = G_tempfile();
 
660
    fcb->null_fd = creat(tempname, 0666);
 
661
    if (fcb->null_fd < 0) {
 
662
        G_free(tempname);
 
663
        G_free(fcb->name);
 
664
        G_free(fcb->mapset);
 
665
        G_free(fcb->temp_name);
 
666
        close(cell_fd);
 
667
        G_fatal_error(_("No temp files available: %s"), strerror(errno));
 
668
    }
 
669
 
 
670
    fcb->null_temp_name = tempname;
 
671
 
 
672
    /* next row to be written (in order) is zero */
 
673
    fcb->null_cur_row = 0;
 
674
 
 
675
    /* allocate null bitstream buffer for writing */
 
676
    fcb->null_bits = Rast__allocate_null_bits(fcb->cellhd.cols);
 
677
 
 
678
    /* init cell stats */
 
679
    /* now works only for int maps */
 
680
    if (fcb->map_type == CELL_TYPE)
 
681
        if ((fcb->want_histogram = R__.want_histogram))
 
682
            Rast_init_cell_stats(&fcb->statf);
 
683
 
 
684
    /* init range and if map is double/float init d/f_range */
 
685
    Rast_init_range(&fcb->range);
 
686
 
 
687
    if (fcb->map_type != CELL_TYPE)
 
688
        Rast_init_fp_range(&fcb->fp_range);
 
689
 
 
690
    /* mark file as open for write */
 
691
    fcb->open_mode = open_mode;
 
692
    fcb->io_error = 0;
 
693
 
 
694
    return fd;
 
695
}
 
696
 
 
697
/*!
 
698
   \brief Set raster map floating-point data format.
 
699
 
 
700
   This controls the storage type for floating-point maps. It affects
 
701
   subsequent calls to G_open_fp_map_new(). The <i>type</i> must be
 
702
   one of FCELL_TYPE (float) or DCELL_TYPE (double). The use of this
 
703
   routine by applications is discouraged since its use would override
 
704
   user preferences.
 
705
 
 
706
   \param type raster data type
 
707
 
 
708
   \return void
 
709
 */
 
710
void Rast_set_fp_type(RASTER_MAP_TYPE map_type)
 
711
{
 
712
    Rast__init();
 
713
 
 
714
    switch (map_type) {
 
715
    case FCELL_TYPE:
 
716
    case DCELL_TYPE:
 
717
        R__.fp_type = map_type;
 
718
        break;
 
719
    default:
 
720
        G_fatal_error(_("Rast_set_fp_type(): can only be called with FCELL_TYPE or DCELL_TYPE"));
 
721
        break;
 
722
    }
 
723
}
 
724
 
 
725
/*!
 
726
   \brief Check if raster map is floating-point
 
727
 
 
728
   Returns true (1) if raster map <i>name</i> in <i>mapset</i>
 
729
   is a floating-point dataset; false(0) otherwise.
 
730
 
 
731
   \param name map name
 
732
   \param mapset mapset name
 
733
 
 
734
   \return 1 floating-point
 
735
   \return 0 int
 
736
 */
 
737
int Rast_map_is_fp(const char *name, const char *mapset)
 
738
{
 
739
    char path[GPATH_MAX];
 
740
    const char *xmapset;
 
741
 
 
742
    xmapset = G_find_raster2(name, mapset);
 
743
    if (!xmapset)
 
744
        G_fatal_error(_("Raster map <%s> not found"),
 
745
                      G_fully_qualified_name(name, mapset));
 
746
 
 
747
    G_file_name(path, "fcell", name, xmapset);
 
748
    if (access(path, 0) == 0)
 
749
        return 1;
 
750
 
 
751
    G_file_name(path, "g3dcell", name, xmapset);
 
752
    if (access(path, 0) == 0)
 
753
        return 1;
 
754
 
 
755
    return 0;
 
756
}
 
757
 
 
758
/*!
 
759
   \brief Determine raster data type
 
760
 
 
761
   Determines if the raster map is floating point or integer. Returns
 
762
   DCELL_TYPE for double maps, FCELL_TYPE for float maps, CELL_TYPE for
 
763
   integer maps, -1 if error has occurred
 
764
 
 
765
   \param name map name 
 
766
   \param mapset mapset where map <i>name</i> lives
 
767
 
 
768
   \return raster data type
 
769
 */
 
770
RASTER_MAP_TYPE Rast_map_type(const char *name, const char *mapset)
 
771
{
 
772
    char path[GPATH_MAX];
 
773
    const char *xmapset;
 
774
 
 
775
    xmapset = G_find_raster2(name, mapset);
 
776
    if (!xmapset) {
 
777
        if (mapset && *mapset)
 
778
            G_fatal_error(_("Raster map <%s> not found in mapset <%s>"),
 
779
                          name, mapset);
 
780
        else
 
781
            G_fatal_error(_("Raster map <%s> not found"), name);
 
782
    }
 
783
 
 
784
    G_file_name(path, "fcell", name, xmapset);
 
785
 
 
786
    if (access(path, 0) == 0)
 
787
        return Rast__check_fp_type(name, xmapset);
 
788
 
 
789
    G_file_name(path, "g3dcell", name, xmapset);
 
790
 
 
791
    if (access(path, 0) == 0)
 
792
        return DCELL_TYPE;
 
793
 
 
794
    return CELL_TYPE;
 
795
}
 
796
 
 
797
/*!
 
798
   \brief Determine raster type from descriptor
 
799
 
 
800
   Determines if the raster map is floating point or integer. Returns
 
801
   DCELL_TYPE for double maps, FCELL_TYPE for float maps, CELL_TYPE for
 
802
   integer maps, -1 if error has occurred
 
803
 
 
804
   \param fd file descriptor
 
805
 
 
806
   \return raster data type
 
807
 */
 
808
RASTER_MAP_TYPE Rast_get_map_type(int fd)
 
809
{
 
810
    struct fileinfo *fcb = &R__.fileinfo[fd];
 
811
 
 
812
    return fcb->map_type;
 
813
}
 
814
 
 
815
/*!
 
816
   \brief Determines whether the floating points cell file has double or float type
 
817
 
 
818
   \param name map name
 
819
   \param mapset mapset where map <i>name</i> lives
 
820
 
 
821
   \return raster type (fcell, dcell)
 
822
 */
 
823
RASTER_MAP_TYPE Rast__check_fp_type(const char *name, const char *mapset)
 
824
{
 
825
    char path[GPATH_MAX];
 
826
    struct Key_Value *format_keys;
 
827
    const char *str, *str1;
 
828
    RASTER_MAP_TYPE map_type;
 
829
    const char *xmapset;
 
830
 
 
831
    xmapset = G_find_raster2(name, mapset);
 
832
    if (!xmapset)
 
833
        G_fatal_error(_("Raster map <%s> not found"),
 
834
                      G_fully_qualified_name(name, mapset));
 
835
 
 
836
    G_file_name_misc(path, "cell_misc", FORMAT_FILE, name, xmapset);
 
837
 
 
838
    if (access(path, 0) != 0)
 
839
        G_fatal_error(_("Unable to find '%s'"), path);
 
840
 
 
841
    format_keys = G_read_key_value_file(path);
 
842
 
 
843
    if ((str = G_find_key_value("type", format_keys)) != NULL) {
 
844
        if (strcmp(str, "double") == 0)
 
845
            map_type = DCELL_TYPE;
 
846
        else if (strcmp(str, "float") == 0)
 
847
            map_type = FCELL_TYPE;
 
848
        else {
 
849
            G_free_key_value(format_keys);
 
850
            G_fatal_error(_("Invalid type: field '%s' in file '%s'"), str, path);
 
851
        }
 
852
    }
 
853
    else {
 
854
        G_free_key_value(format_keys);
 
855
        G_fatal_error(_("Missing type: field in file '%s'"), path);
 
856
    }
 
857
 
 
858
    if ((str1 = G_find_key_value("byte_order", format_keys)) != NULL) {
 
859
        if (strcmp(str1, "xdr") != 0)
 
860
            G_warning(_("Raster map <%s> is not xdr: byte_order: %s"),
 
861
                      name, str);
 
862
        /* here read and translate  byte order if not using xdr */
 
863
    }
 
864
    G_free_key_value(format_keys);
 
865
    return map_type;
 
866
}
 
867
 
 
868
/*!
 
869
   \brief Opens a new raster map
 
870
 
 
871
   Opens a new raster map of type <i>wr_type</i>
 
872
 
 
873
   See warnings and notes for Rast_open_new().
 
874
 
 
875
   Supported data types:
 
876
   - CELL_TYPE
 
877
   - FCELL_TYPE
 
878
   - DCELL_TYPE
 
879
 
 
880
   On CELL_TYPE calls Rast_open_new() otherwise Rast_open_fp_new().
 
881
 
 
882
   \param name map name
 
883
   \param wr_type raster data type
 
884
 
 
885
   \return nonnegative file descriptor (int)
 
886
 */
 
887
int Rast_open_new(const char *name, RASTER_MAP_TYPE wr_type)
 
888
{
 
889
    return open_raster_new(name, OPEN_NEW_COMPRESSED, wr_type);
 
890
}
 
891
 
 
892
/*!
 
893
   \brief Opens a new raster map (uncompressed)
 
894
 
 
895
   See Rast_open_new().
 
896
 
 
897
   \param name map name
 
898
   \param wr_type raster data type
 
899
 
 
900
   \return nonnegative file descriptor (int)
 
901
 */
 
902
int Rast_open_new_uncompressed(const char *name, RASTER_MAP_TYPE wr_type)
 
903
{
 
904
    return open_raster_new(name, OPEN_NEW_UNCOMPRESSED, wr_type);
 
905
}
 
906
 
 
907
/*!
 
908
   \brief Sets quant translation rules for raster map opened for
 
909
   reading.
 
910
 
 
911
   Returned by Rast_open_old(). After calling this function,
 
912
   Rast_get_c_row() and Rast_get_c_row() will use rules defined by q
 
913
   (instead of using rules defined in map's quant file) to convert floats to
 
914
   ints.
 
915
 
 
916
   \param fd file descriptor (cell file)
 
917
   \param q pointer to Quant structure
 
918
 
 
919
   \return void
 
920
 */
 
921
void Rast_set_quant_rules(int fd, struct Quant *q)
 
922
{
 
923
    struct fileinfo *fcb = &R__.fileinfo[fd];
 
924
    CELL cell;
 
925
    DCELL dcell;
 
926
    struct Quant_table *p;
 
927
 
 
928
    if (fcb->open_mode != OPEN_OLD)
 
929
        G_fatal_error(_("Rast_set_quant_rules() can be called only for "
 
930
                        "raster maps opened for reading"));
 
931
 
 
932
    /* copy all info from q to fcb->quant) */
 
933
    Rast_quant_init(&fcb->quant);
 
934
    if (q->truncate_only) {
 
935
        Rast_quant_truncate(&fcb->quant);
 
936
        return;
 
937
    }
 
938
 
 
939
    for (p = &(q->table[q->nofRules - 1]); p >= q->table; p--)
 
940
        Rast_quant_add_rule(&fcb->quant, p->dLow, p->dHigh, p->cLow,
 
941
                            p->cHigh);
 
942
    if (Rast_quant_get_neg_infinite_rule(q, &dcell, &cell) > 0)
 
943
        Rast_quant_set_neg_infinite_rule(&fcb->quant, dcell, cell);
 
944
    if (Rast_quant_get_pos_infinite_rule(q, &dcell, &cell) > 0)
 
945
        Rast_quant_set_pos_infinite_rule(&fcb->quant, dcell, cell);
 
946
}