source: grass/trunk/raster/r.in.gdal/main.c

Last change on this file was 74373, checked in by mmetz, 5 years ago

r.in.gdal: +info on subdatasets, see #798
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 60.2 KB
Line 
1
2/****************************************************************************
3 *
4 * MODULE: r.in.gdal
5 *
6 * AUTHOR(S): Frank Warmerdam (copyright of this file)
7 * Added optional GCP transformation: Markus Neteler 10/2001
8 *
9 * PURPOSE: Imports many GIS/image formats into GRASS utilizing the GDAL
10 * library.
11 *
12 * COPYRIGHT: (C) 2001-2015 by Frank Warmerdam, and the GRASS Development Team
13 *
14 * This program is free software under the GNU General Public
15 * License (>=v2). Read the file COPYING that comes with GRASS
16 * for details.
17 *
18 *****************************************************************************/
19
20#include <stdlib.h>
21#include <stdio.h>
22#include <errno.h>
23#include <unistd.h>
24#include <math.h>
25#include <string.h>
26#include <grass/gis.h>
27#include <grass/raster.h>
28#include <grass/imagery.h>
29#include <grass/gprojects.h>
30#include <grass/glocale.h>
31
32#include <gdal.h>
33#include <cpl_conv.h>
34
35#undef MIN
36#undef MAX
37#define MIN(a,b) ((a) < (b) ? (a) : (b))
38#define MAX(a,b) ((a) > (b) ? (a) : (b))
39
40void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
41 char *outloc, int create_only, int override,
42 int check_only);
43static void ImportBand(GDALRasterBandH hBand, const char *output,
44 struct Ref *group_ref, int *rowmap, int *colmap,
45 int col_offset);
46static void SetupReprojector(const char *pszSrcWKT, const char *pszDstLoc,
47 struct pj_info *iproj, struct pj_info *oproj,
48 struct pj_info *tproj);
49static int dump_rat(GDALRasterBandH hBand, char *outrat, int nBand);
50static void error_handler_ds(void *p);
51static int l1bdriver;
52
53static GDALDatasetH opends(char *dsname, const char **doo, GDALDriverH *hDriver)
54{
55 GDALDatasetH hDS = NULL;
56
57#if GDAL_VERSION_NUM >= 2000000
58 hDS = GDALOpenEx(dsname, GDAL_OF_RASTER | GDAL_OF_READONLY, NULL,
59 doo, NULL);
60#else
61 hDS = GDALOpen(dsname, GA_ReadOnly);
62#endif
63 if (hDS == NULL)
64 G_fatal_error(_("Unable to open datasource <%s>"), dsname);
65 G_add_error_handler(error_handler_ds, hDS);
66
67 *hDriver = GDALGetDatasetDriver(hDS); /* needed for AVHRR data */
68 /* L1B - NOAA/AVHRR data must be treated differently */
69 /* for hDriver names see gdal/frmts/gdalallregister.cpp */
70 G_debug(3, "GDAL Driver: %s", GDALGetDriverShortName(*hDriver));
71 if (strcmp(GDALGetDriverShortName(*hDriver), "L1B") != 0)
72 l1bdriver = 0;
73 else {
74 l1bdriver = 1; /* AVHRR found, needs north south flip */
75 G_warning(_("Input seems to be NOAA/AVHRR data which needs to be "
76 "georeferenced with thin plate spline transformation "
77 "(%s or %s)."), "i.rectify -t", "gdalwarp -tps");
78 }
79
80 return hDS;
81}
82
83/************************************************************************/
84/* main() */
85
86/************************************************************************/
87
88int main(int argc, char *argv[])
89{
90 char *input;
91 char *output;
92 char *title;
93 struct Cell_head cellhd, cur_wind;
94 struct Key_Value *proj_info = NULL, *proj_units = NULL;
95 GDALDatasetH hDS;
96 GDALDriverH hDriver;
97 GDALRasterBandH hBand;
98 double adfGeoTransform[6];
99 int n_bands;
100 int force_imagery = FALSE;
101 int overwrite;
102 int offset = 0;
103 char *suffix;
104 int num_digits = 0;
105 int croptoregion, *rowmapall, *colmapall, *rowmap, *colmap, col_offset;
106 int roff, coff;
107 char **doo;
108
109 struct GModule *module;
110 struct
111 {
112 struct Option *input, *output, *target, *title, *outloc, *band,
113 *memory, *offset, *num_digits, *map_names_file,
114 *rat, *cfg, *doo;
115 } parm;
116 struct Flag *flag_o, *flag_e, *flag_k, *flag_f, *flag_l, *flag_c, *flag_p,
117 *flag_j, *flag_a, *flag_r;
118
119 /* -------------------------------------------------------------------- */
120 /* Initialize. */
121 /* -------------------------------------------------------------------- */
122 G_gisinit(argv[0]);
123
124 module = G_define_module();
125 G_add_keyword(_("raster"));
126 G_add_keyword(_("import"));
127 G_add_keyword(_("create location"));
128 module->description =
129 _("Imports raster data into a GRASS raster map using GDAL library.");
130
131 /* -------------------------------------------------------------------- */
132 /* Setup and fetch parameters. */
133 /* -------------------------------------------------------------------- */
134 parm.input = G_define_standard_option(G_OPT_F_BIN_INPUT);
135 parm.input->description = _("Name of raster file to be imported");
136
137 parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
138
139 parm.band = G_define_option();
140 parm.band->key = "band";
141 parm.band->type = TYPE_INTEGER;
142 parm.band->multiple = YES;
143 parm.band->required = NO;
144 parm.band->description = _("Band(s) to select (default is all bands)");
145 parm.band->guisection = _("Bands");
146
147 parm.memory = G_define_option();
148 parm.memory->key = "memory";
149 parm.memory->type = TYPE_INTEGER;
150 parm.memory->required = NO;
151#if GDAL_VERSION_NUM < 1800
152 parm.memory->options = "0-2047";
153#endif
154 parm.memory->answer = "300";
155 parm.memory->label = _("Maximum memory to be used (in MB)");
156 parm.memory->description = _("Cache size for raster rows");
157
158 parm.target = G_define_option();
159 parm.target->key = "target";
160 parm.target->type = TYPE_STRING;
161 parm.target->required = NO;
162 parm.target->label = _("Name of GCPs target location");
163 parm.target->description =
164 _("Name of location to create or to read projection from for GCPs transformation");
165 parm.target->key_desc = "name";
166 parm.target->guisection = _("Projection");
167
168 parm.title = G_define_option();
169 parm.title->key = "title";
170 parm.title->key_desc = "phrase";
171 parm.title->type = TYPE_STRING;
172 parm.title->required = NO;
173 parm.title->description = _("Title for resultant raster map");
174 parm.title->guisection = _("Metadata");
175
176 parm.offset = G_define_option();
177 parm.offset->key = "offset";
178 parm.offset->type = TYPE_INTEGER;
179 parm.offset->required = NO;
180 parm.offset->answer = "0";
181 parm.offset->label = _("Offset to be added to band numbers");
182 parm.offset->description = _("If 0, no offset is added and the first band is 1");
183 parm.offset->guisection = _("Metadata");
184
185 parm.num_digits = G_define_option();
186 parm.num_digits->key = "num_digits";
187 parm.num_digits->type = TYPE_INTEGER;
188 parm.num_digits->required = NO;
189 parm.num_digits->answer = "0";
190 parm.num_digits->label = _("Zero-padding of band number by filling with leading zeros up to given number");
191 parm.num_digits->description = _("If 0, length will be adjusted to 'offset' number without leading zeros");
192 parm.num_digits->guisection = _("Metadata");
193
194 parm.map_names_file = G_define_standard_option(G_OPT_F_OUTPUT);
195 parm.map_names_file->key = "map_names_file";
196 parm.map_names_file->required = NO;
197 parm.map_names_file->description = _("Name of the output file that contains the imported map names");
198 parm.map_names_file->guisection = _("Metadata");
199
200 parm.outloc = G_define_option();
201 parm.outloc->key = "location";
202 parm.outloc->type = TYPE_STRING;
203 parm.outloc->required = NO;
204 parm.outloc->description = _("Name for new location to create");
205 parm.outloc->key_desc = "name";
206
207 parm.rat = G_define_option();
208 parm.rat->key = "table";
209 parm.rat->type = TYPE_STRING;
210 parm.rat->required = NO;
211 parm.rat->label = _("File prefix for raster attribute tables");
212 parm.rat->description = _("The band number and \".csv\" will be appended to the file prefix");
213 parm.rat->key_desc = "file";
214
215 parm.cfg = G_define_option();
216 parm.cfg->key = "gdal_config";
217 parm.cfg->type = TYPE_STRING;
218 parm.cfg->required = NO;
219 parm.cfg->label = _("GDAL configuration options");
220 parm.cfg->description = _("Comma-separated list of key=value pairs");
221
222 parm.doo = G_define_option();
223 parm.doo->key = "gdal_doo";
224 parm.doo->type = TYPE_STRING;
225 parm.doo->required = NO;
226 parm.doo->label = _("GDAL dataset open options");
227 parm.doo->description = _("Comma-separated list of key=value pairs");
228
229 flag_o = G_define_flag();
230 flag_o->key = 'o';
231 flag_o->label =
232 _("Override projection check (use current location's projection)");
233 flag_o->description =
234 _("Assume that the dataset has same projection as the current location");
235 flag_o->guisection = _("Projection");
236
237 flag_j = G_define_flag();
238 flag_j->key = 'j';
239 flag_j->description =
240 _("Perform projection check only and exit");
241 flag_j->suppress_required = YES;
242 flag_j->guisection = _("Projection");
243 G_option_requires(flag_j, parm.input, NULL);
244
245 flag_e = G_define_flag();
246 flag_e->key = 'e';
247 flag_e->label = _("Extend region extents based on new dataset");
248 flag_e->description =
249 _("Also updates the default region if in the PERMANENT mapset");
250 flag_e->guisection = _("Region");
251
252 flag_f = G_define_flag();
253 flag_f->key = 'f';
254 flag_f->description = _("List supported formats and exit");
255 flag_f->guisection = _("Print");
256 flag_f->suppress_required = YES;
257
258 flag_l = G_define_flag();
259 flag_l->key = 'l';
260 flag_l->description =
261 _("Force Lat/Lon maps to fit into geographic coordinates (90N,S; 180E,W)");
262
263 flag_a = G_define_flag();
264 flag_a->key = 'a';
265 flag_a->label = _("Auto-adjustment for lat/lon");
266 flag_a->description = _("Attempt to fix small precision errors in resolution and extents");
267
268 flag_k = G_define_flag();
269 flag_k->key = 'k';
270 flag_k->description =
271 _("Keep band numbers instead of using band color names");
272 flag_k->guisection = _("Bands");
273
274 flag_c = G_define_flag();
275 flag_c->key = 'c';
276 flag_c->description =
277 _("Create the location specified by the \"location\" parameter and exit."
278 " Do not import the raster file.");
279
280 flag_r = G_define_flag();
281 flag_r->key = 'r';
282 flag_r->description = _("Limit import to the current region");
283 flag_r->guisection = _("Region");
284
285 flag_p = G_define_flag();
286 flag_p->key = 'p';
287 flag_p->description = _("Print number of bands and exit");
288 flag_p->guisection = _("Print");
289 flag_p->suppress_required = YES;
290 G_option_requires(flag_p, parm.input, NULL);
291
292 /* 1. list supported formats */
293 /* 2. print input bands */
294 /* 3. open input */
295 /* 4. check projection / create location */
296 /* 5. open output */
297
298 /* The parser checks if the map already exists in current mapset, this is
299 * wrong if location options is used, so we switch out the check and do it
300 * in the module after the parser */
301 overwrite = G_check_overwrite(argc, argv);
302
303 if (G_parser(argc, argv))
304 exit(EXIT_FAILURE);
305
306 /* -------------------------------------------------------------------- */
307 /* Fire up the engines. */
308 /* -------------------------------------------------------------------- */
309 GDALAllRegister();
310
311 /* -------------------------------------------------------------------- */
312 /* List supported formats and exit. */
313 /* code from GDAL 1.2.5 gcore/gdal_misc.cpp */
314 /* Copyright (c) 1999, Frank Warmerdam */
315 /* -------------------------------------------------------------------- */
316 if (flag_f->answer) {
317 int iDr;
318
319 G_message(_("Supported formats:"));
320 for (iDr = 0; iDr < GDALGetDriverCount(); iDr++) {
321 const char *pszRWFlag;
322
323 hDriver = GDALGetDriver(iDr);
324
325#ifdef GDAL_DCAP_RASTER
326 /* Starting with GDAL 2.0, vector drivers can also be returned */
327 /* Only keep raster drivers */
328 if (!GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, NULL))
329 continue;
330#endif
331
332 if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, NULL))
333 pszRWFlag = "rw+";
334 else if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, NULL))
335 pszRWFlag = "rw";
336 else
337 pszRWFlag = "ro";
338
339 fprintf(stdout, " %s (%s): %s\n",
340 GDALGetDriverShortName(hDriver),
341 pszRWFlag, GDALGetDriverLongName(hDriver));
342 }
343 exit(EXIT_SUCCESS);
344 }
345
346 input = parm.input->answer;
347
348 output = parm.output->answer;
349
350 offset = atoi(parm.offset->answer);
351
352 num_digits = atoi(parm.num_digits->answer);
353
354 if ((title = parm.title->answer))
355 G_strip(title);
356
357 /* -------------------------------------------------------------------- */
358 /* Do some additional parameter validation. */
359 /* -------------------------------------------------------------------- */
360 if (parm.target->answer && parm.outloc->answer
361 && strcmp(parm.target->answer, parm.outloc->answer) == 0) {
362 G_fatal_error(_("You have to specify a target location different from output location"));
363 }
364
365 if (flag_c->answer && parm.outloc->answer == NULL ) {
366 G_fatal_error(_("You need to specify valid location name."));
367 }
368
369 if (flag_l->answer && G_projection() != PROJECTION_LL)
370 G_fatal_error(_("The '-l' flag only works in Lat/Lon locations"));
371
372 if (num_digits < 0)
373 G_fatal_error(_("The number of digits for band numbering must be equal or greater than 0"));
374
375 /* Allocate the suffix string */
376 if (num_digits > 0) {
377 suffix = G_calloc( num_digits + 1, sizeof(char));
378 }
379 else {
380 /* Band number length should not exceed 64 digits */
381 suffix = G_calloc(65, sizeof(char));
382 }
383
384 croptoregion = flag_r->answer;
385 if (flag_r->answer && parm.outloc->answer) {
386 G_warning(_("Disabling '-r' flag for new location"));
387 croptoregion = 0;
388 }
389
390 /* default GDAL memory cache size appears to be only 40 MiB, slowing down r.in.gdal */
391 if (parm.memory->answer && *parm.memory->answer) {
392 G_verbose_message(_("Using memory cache size: %s MiB"), parm.memory->answer);
393 CPLSetConfigOption("GDAL_CACHEMAX", parm.memory->answer);
394 }
395
396 /* GDAL configuration options */
397 if (parm.cfg->answer) {
398 char **tokens, *tok, *key, *value;
399 int i, ntokens;
400
401 tokens = G_tokenize(parm.cfg->answer, ",");
402 ntokens = G_number_of_tokens(tokens);
403 for (i = 0; i < ntokens; i++) {
404 G_debug(1, "%d=[%s]", i, tokens[i]);
405 tok = G_store(tokens[i]);
406 G_squeeze(tok);
407 key = tok;
408 value = strstr(tok, "=");
409 if (value) {
410 *value = '\0';
411 value++;
412 CPLSetConfigOption(key, value);
413 }
414 G_free(tok);
415 }
416 G_free_tokens(tokens);
417 }
418
419 /* GDAL dataset open options */
420 doo = NULL;
421 if (parm.doo->answer) {
422 char **tokens;
423 int i, ntokens;
424
425 tokens = G_tokenize(parm.doo->answer, ",");
426 ntokens = G_number_of_tokens(tokens);
427 doo = G_malloc(sizeof(char *) * (ntokens + 1));
428 for (i = 0; i < ntokens; i++) {
429 G_debug(1, "%d=[%s]", i, tokens[i]);
430 doo[i] = G_store(tokens[i]);
431 }
432 G_free_tokens(tokens);
433 doo[ntokens] = NULL;
434 }
435
436 /* -------------------------------------------------------------------- */
437 /* Open the dataset. */
438 /* -------------------------------------------------------------------- */
439
440 hDS = opends(input, (const char **)doo, &hDriver);
441
442 /* does the driver support subdatasets? */
443 /* test for capability GDAL_DMD_SUBDATASETS */
444
445 /* does the dataset include subdatasets? */
446 {
447 char **sds = GDALGetMetadata(hDS, "SUBDATASETS");
448
449 if (sds && *sds) {
450 int i;
451
452 G_warning(_("Input contains subdatasets which may need to "
453 "be imported separately by name:"));
454 /* list subdatasets */
455 for (i = 0; sds[i] != NULL; i++) {
456 char *sdsi = G_store(sds[i]);
457 char *p = sdsi;
458
459 while (*p != '=' && *p != '\0')
460 p++;
461 if (*p == '=')
462 p++;
463 if ((i & 1) == 0) {
464 fprintf(stderr, "Subdataset %d:\n", i + 1);
465 fprintf(stderr, " Name: %s\n", p);
466 }
467 else {
468 char *sdsdim, *sdsdesc, *sdstype;
469 int sdsdlen;
470
471 sdsdim = sdsdesc = sdstype = NULL;
472 sdsdlen = strlen(sdsi);
473 while (*p != '[' && *p != '\0')
474 p++;
475 if (*p == '[') {
476 p++;
477 sdsdim = p;
478 while (*p != ']' && *p != '\0')
479 p++;
480 if (*p == ']') {
481 *p = '\0';
482 p++;
483 }
484 }
485 while (*p == ' ' && *p != '\0')
486 p++;
487 sdsdesc = p;
488 if (*p != '\0') {
489 p = sdsi + sdsdlen - 1;
490
491 if (*p == ')') {
492 *p = '\0';
493 p--;
494 while (*p != '(')
495 p--;
496 if (*p == '(') {
497 sdstype = p + 1;
498 p--;
499 while (*p == ' ')
500 p--;
501 }
502 p++;
503 if (*p == ' ')
504 *p = '\0';
505 }
506 }
507 if (sdsdesc && *sdsdesc)
508 fprintf(stderr, " Description: %s\n", sdsdesc);
509 if (sdsdim && *sdsdim)
510 fprintf(stderr, " Dimension: %s\n", sdsdim);
511 if (sdstype && *sdstype)
512 fprintf(stderr, " Data type: %s\n", sdstype);
513 }
514 G_free(sdsi);
515 }
516 }
517 }
518
519 if (GDALGetRasterCount(hDS) == 0) {
520 G_fatal_error(_("No raster bands found in <%s>"), input);
521 }
522
523 if (flag_p->answer) {
524 /* print number of bands */
525 fprintf(stdout, "%d\n", GDALGetRasterCount(hDS));
526 GDALClose(hDS);
527 exit(EXIT_SUCCESS);
528 }
529
530 if (output && !parm.outloc->answer &&
531 GDALGetRasterCount(hDS) == 1) { /* Check if the map exists */
532 if (G_find_raster2(output, G_mapset())) {
533 if (overwrite)
534 G_warning(_("Raster map <%s> already exists and will be overwritten"),
535 output);
536 else
537 G_fatal_error(_("Raster map <%s> already exists"), output);
538 }
539 }
540
541 /* zero cell header */
542 G_zero(&cellhd, sizeof(struct Cell_head));
543
544 /* -------------------------------------------------------------------- */
545 /* Set up the window representing the data we have. */
546 /* -------------------------------------------------------------------- */
547 G_debug(3, "GDAL size: row = %d, col = %d", GDALGetRasterYSize(hDS),
548 GDALGetRasterXSize(hDS));
549 cellhd.rows = GDALGetRasterYSize(hDS);
550 cellhd.rows3 = GDALGetRasterYSize(hDS);
551 cellhd.cols = GDALGetRasterXSize(hDS);
552 cellhd.cols3 = GDALGetRasterXSize(hDS);
553
554 if (GDALGetGeoTransform(hDS, adfGeoTransform) == CE_None) {
555 if (adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 ||
556 adfGeoTransform[1] <= 0.0 || adfGeoTransform[5] >= 0.0) {
557 G_debug(0, "adfGeoTransform[2] %g", adfGeoTransform[2]);
558 G_debug(0, "adfGeoTransform[4] %g", adfGeoTransform[4]);
559 G_debug(0, "adfGeoTransform[1] %g", adfGeoTransform[1]);
560 G_debug(0, "adfGeoTransform[5] %g", adfGeoTransform[5]);
561 G_fatal_error(_("Input raster map is flipped or rotated - cannot import. "
562 "You may use 'gdalwarp' to transform the map to North-up."));
563 }
564 cellhd.north = adfGeoTransform[3];
565 cellhd.ns_res = fabs(adfGeoTransform[5]);
566 cellhd.ns_res3 = fabs(adfGeoTransform[5]);
567 cellhd.south = cellhd.north - cellhd.ns_res * cellhd.rows;
568
569 cellhd.west = adfGeoTransform[0];
570 cellhd.ew_res = fabs(adfGeoTransform[1]);
571 cellhd.ew_res3 = fabs(adfGeoTransform[1]);
572 cellhd.east = cellhd.west + cellhd.cols * cellhd.ew_res;
573
574 cellhd.top = 1.;
575 cellhd.bottom = 0.;
576 cellhd.tb_res = 1.;
577 cellhd.depths = 1;
578 }
579 else {
580 if (croptoregion) {
581 G_fatal_error(_("Unable to fetch the affine transformation coefficients. "
582 "Flag -%c cannot be used in this case."), flag_r->key);
583 }
584 cellhd.north = cellhd.rows;
585 cellhd.south = 0.0;
586 cellhd.ns_res = 1.0;
587 cellhd.ns_res3 = 1.0;
588 cellhd.west = 0.0;
589 cellhd.east = cellhd.cols;
590 cellhd.ew_res = 1.0;
591 cellhd.ew_res3 = 1.0;
592 cellhd.top = 1.;
593 cellhd.bottom = 0.;
594 cellhd.tb_res = 1.;
595 cellhd.depths = 1;
596 }
597
598 /* constrain to geographic coords */
599 if (flag_l->answer && G_projection() == PROJECTION_LL) {
600 if (cellhd.north > 90.) cellhd.north = 90.;
601 if (cellhd.south < -90.) cellhd.south = -90.;
602 if (cellhd.east > 360.) cellhd.east = 180.;
603 if (cellhd.west < -180.) cellhd.west = -180.;
604 cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
605 cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
606 cellhd.ew_res3 = cellhd.ew_res;
607 cellhd.ns_res3 = cellhd.ns_res;
608
609 G_warning(_("Map bounds have been constrained to geographic "
610 "coordinates. You will almost certainly want to check "
611 "map bounds and resolution with r.info and reset them "
612 "with r.region before going any further."));
613 }
614
615 check_projection(&cellhd, hDS, parm.outloc->answer, flag_c->answer,
616 flag_o->answer, flag_j->answer);
617
618 /* -------------------------------------------------------------------- */
619 /* Set the active window to match the available data. */
620 /* -------------------------------------------------------------------- */
621 if (flag_a->answer && cellhd.proj == PROJECTION_LL) {
622 G_adjust_Cell_head(&cellhd, 1, 1);
623 G_adjust_window_ll(&cellhd);
624 }
625
626 roff = coff = 0;
627 col_offset = 0;
628 if (croptoregion) {
629 int row, col, first, last;
630
631 Rast_get_window(&cur_wind);
632 Rast_align_window(&cur_wind, &cellhd);
633
634 rowmapall = G_malloc(sizeof(int) * cur_wind.rows);
635 colmapall = G_malloc(sizeof(int) * cur_wind.cols);
636
637 /* window mapping */
638 first = -1;
639 last = cur_wind.rows - 1;
640 for (row = 0; row < cur_wind.rows; row++) {
641 /* get GDAL row for cur wind row */
642 double north = Rast_row_to_northing(row + 0.5, &cur_wind);
643
644 rowmapall[row] = Rast_northing_to_row(north, &cellhd);
645 if (rowmapall[row] < 0 || rowmapall[row] >= cellhd.rows)
646 rowmapall[row] = -1;
647 else {
648 if (first < 0)
649 first = row;
650 last = row;
651 }
652 }
653 if (first == -1)
654 G_fatal_error(_("Input raster does not overlap current "
655 "computational region. Nothing to import."));
656
657 G_debug(1, "first row in cur wind %d, first row in source %d", first, rowmapall[first]);
658 rowmap = &rowmapall[first];
659 /* crop window */
660 if (first != 0 || last != cur_wind.rows - 1) {
661 G_debug(1, "Cropping NS extents");
662 cur_wind.north -= first * cur_wind.ns_res;
663 cur_wind.south += (cur_wind.rows - 1 - last) * cur_wind.ns_res;
664 cur_wind.rows = last - first + 1;
665 }
666 roff = (cellhd.north - cur_wind.north + cellhd.ns_res / 2.) / cellhd.ns_res;
667
668 first = -1;
669 last = cur_wind.cols - 1;
670 for (col = 0; col < cur_wind.cols; col++) {
671 /* get GDAL col for cur wind col */
672 double east = Rast_col_to_easting(col + 0.5, &cur_wind);
673
674 colmapall[col] = Rast_easting_to_col(east, &cellhd);
675 if (colmapall[col] < 0 || colmapall[col] >= cellhd.cols)
676 colmapall[col] = -1;
677 else {
678 if (first < 0)
679 first = col;
680 last = col;
681 }
682 }
683 if (first == -1)
684 G_fatal_error(_("Input raster does not overlap current "
685 "computational region. Nothing to import."));
686
687 G_debug(1, "first col in cur wind %d, first col in source %d", first, colmapall[first]);
688 col_offset = colmapall[first];
689 colmap = &colmapall[first];
690 /* crop window */
691 if (first != 0 || last != cur_wind.cols - 1) {
692 G_debug(1, "Cropping EW extents");
693 cur_wind.west += first * cur_wind.ew_res;
694 cur_wind.east -= (cur_wind.cols - 1 - last) * cur_wind.ew_res;
695 cur_wind.cols = last - first + 1;
696 }
697 coff = (cur_wind.west - cellhd.west + cellhd.ew_res / 2.) / cellhd.ew_res;
698
699 cellhd = cur_wind;
700 }
701 else {
702 int row, col;
703
704 rowmap = G_malloc(sizeof(int) * cellhd.rows);
705 colmap = G_malloc(sizeof(int) * cellhd.cols);
706
707 for (row = 0; row < cellhd.rows; row++)
708 rowmap[row] = row;
709 for (col = 0; col < cellhd.cols; col++)
710 colmap[col] = col;
711 }
712 Rast_set_window(&cellhd);
713
714 /* -------------------------------------------------------------------- */
715 /* Do we want to generate a simple raster, or an imagery group? */
716 /* -------------------------------------------------------------------- */
717 n_bands = 0;
718 if (parm.band->answer != NULL) {
719 while (parm.band->answers[n_bands])
720 n_bands++;
721 }
722
723 if (GDALGetRasterCount(hDS) > 1 && n_bands != 1) {
724 G_message(_("Importing %d raster bands..."),
725 (n_bands > 1 ? n_bands : GDALGetRasterCount(hDS)));
726 }
727
728 if ((GDALGetRasterCount(hDS) > 1 && n_bands != 1)
729 || GDALGetGCPCount(hDS) > 0)
730 force_imagery = TRUE;
731
732 /* -------------------------------------------------------------------- */
733 /* Simple case. Import a single band as a raster cell. */
734 /* -------------------------------------------------------------------- */
735 if (!force_imagery) {
736 int nBand = 1;
737
738 if (parm.band->answer != NULL)
739 nBand = atoi(parm.band->answers[0]);
740
741 hBand = GDALGetRasterBand(hDS, nBand);
742 if (hBand == NULL) {
743 G_fatal_error(_("Selected band (%d) does not exist"), nBand);
744 }
745
746 ImportBand(hBand, output, NULL, rowmap, colmap, col_offset);
747 if (parm.rat->answer)
748 dump_rat(hBand, parm.rat->answer, nBand);
749
750 if (title)
751 Rast_put_cell_title(output, title);
752 }
753
754 /* -------------------------------------------------------------------- */
755 /* Complete case, import a set of raster bands as an imagery */
756 /* group. */
757 /* -------------------------------------------------------------------- */
758 else {
759 struct Ref ref;
760 char szBandName[1024];
761 int nBand = 0;
762 char colornamebuf[512], colornamebuf2[512];
763 FILE *map_names_file = NULL;
764
765 if(parm.map_names_file->answer) {
766 map_names_file = fopen(parm.map_names_file->answer, "w");
767 if(map_names_file == NULL) {
768 G_fatal_error(_("Unable to open the map names output text file"));
769 }
770 }
771
772 I_init_group_ref(&ref);
773
774 colornamebuf2[0] = '\0';
775
776 n_bands = 0;
777 while (TRUE) {
778 if (parm.band->answer != NULL) {
779 if (parm.band->answers[n_bands] == NULL)
780 break;
781 nBand = atoi(parm.band->answers[n_bands++]);
782 }
783 else {
784 if (nBand >= GDALGetRasterCount(hDS))
785 break;
786 nBand++;
787 }
788
789 /* Generate the suffix */
790 if(num_digits > 0) {
791 G_snprintf(suffix, num_digits + 1, "%0*d", num_digits, nBand + offset);
792 } else {
793 G_snprintf(suffix, 65, "%d", nBand + offset);
794 }
795
796 G_debug(3, "Import raster band %d", nBand);
797 hBand = GDALGetRasterBand(hDS, nBand);
798 if (!hBand)
799 G_fatal_error(_("Unable to get raster band number %d"), nBand);
800 if (!flag_k->answer) {
801 /* use channel color names if present: */
802 strcpy(colornamebuf,
803 GDALGetColorInterpretationName
804 (GDALGetRasterColorInterpretation(hBand)));
805
806 /* check: two channels with identical name ? */
807 if (strcmp(colornamebuf, colornamebuf2) == 0)
808 sprintf(colornamebuf, "%s", suffix);
809 else
810 strcpy(colornamebuf2, colornamebuf);
811
812 /* avoid bad color names; in case of 'Gray' often all channels are named 'Gray' */
813 if (strcmp(colornamebuf, "Undefined") == 0 ||
814 strcmp(colornamebuf, "Gray") == 0)
815 sprintf(szBandName, "%s.%s", output, suffix);
816 else {
817 G_tolcase(colornamebuf);
818 sprintf(szBandName, "%s.%s", output, colornamebuf);
819 }
820 }
821 else
822 sprintf(szBandName, "%s.%s", output, suffix);
823
824 if (!parm.outloc->answer) { /* Check if the map exists */
825 if (G_find_raster2(szBandName, G_mapset())) {
826 if (overwrite)
827 G_warning(_("Raster map <%s> already exists and will be overwritten"),
828 szBandName);
829 else
830 G_fatal_error(_("Raster map <%s> already exists"), szBandName);
831 }
832 }
833
834 ImportBand(hBand, szBandName, &ref, rowmap, colmap, col_offset);
835
836 if(map_names_file)
837 fprintf(map_names_file, "%s\n", szBandName);
838
839 if (title)
840 Rast_put_cell_title(szBandName, title);
841 }
842
843 if(map_names_file)
844 fclose(map_names_file);
845
846 I_put_group_ref(output, &ref);
847 I_free_group_ref(&ref);
848
849 /* make this group the current group */
850 I_put_group(output);
851
852 /* -------------------------------------------------------------------- */
853 /* Output GCPs if present, we can only do this when writing an */
854 /* imagery group. */
855 /* -------------------------------------------------------------------- */
856 if (GDALGetGCPCount(hDS) > 0) {
857 struct Control_Points sPoints;
858 const GDAL_GCP *pasGCPs = GDALGetGCPs(hDS);
859 int iGCP;
860 struct pj_info iproj, /* input map proj parameters */
861 oproj, /* output map proj parameters */
862 tproj; /* transformation parameters */
863 int create_target;
864 struct Cell_head gcpcellhd;
865 double emin, emax, nmin, nmax;
866
867 G_zero(&gcpcellhd, sizeof(struct Cell_head));
868
869 sPoints.count = GDALGetGCPCount(hDS);
870 sPoints.e1 =
871 (double *)G_malloc(sizeof(double) * sPoints.count * 4);
872 sPoints.n1 = sPoints.e1 + sPoints.count;
873 sPoints.e2 = sPoints.e1 + 2 * sPoints.count;
874 sPoints.n2 = sPoints.e1 + 3 * sPoints.count;
875 sPoints.status = (int *)G_malloc(sizeof(int) * sPoints.count);
876
877 G_message(_("Copying %d GCPS in points file for <%s>"),
878 sPoints.count, output);
879 if (GDALGetGCPProjection(hDS) != NULL
880 && strlen(GDALGetGCPProjection(hDS)) > 0) {
881 G_message("%s\n"
882 "--------------------------------------------\n"
883 "%s\n"
884 "--------------------------------------------",
885 _("GCPs have the following OpenGIS WKT Coordinate System:"),
886 GDALGetGCPProjection(hDS));
887 }
888
889 create_target = 0;
890 if (parm.target->answer) {
891 char target_mapset[GMAPSET_MAX];
892
893 /* does the target location exist? */
894 G_create_alt_env();
895 G_setenv_nogisrc("LOCATION_NAME", parm.target->answer);
896 sprintf(target_mapset, "PERMANENT"); /* must exist */
897
898 if (G_mapset_permissions(target_mapset) == -1) {
899 /* create target location later */
900 create_target = 1;
901 }
902 G_switch_env();
903 }
904
905 if (parm.target->answer && !create_target) {
906 SetupReprojector(GDALGetGCPProjection(hDS),
907 parm.target->answer, &iproj, &oproj, &tproj);
908 G_message(_("Re-projecting GCPs table:"));
909 G_message(_("* Input projection for GCP table: %s"),
910 iproj.proj);
911 G_message(_("* Output projection for GCP table: %s"),
912 oproj.proj);
913 }
914
915 emin = emax = pasGCPs[0].dfGCPX;
916 nmin = nmax = pasGCPs[0].dfGCPY;
917
918 for (iGCP = 0; iGCP < sPoints.count; iGCP++) {
919 sPoints.e1[iGCP] = pasGCPs[iGCP].dfGCPPixel + coff;
920 /* GDAL lines from N to S -> GRASS Y from S to N */
921 sPoints.n1[iGCP] = cellhd.rows - pasGCPs[iGCP].dfGCPLine + roff;
922
923 sPoints.e2[iGCP] = pasGCPs[iGCP].dfGCPX; /* target */
924 sPoints.n2[iGCP] = pasGCPs[iGCP].dfGCPY;
925 sPoints.status[iGCP] = 1;
926
927 /* If desired, do GCPs transformation to other projection */
928 if (parm.target->answer) {
929 /* re-project target GCPs */
930 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
931 &(sPoints.e2[iGCP]),
932 &(sPoints.n2[iGCP]), NULL) < 0)
933 G_fatal_error(_("Error in %s (can't "
934 "re-project GCP %i)"),
935 "GPJ_transform()", iGCP);
936 }
937
938 /* figure out legal e, w, n, s values for new target location */
939 if (create_target) {
940 if (emin > pasGCPs[iGCP].dfGCPX)
941 emin = pasGCPs[iGCP].dfGCPX;
942 if (emax < pasGCPs[iGCP].dfGCPX)
943 emax = pasGCPs[iGCP].dfGCPX;
944 if (nmin > pasGCPs[iGCP].dfGCPY)
945 nmin = pasGCPs[iGCP].dfGCPY;
946 if (nmax < pasGCPs[iGCP].dfGCPY)
947 nmax = pasGCPs[iGCP].dfGCPY;
948 }
949 } /* for all GCPs */
950
951 I_put_control_points(output, &sPoints);
952 if (create_target) {
953 /* create target location */
954 if (GPJ_wkt_to_grass(&gcpcellhd, &proj_info,
955 &proj_units, GDALGetGCPProjection(hDS), 0) < 0) {
956 G_warning(_("Unable to convert input map projection to GRASS "
957 "format; cannot create new location."));
958 }
959 else {
960 gcpcellhd.west = emin;
961 gcpcellhd.east = emax;
962 gcpcellhd.south = nmin;
963 gcpcellhd.north = nmax;
964 gcpcellhd.rows = GDALGetRasterYSize(hDS);
965 gcpcellhd.cols = GDALGetRasterXSize(hDS);
966 gcpcellhd.ns_res = 1.0;
967 gcpcellhd.ns_res3 = 1.0;
968 gcpcellhd.ew_res = 1.0;
969 gcpcellhd.ew_res3 = 1.0;
970 gcpcellhd.top = 1.;
971 gcpcellhd.bottom = 0.;
972 gcpcellhd.tb_res = 1.;
973 gcpcellhd.depths = 1;
974
975 G_adjust_Cell_head(&gcpcellhd, 1, 1);
976
977 G_create_alt_env();
978 if (0 != G_make_location(parm.target->answer, &gcpcellhd,
979 proj_info, proj_units)) {
980 G_fatal_error(_("Unable to create new location <%s>"),
981 parm.target->answer);
982 }
983 /* switch back to import location */
984 G_switch_env();
985
986 G_message(_("Location <%s> created"), parm.target->answer);
987 /* set the group's target */
988 I_put_target(output, parm.target->answer, "PERMANENT");
989 G_message(_("The target for the output group <%s> has been set to "
990 "location <%s>, mapset <PERMANENT>."),
991 output, parm.target->answer);
992 }
993 }
994
995 G_free(sPoints.e1);
996 G_free(sPoints.status);
997 }
998 }
999
1000 /* close the GDALDataset to avoid segfault in libgdal */
1001 GDALClose(hDS);
1002
1003 /* -------------------------------------------------------------------- */
1004 /* Extend current window based on dataset. */
1005 /* -------------------------------------------------------------------- */
1006 if (flag_e->answer && !croptoregion) {
1007 if (strcmp(G_mapset(), "PERMANENT") == 0)
1008 /* fixme: expand WIND and DEFAULT_WIND independently. (currently
1009 WIND gets forgotten and DEFAULT_WIND is expanded for both) */
1010 G_get_default_window(&cur_wind);
1011 else
1012 G_get_window(&cur_wind);
1013
1014 cur_wind.north = MAX(cur_wind.north, cellhd.north);
1015 cur_wind.south = MIN(cur_wind.south, cellhd.south);
1016 cur_wind.west = MIN(cur_wind.west, cellhd.west);
1017 cur_wind.east = MAX(cur_wind.east, cellhd.east);
1018
1019 cur_wind.rows = (int)ceil((cur_wind.north - cur_wind.south)
1020 / cur_wind.ns_res);
1021 cur_wind.south = cur_wind.north - cur_wind.rows * cur_wind.ns_res;
1022
1023 cur_wind.cols = (int)ceil((cur_wind.east - cur_wind.west)
1024 / cur_wind.ew_res);
1025 cur_wind.east = cur_wind.west + cur_wind.cols * cur_wind.ew_res;
1026
1027 if (strcmp(G_mapset(), "PERMANENT") == 0) {
1028 G_put_element_window(&cur_wind, "", "DEFAULT_WIND");
1029 G_message(_("Default region for this location updated"));
1030 }
1031 G_put_window(&cur_wind);
1032 G_message(_("Region for the current mapset updated"));
1033 }
1034
1035 exit(EXIT_SUCCESS);
1036}
1037
1038/************************************************************************/
1039/* SetupReprojector() */
1040
1041/************************************************************************/
1042
1043static void SetupReprojector(const char *pszSrcWKT, const char *pszDstLoc,
1044 struct pj_info *iproj, struct pj_info *oproj,
1045 struct pj_info *tproj)
1046{
1047 struct Cell_head cellhd;
1048 struct Key_Value *proj_info = NULL, *proj_units = NULL;
1049 char errbuf[1024];
1050 int permissions;
1051 char target_mapset[GMAPSET_MAX];
1052 struct Key_Value *out_proj_info, /* projection information of */
1053 *out_unit_info; /* input and output mapsets */
1054
1055 /* -------------------------------------------------------------------- */
1056 /* Translate GCP WKT coordinate system into GRASS format. */
1057 /* -------------------------------------------------------------------- */
1058 GPJ_wkt_to_grass(&cellhd, &proj_info, &proj_units, pszSrcWKT, 0);
1059
1060 if (pj_get_kv(iproj, proj_info, proj_units) < 0)
1061 G_fatal_error(_("Unable to translate projection key values of input GCPs"));
1062
1063 /* -------------------------------------------------------------------- */
1064 /* Get the projection of the target location. */
1065 /* -------------------------------------------------------------------- */
1066
1067 /* Change to user defined target location for GCPs transformation */
1068 G_create_alt_env();
1069 G_setenv_nogisrc("LOCATION_NAME", (char *)pszDstLoc);
1070 sprintf(target_mapset, "PERMANENT"); /* to find PROJ_INFO */
1071
1072 permissions = G_mapset_permissions(target_mapset);
1073 if (permissions >= 0) {
1074
1075 /* Get projection info from target location */
1076 if ((out_proj_info = G_get_projinfo()) == NULL)
1077 G_fatal_error(_("Unable to get projection info of target location"));
1078 if ((out_unit_info = G_get_projunits()) == NULL)
1079 G_fatal_error(_("Unable to get projection units of target location"));
1080 if (pj_get_kv(oproj, out_proj_info, out_unit_info) < 0)
1081 G_fatal_error(_("Unable to get projection key values of target location"));
1082 tproj->def = NULL;
1083 if (GPJ_init_transform(iproj, oproj, tproj) < 0)
1084 G_fatal_error(_("Unable to initialize coordinate transformation"));
1085 }
1086 else { /* can't access target mapset */
1087 /* access to mapset PERMANENT in target location is not required */
1088 sprintf(errbuf, _("Mapset <%s> in target location <%s> - "),
1089 target_mapset, pszDstLoc);
1090 strcat(errbuf, permissions == 0 ? _("permission denied")
1091 : _("not found"));
1092 G_fatal_error("%s", errbuf);
1093 } /* permission check */
1094
1095 /* And switch back to original location */
1096 G_switch_env();
1097}
1098
1099
1100/************************************************************************/
1101/* ImportBand() */
1102
1103/************************************************************************/
1104
1105static void ImportBand(GDALRasterBandH hBand, const char *output,
1106 struct Ref *group_ref, int *rowmap, int *colmap,
1107 int col_offset)
1108{
1109 RASTER_MAP_TYPE data_type;
1110 GDALDataType eGDT, eRawGDT;
1111 int row, nrows, ncols, complex;
1112 int ncols_gdal, map_cols, use_cell_gdal;
1113 int cf, cfR, cfI, bNoDataEnabled;
1114 int indx;
1115 void *cell, *cellReal, *cellImg;
1116 void *cell_gdal;
1117 void *bufComplex;
1118 double dfNoData;
1119 char outputReal[GNAME_MAX], outputImg[GNAME_MAX];
1120 char *nullFlags = NULL;
1121 struct History history;
1122 char **GDALmetadata;
1123 int have_colors = 0;
1124 GDALRasterAttributeTableH gdal_rat;
1125
1126 G_message(_("Importing raster map <%s>..."), output);
1127
1128 /* -------------------------------------------------------------------- */
1129 /* Select a cell type for the new cell. */
1130 /* -------------------------------------------------------------------- */
1131 eRawGDT = GDALGetRasterDataType(hBand);
1132 complex = FALSE;
1133
1134 switch (eRawGDT) {
1135 case GDT_Float32:
1136 data_type = FCELL_TYPE;
1137 eGDT = GDT_Float32;
1138 break;
1139
1140 case GDT_Float64:
1141 data_type = DCELL_TYPE;
1142 eGDT = GDT_Float64;
1143 break;
1144
1145 case GDT_Byte:
1146 data_type = CELL_TYPE;
1147 eGDT = GDT_Int32;
1148 Rast_set_cell_format(0);
1149 break;
1150
1151 case GDT_Int16:
1152 case GDT_UInt16:
1153 data_type = CELL_TYPE;
1154 eGDT = GDT_Int32;
1155 Rast_set_cell_format(1);
1156 break;
1157
1158 /* TODO: complex */
1159
1160 default:
1161 data_type = CELL_TYPE;
1162 eGDT = GDT_Int32;
1163 Rast_set_cell_format(3);
1164 break;
1165 }
1166
1167 ncols = Rast_window_cols();
1168 ncols_gdal = GDALGetRasterBandXSize(hBand);
1169 nrows = Rast_window_rows();
1170
1171 /* -------------------------------------------------------------------- */
1172 /* Do we have a null value? */
1173 /* -------------------------------------------------------------------- */
1174 map_cols = 0;
1175 use_cell_gdal = 1;
1176
1177 for (indx = 0; indx < ncols; indx++) {
1178 if (indx != colmap[indx] - col_offset) {
1179 map_cols = 1;
1180 use_cell_gdal = 0;
1181 }
1182 if (colmap[indx] < 0) {
1183 nullFlags = (char *)G_malloc(sizeof(char) * ncols);
1184 memset(nullFlags, 0, ncols);
1185 map_cols = 1;
1186 use_cell_gdal = 0;
1187 break;
1188 }
1189 }
1190 G_debug(1, "need column mapping: %d", map_cols);
1191 G_debug(1, "use cell_gdal: %d", use_cell_gdal);
1192 dfNoData = GDALGetRasterNoDataValue(hBand, &bNoDataEnabled);
1193 if (bNoDataEnabled && !nullFlags) {
1194 nullFlags = (char *)G_malloc(sizeof(char) * ncols);
1195 memset(nullFlags, 0, ncols);
1196 }
1197
1198 /* -------------------------------------------------------------------- */
1199 /* Create the new raster(s) */
1200 /* -------------------------------------------------------------------- */
1201 if (complex) {
1202 sprintf(outputReal, "%s.real", output);
1203 cfR = Rast_open_new(outputReal, data_type);
1204 sprintf(outputImg, "%s.imaginary", output);
1205
1206 cfI = Rast_open_new(outputImg, data_type);
1207
1208 cellReal = Rast_allocate_buf(data_type);
1209 cellImg = Rast_allocate_buf(data_type);
1210 if (eGDT == GDT_Float64)
1211 bufComplex = (double *)G_malloc(sizeof(double) * ncols_gdal * 2);
1212 else
1213 bufComplex = (float *)G_malloc(sizeof(float) * ncols_gdal * 2);
1214
1215 if (group_ref != NULL) {
1216 I_add_file_to_group_ref(outputReal, G_mapset(), group_ref);
1217 I_add_file_to_group_ref(outputImg, G_mapset(), group_ref);
1218 }
1219 }
1220 else {
1221 cf = Rast_open_new(output, data_type);
1222
1223 if (group_ref != NULL)
1224 I_add_file_to_group_ref((char *)output, G_mapset(), group_ref);
1225
1226 cell_gdal = G_malloc(Rast_cell_size(data_type) * ncols_gdal);
1227 if (use_cell_gdal)
1228 cell = (char *)cell_gdal + Rast_cell_size(data_type) * col_offset;
1229 else
1230 cell = Rast_allocate_buf(data_type);
1231 }
1232
1233 /* -------------------------------------------------------------------- */
1234 /* Write the raster one scanline at a time. */
1235 /* We have to distinguish some cases due to the different */
1236 /* coordinate system orientation of GDAL and GRASS for xy data */
1237 /* -------------------------------------------------------------------- */
1238
1239 /* special cases first */
1240 if (complex) { /* CEOS SAR et al.: import flipped to match GRASS coordinates */
1241 for (row = 0; row < nrows; row++) {
1242 if (rowmap[row] < 0)
1243 G_fatal_error(_("Invalid row"));
1244
1245 G_percent(row, nrows, 2);
1246
1247 GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
1248 bufComplex, ncols_gdal, 1, eGDT, 0, 0);
1249
1250 for (indx = ncols - 1; indx >= 0; indx--) { /* CEOS: flip east-west during import - MN */
1251 if (eGDT == GDT_Int32) {
1252 if (colmap[indx] < 0) {
1253 Rast_set_c_null_value(&(((CELL *) cellReal)[ncols - indx]), 1);
1254 Rast_set_c_null_value(&(((CELL *) cellImg)[ncols - indx]), 1);
1255 }
1256 else {
1257 ((CELL *) cellReal)[ncols - indx] =
1258 ((GInt32 *) bufComplex)[colmap[indx] * 2];
1259 ((CELL *) cellImg)[ncols - indx] =
1260 ((GInt32 *) bufComplex)[colmap[indx] * 2 + 1];
1261 }
1262 }
1263 else if (eGDT == GDT_Float32) {
1264 if (colmap[indx] < 0) {
1265 Rast_set_f_null_value(&(((FCELL *) cellReal)[ncols - indx]), 1);
1266 Rast_set_f_null_value(&(((FCELL *) cellImg)[ncols - indx]), 1);
1267 }
1268 else {
1269 ((FCELL *)cellReal)[ncols - indx] =
1270 ((float *)bufComplex)[colmap[indx] * 2];
1271 ((FCELL *)cellImg)[ncols - indx] =
1272 ((float *)bufComplex)[colmap[indx] * 2 + 1];
1273 }
1274 }
1275 else if (eGDT == GDT_Float64) {
1276 if (colmap[indx] < 0) {
1277 Rast_set_d_null_value(&(((DCELL *) cellReal)[ncols - indx]), 1);
1278 Rast_set_d_null_value(&(((DCELL *) cellImg)[ncols - indx]), 1);
1279 }
1280 else {
1281 ((DCELL *)cellReal)[ncols - indx] =
1282 ((double *)bufComplex)[colmap[indx] * 2];
1283 ((DCELL *)cellImg)[ncols - indx] =
1284 ((double *)bufComplex)[colmap[indx] * 2 + 1];
1285 }
1286 }
1287 }
1288 Rast_put_row(cfR, cellReal, data_type);
1289 Rast_put_row(cfI, cellImg, data_type);
1290 }
1291 } /* end of complex */
1292 else if (l1bdriver) { /* AVHRR */
1293 /* AVHRR - read from south to north to match GCPs */
1294 /* AVHRR - as for other formats, read from north to south to match GCPs
1295 * MM 2013 with gdal 1.10 */
1296 for (row = 0; row < nrows; row++) {
1297 if (rowmap[row] < 0)
1298 G_fatal_error(_("Invalid row"));
1299
1300 G_percent(row, nrows, 2);
1301
1302 GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
1303 cell_gdal, ncols_gdal, 1, eGDT, 0, 0);
1304
1305 if (nullFlags != NULL) {
1306 memset(nullFlags, 0, ncols);
1307
1308 if (eGDT == GDT_Int32) {
1309 for (indx = 0; indx < ncols; indx++) {
1310 if (colmap[indx] < 0)
1311 nullFlags[indx] = 1;
1312 else if (bNoDataEnabled &&
1313 ((CELL *) cell_gdal)[colmap[indx]] == (GInt32) dfNoData) {
1314 nullFlags[indx] = 1;
1315 }
1316 else
1317 ((CELL *)cell)[indx] = ((CELL *)cell_gdal)[colmap[indx]];
1318 }
1319 }
1320 else if (eGDT == GDT_Float32) {
1321 for (indx = 0; indx < ncols; indx++) {
1322 if (colmap[indx] < 0)
1323 nullFlags[indx] = 1;
1324 else if (bNoDataEnabled &&
1325 ((FCELL *)cell_gdal)[colmap[indx]] == (float)dfNoData) {
1326 nullFlags[indx] = 1;
1327 }
1328 else
1329 ((FCELL *)cell)[indx] = ((FCELL *)cell_gdal)[colmap[indx]];
1330 }
1331 }
1332 else if (eGDT == GDT_Float64) {
1333 for (indx = 0; indx < ncols; indx++) {
1334 if (colmap[indx] < 0)
1335 nullFlags[indx] = 1;
1336 else if (bNoDataEnabled &&
1337 ((DCELL *)cell_gdal)[colmap[indx]] == dfNoData) {
1338 nullFlags[indx] = 1;
1339 }
1340 else
1341 ((DCELL *)cell)[indx] = ((DCELL *)cell_gdal)[colmap[indx]];
1342 }
1343 }
1344
1345 Rast_insert_null_values(cell, nullFlags, ncols, data_type);
1346 }
1347 else if (map_cols) {
1348 if (eGDT == GDT_Int32) {
1349 for (indx = 0; indx < ncols; indx++) {
1350 ((CELL *)cell)[indx] = ((CELL *)cell_gdal)[colmap[indx]];
1351 }
1352 }
1353 else if (eGDT == GDT_Float32) {
1354 for (indx = 0; indx < ncols; indx++) {
1355 ((FCELL *)cell)[indx] = ((FCELL *)cell_gdal)[colmap[indx]];
1356 }
1357 }
1358 else if (eGDT == GDT_Float64) {
1359 for (indx = 0; indx < ncols; indx++) {
1360 ((DCELL *)cell)[indx] = ((DCELL *)cell_gdal)[colmap[indx]];
1361 }
1362 }
1363 }
1364
1365 Rast_put_row(cf, cell, data_type);
1366 }
1367 } /* end AVHRR */
1368 else {
1369 /* default */
1370 for (row = 0; row < nrows; row++) {
1371 if (rowmap[row] < 0)
1372 G_fatal_error(_("Invalid row"));
1373
1374 G_percent(row, nrows, 2);
1375
1376 GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
1377 cell_gdal, ncols_gdal, 1, eGDT, 0, 0);
1378
1379 if (nullFlags != NULL) {
1380 memset(nullFlags, 0, ncols);
1381
1382 if (eGDT == GDT_Int32) {
1383 for (indx = 0; indx < ncols; indx++) {
1384 if (colmap[indx] < 0)
1385 nullFlags[indx] = 1;
1386 else if (bNoDataEnabled &&
1387 ((CELL *) cell_gdal)[colmap[indx]] == (GInt32) dfNoData) {
1388 nullFlags[indx] = 1;
1389 }
1390 else
1391 ((CELL *)cell)[indx] = ((CELL *)cell_gdal)[colmap[indx]];
1392 }
1393 }
1394 else if (eGDT == GDT_Float32) {
1395 for (indx = 0; indx < ncols; indx++) {
1396 if (colmap[indx] < 0)
1397 nullFlags[indx] = 1;
1398 else if (bNoDataEnabled &&
1399 ((FCELL *)cell_gdal)[colmap[indx]] == (float)dfNoData) {
1400 nullFlags[indx] = 1;
1401 }
1402 else
1403 ((FCELL *)cell)[indx] = ((FCELL *)cell_gdal)[colmap[indx]];
1404 }
1405 }
1406 else if (eGDT == GDT_Float64) {
1407 for (indx = 0; indx < ncols; indx++) {
1408 if (colmap[indx] < 0)
1409 nullFlags[indx] = 1;
1410 else if (bNoDataEnabled &&
1411 ((DCELL *)cell_gdal)[colmap[indx]] == dfNoData) {
1412 nullFlags[indx] = 1;
1413 }
1414 else
1415 ((DCELL *)cell)[indx] = ((DCELL *)cell_gdal)[colmap[indx]];
1416 }
1417 }
1418
1419 Rast_insert_null_values(cell, nullFlags, ncols, data_type);
1420 }
1421 else if (map_cols) {
1422 if (eGDT == GDT_Int32) {
1423 for (indx = 0; indx < ncols; indx++) {
1424 ((CELL *)cell)[indx] = ((CELL *)cell_gdal)[colmap[indx]];
1425 }
1426 }
1427 else if (eGDT == GDT_Float32) {
1428 for (indx = 0; indx < ncols; indx++) {
1429 ((FCELL *)cell)[indx] = ((FCELL *)cell_gdal)[colmap[indx]];
1430 }
1431 }
1432 else if (eGDT == GDT_Float64) {
1433 for (indx = 0; indx < ncols; indx++) {
1434 ((DCELL *)cell)[indx] = ((DCELL *)cell_gdal)[colmap[indx]];
1435 }
1436 }
1437 }
1438
1439 Rast_put_row(cf, cell, data_type);
1440 }
1441 }
1442 G_percent(1, 1, 1);
1443
1444 /* -------------------------------------------------------------------- */
1445 /* Cleanup */
1446 /* -------------------------------------------------------------------- */
1447 if (complex) {
1448 G_debug(1, "Creating support files for %s", outputReal);
1449 Rast_close(cfR);
1450 Rast_short_history((char *)outputReal, "raster", &history);
1451 Rast_command_history(&history);
1452 Rast_write_history((char *)outputReal, &history);
1453
1454 G_debug(1, "Creating support files for %s", outputImg);
1455 Rast_close(cfI);
1456 Rast_short_history((char *)outputImg, "raster", &history);
1457 Rast_command_history(&history);
1458 Rast_write_history((char *)outputImg, &history);
1459
1460 G_free(bufComplex);
1461 G_free(cellReal);
1462 G_free(cellImg);
1463 }
1464 else {
1465 G_debug(1, "Creating support files for %s", output);
1466 Rast_close(cf);
1467 Rast_short_history((char *)output, "raster", &history);
1468 Rast_command_history(&history);
1469 Rast_write_history((char *)output, &history);
1470
1471 if (!use_cell_gdal)
1472 G_free(cell);
1473 G_free(cell_gdal);
1474 }
1475
1476 if (nullFlags != NULL)
1477 G_free(nullFlags);
1478
1479 /* -------------------------------------------------------------------- */
1480 /* Transfer colormap, if there is one. */
1481 /* prefer color rules over color tables, search: */
1482 /* 1. GRASS color rules in metadata */
1483 /* 2. Raster attribute table with color rules */
1484 /* 3. Raster color table */
1485 /* -------------------------------------------------------------------- */
1486
1487 /* GRASS color rules in metadata? */
1488 GDALmetadata = GDALGetMetadata(hBand, "");
1489
1490 if (GDALmetadata) {
1491 struct Colors colors;
1492 DCELL val1, val2;
1493 int r1, g1, b1, r2, g2, b2;
1494
1495 Rast_init_colors(&colors);
1496
1497 while (GDALmetadata && GDALmetadata[0]) {
1498 G_debug(2, "%s", GDALmetadata[0]);
1499
1500 if (!strncmp("COLOR_TABLE_RULE_RGB_", GDALmetadata[0], 21)) {
1501 char *p;
1502
1503 for (p = GDALmetadata[0]; *p != '=' && *p != '\0'; p++);
1504
1505 if (*p == '=') {
1506 p++;
1507 }
1508 if (p && *p != '\0') {
1509 if (sscanf(p, "%lf %lf %d %d %d %d %d %d",
1510 &val1, &val2, &r1, &g1, &b1, &r2, &g2, &b2) == 8) {
1511
1512 Rast_add_d_color_rule(&val1, r1, g1, b1,
1513 &val2, r2, g2, b2,
1514 &colors);
1515 have_colors = 1;
1516 }
1517 }
1518 }
1519 GDALmetadata++;
1520 }
1521 if (have_colors)
1522 Rast_write_colors((char *)output, G_mapset(), &colors);
1523
1524 Rast_free_colors(&colors);
1525 }
1526
1527 /* colors in raster attribute table? */
1528 gdal_rat = GDALGetDefaultRAT(hBand);
1529 if (!have_colors && gdal_rat != NULL) {
1530 nrows = GDALRATGetRowCount(gdal_rat);
1531 ncols = GDALRATGetColumnCount(gdal_rat);
1532
1533 if (nrows > 0 && ncols > 0) {
1534 int minc, maxc, minmaxc;
1535 int rc, gc, bc, rminc, rmaxc, gminc, gmaxc, bminc, bmaxc;
1536 GDALRATFieldUsage field_use;
1537 struct Colors colors;
1538 DCELL val1, val2;
1539 double r1, g1, b1, r2, g2, b2;
1540 int cf;
1541
1542 Rast_init_colors(&colors);
1543
1544 minc = maxc = minmaxc = -1;
1545 rc = gc = bc = rminc = rmaxc = gminc = gmaxc = bminc = bmaxc = -1;
1546
1547 for (indx = 0; indx < ncols; indx++) {
1548 field_use = GDALRATGetUsageOfCol(gdal_rat, indx);
1549
1550 if (field_use == GFU_Min)
1551 minc = indx;
1552 else if (field_use == GFU_Max)
1553 maxc = indx;
1554 else if (field_use == GFU_MinMax)
1555 minmaxc = indx;
1556 else if (field_use == GFU_Red)
1557 rc = indx;
1558 else if (field_use == GFU_Green)
1559 gc = indx;
1560 else if (field_use == GFU_Blue)
1561 bc = indx;
1562 else if (field_use == GFU_RedMin)
1563 rminc = indx;
1564 else if (field_use == GFU_GreenMin)
1565 gminc = indx;
1566 else if (field_use == GFU_BlueMin)
1567 bminc = indx;
1568 else if (field_use == GFU_RedMax)
1569 rmaxc = indx;
1570 else if (field_use == GFU_GreenMax)
1571 gmaxc = indx;
1572 else if (field_use == GFU_BlueMax)
1573 bmaxc = indx;
1574 }
1575
1576 /* guess color range 0, 1 or 0, 255 */
1577
1578 if (minc >= 0 && maxc >= 0 && rminc >= 0 && rmaxc >= 0 &&
1579 gminc >= 0 && gmaxc >= 0 && bminc >= 0 && bmaxc >= 0) {
1580
1581 cf = 1;
1582
1583 /* analyze color rules */
1584 for (indx = 0; indx < nrows; indx++) {
1585 val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minc);
1586 val2 = GDALRATGetValueAsDouble(gdal_rat, indx, maxc);
1587
1588 r1 = GDALRATGetValueAsDouble(gdal_rat, indx, rminc);
1589 g1 = GDALRATGetValueAsDouble(gdal_rat, indx, gminc);
1590 b1 = GDALRATGetValueAsDouble(gdal_rat, indx, bminc);
1591
1592 r2 = GDALRATGetValueAsDouble(gdal_rat, indx, rmaxc);
1593 g2 = GDALRATGetValueAsDouble(gdal_rat, indx, gmaxc);
1594 b2 = GDALRATGetValueAsDouble(gdal_rat, indx, bmaxc);
1595
1596 if (r1 > 0.0 && r1 < 1.0)
1597 cf = 255;
1598 else if (cf == 255 && r1 > 1.0) {
1599 cf = 0;
1600 break;
1601 }
1602
1603 if (g1 > 0.0 && g1 < 1.0)
1604 cf = 255;
1605 else if (cf == 255 && g1 > 1.0) {
1606 cf = 0;
1607 break;
1608 }
1609
1610 if (b1 > 0.0 && b1 < 1.0)
1611 cf = 255;
1612 else if (cf == 255 && b1 > 1.0) {
1613 cf = 0;
1614 break;
1615 }
1616
1617 if (r2 > 0.0 && r2 < 1.0)
1618 cf = 255;
1619 else if (cf == 255 && r2 > 1.0) {
1620 cf = 0;
1621 break;
1622 }
1623
1624 if (g2 > 0.0 && g2 < 1.0)
1625 cf = 255;
1626 else if (cf == 255 && g2 > 1.0) {
1627 cf = 0;
1628 break;
1629 }
1630
1631 if (b2 > 0.0 && b2 < 1.0)
1632 cf = 255;
1633 else if (cf == 255 && b2 > 1.0) {
1634 cf = 0;
1635 break;
1636 }
1637 }
1638
1639 if (cf == 0)
1640 G_warning(_("Inconsistent color rules in RAT"));
1641 else {
1642 /* fetch color rules */
1643 for (indx = 0; indx < nrows; indx++) {
1644 val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minc);
1645 val2 = GDALRATGetValueAsDouble(gdal_rat, indx, maxc);
1646
1647 r1 = GDALRATGetValueAsDouble(gdal_rat, indx, rminc);
1648 g1 = GDALRATGetValueAsDouble(gdal_rat, indx, gminc);
1649 b1 = GDALRATGetValueAsDouble(gdal_rat, indx, bminc);
1650
1651 r2 = GDALRATGetValueAsDouble(gdal_rat, indx, rmaxc);
1652 g2 = GDALRATGetValueAsDouble(gdal_rat, indx, gmaxc);
1653 b2 = GDALRATGetValueAsDouble(gdal_rat, indx, bmaxc);
1654
1655 Rast_add_d_color_rule(&val1, r1 * cf, g1 * cf, b1 * cf,
1656 &val2, r2 * cf, g2 * cf, b2 * cf,
1657 &colors);
1658 }
1659 }
1660 }
1661 else if (minmaxc >= 0 && rc >= 0 && gc >= 0 && bc >= 0) {
1662
1663 cf = 1;
1664
1665 /* analyze color table */
1666 for (indx = 0; indx < nrows; indx++) {
1667 val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minmaxc);
1668
1669 r1 = GDALRATGetValueAsDouble(gdal_rat, indx, rc);
1670 g1 = GDALRATGetValueAsDouble(gdal_rat, indx, gc);
1671 b1 = GDALRATGetValueAsDouble(gdal_rat, indx, bc);
1672
1673
1674 if (r1 > 0.0 && r1 < 1.0)
1675 cf = 255;
1676 else if (cf == 255 && r1 > 1.0) {
1677 cf = 0;
1678 break;
1679 }
1680
1681 if (g1 > 0.0 && g1 < 1.0)
1682 cf = 255;
1683 else if (cf == 255 && g1 > 1.0) {
1684 cf = 0;
1685 break;
1686 }
1687
1688 if (b1 > 0.0 && b1 < 1.0)
1689 cf = 255;
1690 else if (cf == 255 && b1 > 1.0) {
1691 cf = 0;
1692 break;
1693 }
1694 }
1695
1696 if (cf == 0)
1697 G_warning(_("Inconsistent color rules in RAT"));
1698 else {
1699 /* fetch color table */
1700 for (indx = 0; indx < nrows; indx++) {
1701 val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minmaxc);
1702
1703 r1 = GDALRATGetValueAsDouble(gdal_rat, indx, rc);
1704 g1 = GDALRATGetValueAsDouble(gdal_rat, indx, gc);
1705 b1 = GDALRATGetValueAsDouble(gdal_rat, indx, bc);
1706
1707 Rast_set_d_color(val1, r1 * cf, g1 * cf, b1 * cf, &colors);
1708 }
1709 }
1710 }
1711
1712 have_colors = Rast_colors_count(&colors) > 0;
1713
1714 if (have_colors)
1715 Rast_write_colors((char *)output, G_mapset(), &colors);
1716
1717 Rast_free_colors(&colors);
1718 }
1719 }
1720
1721 /* colors in raster color table? */
1722
1723 if (!have_colors && !complex && GDALGetRasterColorTable(hBand) != NULL) {
1724 GDALColorTableH hCT;
1725 struct Colors colors;
1726 int iColor;
1727
1728 G_debug(1, "Copying color table for %s", output);
1729
1730 hCT = GDALGetRasterColorTable(hBand);
1731
1732 Rast_init_colors(&colors);
1733 for (iColor = 0; iColor < GDALGetColorEntryCount(hCT); iColor++) {
1734 GDALColorEntry sEntry;
1735
1736 GDALGetColorEntryAsRGB(hCT, iColor, &sEntry);
1737 if (sEntry.c4 == 0)
1738 continue;
1739
1740 Rast_set_c_color(iColor, sEntry.c1, sEntry.c2, sEntry.c3, &colors);
1741 }
1742
1743 Rast_write_colors((char *)output, G_mapset(), &colors);
1744 Rast_free_colors(&colors);
1745 have_colors = 1;
1746 }
1747 if (!have_colors) { /* no color table present */
1748
1749 /* types are defined in GDAL: ./core/gdal.h */
1750 if ((GDALGetRasterDataType(hBand) == GDT_Byte)) {
1751 /* found 0..255 data: we set to grey scale: */
1752 struct Colors colors;
1753
1754 G_verbose_message(_("Setting grey color table for <%s> (8bit, full range)"),
1755 output);
1756
1757 Rast_init_colors(&colors);
1758 Rast_make_grey_scale_colors(&colors, 0, 255); /* full range */
1759 Rast_write_colors((char *)output, G_mapset(), &colors);
1760 Rast_free_colors(&colors);
1761 }
1762 if ((GDALGetRasterDataType(hBand) == GDT_UInt16)) {
1763 /* found 0..65535 data: we set to grey scale: */
1764 struct Colors colors;
1765 struct Range range;
1766 CELL min, max;
1767
1768 G_verbose_message(_("Setting grey color table for <%s> (16bit, image range)"),
1769 output);
1770 Rast_read_range((char *)output, G_mapset(), &range);
1771 Rast_get_range_min_max(&range, &min, &max);
1772
1773 Rast_init_colors(&colors);
1774 Rast_make_grey_scale_colors(&colors, min, max); /* image range */
1775 Rast_write_colors((char *)output, G_mapset(), &colors);
1776 Rast_free_colors(&colors);
1777 }
1778 }
1779
1780 /* categories in raster attribute table? */
1781
1782 if (gdal_rat != NULL) {
1783 nrows = GDALRATGetRowCount(gdal_rat);
1784 ncols = GDALRATGetColumnCount(gdal_rat);
1785
1786 if (nrows > 0 && ncols > 0) {
1787 int minc, maxc, minmaxc, namec;
1788 GDALRATFieldUsage field_use;
1789 DCELL val1, val2;
1790 struct Categories cats;
1791 const char *label;
1792
1793 minc = maxc = minmaxc = namec = -1;
1794 for (indx = 0; indx < ncols; indx++) {
1795 field_use = GDALRATGetUsageOfCol(gdal_rat, indx);
1796
1797 if (field_use == GFU_Min)
1798 minc = indx;
1799 else if (field_use == GFU_Max)
1800 maxc = indx;
1801 else if (field_use == GFU_MinMax)
1802 minmaxc = indx;
1803 else if (field_use == GFU_Name)
1804 namec = indx;
1805 }
1806
1807 if (namec >= 0 && minmaxc >= 0) {
1808 Rast_init_cats("", &cats);
1809
1810 /* fetch labels */
1811 for (indx = 0; indx < nrows; indx++) {
1812 val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minmaxc);
1813 val2 = val1;
1814 label = GDALRATGetValueAsString(gdal_rat, indx, namec);
1815
1816 if (label)
1817 Rast_set_d_cat(&val1, &val2, label, &cats);
1818 }
1819 Rast_write_cats(output, &cats);
1820
1821 Rast_free_cats(&cats);
1822 }
1823 else if (namec >= 0 && minc >= 0 && maxc >= 0) {
1824 Rast_init_cats("", &cats);
1825
1826 /* fetch labels */
1827 for (indx = 0; indx < nrows; indx++) {
1828 val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minc);
1829 val2 = GDALRATGetValueAsDouble(gdal_rat, indx, maxc);
1830 label = GDALRATGetValueAsString(gdal_rat, indx, namec);
1831
1832 if (label)
1833 Rast_set_d_cat(&val1, &val2, label, &cats);
1834 }
1835 Rast_write_cats(output, &cats);
1836
1837 Rast_free_cats(&cats);
1838 }
1839 }
1840 }
1841
1842 return;
1843}
1844
1845static int dump_rat(GDALRasterBandH hBand, char *outrat, int nBand)
1846{
1847 int row, col, nrows, ncols;
1848 const char *field_name;
1849 GDALRATFieldUsage field_use;
1850 GDALRATFieldType *field_type;
1851 GDALRasterAttributeTableH gdal_rat;
1852 FILE *fp;
1853 char fname[GNAME_MAX];
1854
1855 if ((gdal_rat = GDALGetDefaultRAT(hBand)) == NULL)
1856 return 0;
1857
1858 nrows = GDALRATGetRowCount(gdal_rat);
1859 ncols = GDALRATGetColumnCount(gdal_rat);
1860
1861 if (nrows == 0 || ncols == 0)
1862 return 0;
1863
1864 field_type = G_malloc(ncols * sizeof(GDALRATFieldType));
1865
1866 G_snprintf(fname, GNAME_MAX, "%s_%d.csv", outrat, nBand);
1867 if (!(fp = fopen(fname, "w"))) {
1868 int err = errno;
1869
1870 G_fatal_error(_("Unable to open file <%s>: %s."),
1871 fname, strerror(err));
1872 }
1873
1874 /* dump column names and usage */
1875 for (col = 0; col < ncols; col++) {
1876
1877 if (col)
1878 fprintf(fp, "|");
1879
1880 field_name = GDALRATGetNameOfCol(gdal_rat, col);
1881 fprintf(fp, "%s", field_name);
1882 field_use = GDALRATGetUsageOfCol(gdal_rat, col);
1883
1884 if (field_use == GFU_Generic)
1885 fprintf(fp, " (General purpose field)");
1886 else if (field_use == GFU_PixelCount)
1887 fprintf(fp, " (Histogram pixel count)");
1888 else if (field_use == GFU_Name)
1889 fprintf(fp, " (Class name)");
1890 else if (field_use == GFU_Min)
1891 fprintf(fp, " (Class range minimum)");
1892 else if (field_use == GFU_Max)
1893 fprintf(fp, " (Class range maximum)");
1894 else if (field_use == GFU_MinMax)
1895 fprintf(fp, " (Class value (min=max))");
1896 else if (field_use == GFU_Red)
1897 fprintf(fp, " (Red class color (0-255))");
1898 else if (field_use == GFU_Green)
1899 fprintf(fp, " (Green class color (0-255))");
1900 else if (field_use == GFU_Blue)
1901 fprintf(fp, " (Blue class color (0-255))");
1902 else if (field_use == GFU_Alpha)
1903 fprintf(fp, " (Alpha (0=transparent,255=opaque))");
1904 else if (field_use == GFU_RedMin)
1905 fprintf(fp, " (Color Range Red Minimum)");
1906 else if (field_use == GFU_GreenMin)
1907 fprintf(fp, " (Color Range Green Minimum)");
1908 else if (field_use == GFU_BlueMin)
1909 fprintf(fp, " (Color Range Blue Minimum)");
1910 else if (field_use == GFU_AlphaMin)
1911 fprintf(fp, " (Color Range Alpha Minimum)");
1912 else if (field_use == GFU_RedMax)
1913 fprintf(fp, " (Color Range Red Maximum)");
1914 else if (field_use == GFU_GreenMax)
1915 fprintf(fp, " (Color Range Green Maximum)");
1916 else if (field_use == GFU_BlueMax)
1917 fprintf(fp, " (Color Range Blue Maximum)");
1918 else if (field_use == GFU_AlphaMax)
1919 fprintf(fp, " (Color Range Alpha Maximum)");
1920 else if (field_use == GFU_MaxCount)
1921 fprintf(fp, " (Maximum GFU value)");
1922 else
1923 fprintf(fp, " (Unknown)");
1924
1925 /* remember column type */
1926 field_type[col] = GDALRATGetTypeOfCol(gdal_rat, col);
1927 }
1928 fprintf(fp, "\n");
1929
1930 /* dump entries */
1931 for (row = 0; row < nrows; row++) {
1932
1933 for (col = 0; col < ncols; col++) {
1934 if (col)
1935 fprintf(fp, "|");
1936 if (field_type[col] == GFT_Integer)
1937 fprintf(fp, "%d", GDALRATGetValueAsInt(gdal_rat, row, col));
1938 else if (field_type[col] == GFT_Real)
1939 fprintf(fp, "%.15g", GDALRATGetValueAsDouble(gdal_rat, row, col));
1940 else
1941 fprintf(fp, "%s", GDALRATGetValueAsString(gdal_rat, row, col));
1942 }
1943 fprintf(fp, "\n");
1944 }
1945 fclose(fp);
1946
1947 return 1;
1948}
1949
1950void error_handler_ds(void *p)
1951{
1952 GDALDatasetH hDS = (GDALDatasetH) p;
1953 GDALClose(hDS);
1954}
Note: See TracBrowser for help on using the repository browser.