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