1
2 /****************************************************************************
3 *
4 * MODULE: r.out.gdal
5 * AUTHOR(S): Vytautas Vebra <olivership@gmail.com>
6 * PURPOSE: Exports GRASS raster to GDAL suported formats;
7 * based on GDAL library.
8 *
9 * COPYRIGHT: (C) 2006-2009 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * License (>=v2). Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 *****************************************************************************/
16
17 #include <grass/gis.h>
18 #include <grass/raster.h>
19 #include <grass/glocale.h>
20
21 #include <gdal.h>
22 #include <cpl_string.h>
23 #include <cpl_port.h>
24 #include "local_proto.h"
25
26 int exact_range_check(double, double, GDALDataType, const char *);
27
28 /* exact check for each band
29 * returns 0 on success
30 * -1 if given nodata value was present in data
31 * -2 if selected GDAL datatype could not hold all values
32 * */
exact_checks(GDALDataType export_datatype,const char * name,const char * mapset,struct Cell_head * cellhead,RASTER_MAP_TYPE maptype,double nodataval,const char * nodatakey,int default_nodataval)33 int exact_checks(GDALDataType export_datatype,
34 const char *name, const char *mapset,
35 struct Cell_head *cellhead, RASTER_MAP_TYPE maptype,
36 double nodataval, const char *nodatakey,
37 int default_nodataval)
38 {
39 double dfCellMin;
40 double dfCellMax;
41 int fd;
42 int cols = cellhead->cols;
43 int rows = cellhead->rows;
44 int ret = 0;
45
46 /* Open GRASS raster */
47 fd = Rast_open_old(name, mapset);
48
49 /* Create GRASS raster buffer */
50 void *bufer = Rast_allocate_buf(maptype);
51
52 if (bufer == NULL) {
53 G_warning(_("Unable to allocate buffer for reading raster map"));
54 return -1;
55 }
56
57 /* the following routine must be kept identical to export_band */
58
59 /* Copy data form GRASS raster to GDAL raster */
60 int row, col;
61 int n_nulls = 0, nodatavalmatch = 0;
62
63 dfCellMin = TYPE_FLOAT64_MAX;
64 dfCellMax = TYPE_FLOAT64_MIN;
65
66 /* Better use selected GDAL datatype instead of
67 * the best match with GRASS raster map types ? */
68
69 if (maptype == FCELL_TYPE) {
70
71 FCELL fnullval = (FCELL) nodataval;
72
73 G_debug(1, "FCELL nodata val: %f", fnullval);
74
75 for (row = 0; row < rows; row++) {
76
77 Rast_get_row(fd, bufer, row, maptype);
78 for (col = 0; col < cols; col++) {
79 FCELL fval = ((FCELL *) bufer)[col];
80
81 if (Rast_is_f_null_value(&fval)) {
82 ((FCELL *) bufer)[col] = fnullval;
83 n_nulls++;
84 }
85 else {
86 if (fval == fnullval) {
87 nodatavalmatch = 1;
88 }
89 if (!CPLIsInf(fval)) {
90 /* ignore inf */
91 if (dfCellMin > fval)
92 dfCellMin = fval;
93 if (dfCellMax < fval)
94 dfCellMax = fval;
95 }
96 }
97 }
98 G_percent(row + 1, rows, 2);
99 }
100 }
101 else if (maptype == DCELL_TYPE) {
102
103 DCELL dnullval = (DCELL) nodataval;
104
105 G_debug(1, "DCELL nodata val: %f", dnullval);
106
107 for (row = 0; row < rows; row++) {
108
109 Rast_get_row(fd, bufer, row, maptype);
110 for (col = 0; col < cols; col++) {
111 DCELL dval = ((DCELL *) bufer)[col];
112
113 if (Rast_is_d_null_value(&dval)) {
114 ((DCELL *) bufer)[col] = dnullval;
115 n_nulls++;
116 }
117 else {
118 if (dval == dnullval) {
119 nodatavalmatch = 1;
120 }
121 if (!CPLIsInf(dval)) {
122 /* ignore inf */
123 if (dfCellMin > dval)
124 dfCellMin = dval;
125 if (dfCellMax < dval)
126 dfCellMax = dval;
127 }
128 }
129 }
130 G_percent(row + 1, rows, 2);
131 }
132 }
133 else {
134
135 CELL inullval = (CELL) nodataval;
136
137 G_debug(1, "CELL nodata val: %d", inullval);
138
139 for (row = 0; row < rows; row++) {
140
141 Rast_get_row(fd, bufer, row, maptype);
142 for (col = 0; col < cols; col++) {
143 if (Rast_is_c_null_value(&((CELL *) bufer)[col])) {
144 ((CELL *) bufer)[col] = inullval;
145 n_nulls++;
146 }
147 else {
148 if (((CELL *) bufer)[col] == inullval) {
149 nodatavalmatch = 1;
150 }
151 if (dfCellMin > ((CELL *) bufer)[col])
152 dfCellMin = ((CELL *) bufer)[col];
153 if (dfCellMax < ((CELL *) bufer)[col])
154 dfCellMax = ((CELL *) bufer)[col];
155 }
156 }
157 G_percent(row + 1, rows, 2);
158 }
159 }
160 G_debug(1, "min %g max %g", dfCellMin, dfCellMax);
161
162 /* can the GDAL datatype hold the data range to be exported ? */
163 /* f-flag does not override */
164 if (exact_range_check(dfCellMin, dfCellMax, export_datatype, name)) {
165 G_warning("Raster export results in data loss.");
166 ret = -2;
167 }
168 G_message(_("Using GDAL data type <%s>"), GDALGetDataTypeName(export_datatype));
169
170 /* a default nodata value was used and NULL cells were present */
171 if (n_nulls && default_nodataval) {
172 if (maptype == CELL_TYPE)
173 G_important_message(_("Input raster map contains cells with NULL-value (no-data). "
174 "The value %d will be used to represent no-data values in the input map. "
175 "You can specify a nodata value with the %s option."),
176 (int)nodataval, nodatakey);
177 else
178 G_important_message(_("Input raster map contains cells with NULL-value (no-data). "
179 "The value %g will be used to represent no-data values in the input map. "
180 "You can specify a nodata value with the %s option."),
181 nodataval, nodatakey);
182 }
183
184 /* the nodata value was present in the exported data */
185 if (nodatavalmatch && n_nulls) {
186 /* default nodataval didn't work */
187 if (default_nodataval) {
188 G_warning(_("The default nodata value is present in raster"
189 "band <%s> and would lead to data loss. Please specify a "
190 "custom nodata value with the %s parameter."),
191 name, nodatakey);
192 }
193 /* user-specified nodataval didn't work */
194 else {
195 G_warning(_("The user given nodata value %g is present in raster"
196 "band <%s> and would lead to data loss. Please specify a "
197 "different nodata value with the %s parameter."),
198 nodataval, name, nodatakey);
199 }
200 ret = -1;
201 }
202
203 Rast_close(fd);
204
205 G_free(bufer);
206
207 return ret;
208 }
209
210 /* actual raster band export
211 * returns 0 on success
212 * -1 on raster data read/write error
213 * */
export_band(GDALDatasetH hMEMDS,int band,const char * name,const char * mapset,struct Cell_head * cellhead,RASTER_MAP_TYPE maptype,double nodataval,int suppress_main_colortable,int no_metadata,int writenodata)214 int export_band(GDALDatasetH hMEMDS, int band,
215 const char *name, const char *mapset,
216 struct Cell_head *cellhead, RASTER_MAP_TYPE maptype,
217 double nodataval, int suppress_main_colortable,
218 int no_metadata, int writenodata)
219 {
220 struct Colors sGrassColors;
221 GDALColorTableH hCT;
222 int iColor;
223 int bHaveMinMax;
224 double dfCellMin;
225 double dfCellMax;
226 struct FPRange sRange;
227 int fd;
228 int cols = cellhead->cols;
229 int rows = cellhead->rows;
230 int ret = 0;
231 char value[200];
232 char *units;
233
234 /* Open GRASS raster */
235 fd = Rast_open_old(name, mapset);
236
237 /* Get raster band */
238 GDALRasterBandH hBand = GDALGetRasterBand(hMEMDS, band);
239
240 if (hBand == NULL) {
241 G_warning(_("Unable to get raster band"));
242 return -1;
243 }
244
245 if (!no_metadata)
246 GDALSetDescription(hBand, name);
247
248 units = Rast_read_units(name, mapset);
249 if (units)
250 GDALSetRasterUnitType(hBand, units);
251
252 /* Get min/max values. */
253 if (Rast_read_fp_range(name, mapset, &sRange) == -1) {
254 bHaveMinMax = FALSE;
255 }
256 else {
257 bHaveMinMax = TRUE;
258 Rast_get_fp_range_min_max(&sRange, &dfCellMin, &dfCellMax);
259 }
260
261 /* use default color rules if no color rules are given */
262 if (Rast_read_colors(name, mapset, &sGrassColors) >= 0) {
263 int maxcolor, i;
264 CELL min, max;
265 char key[200];
266 int rcount;
267
268 Rast_get_c_color_range(&min, &max, &sGrassColors);
269 if (bHaveMinMax) {
270 if (max < dfCellMax) {
271 maxcolor = max;
272 }
273 else {
274 maxcolor = (int)ceil(dfCellMax);
275 }
276 if (maxcolor > GRASS_MAX_COLORS) {
277 maxcolor = GRASS_MAX_COLORS;
278 G_warning("Too many values, color table cut to %d entries",
279 maxcolor);
280 }
281 }
282 else {
283 if (max < GRASS_MAX_COLORS) {
284 maxcolor = max;
285 }
286 else {
287 maxcolor = GRASS_MAX_COLORS;
288 G_warning("Too many values, color table set to %d entries",
289 maxcolor);
290 }
291 }
292
293 rcount = Rast_colors_count(&sGrassColors);
294
295 G_debug(3, "dfCellMin: %f, dfCellMax: %f, maxcolor: %d", dfCellMin,
296 dfCellMax, maxcolor);
297
298 if (!suppress_main_colortable) {
299 hCT = GDALCreateColorTable(GPI_RGB);
300
301 for (iColor = 0; iColor <= maxcolor; iColor++) {
302 int nRed, nGreen, nBlue;
303 GDALColorEntry sColor;
304
305 if (Rast_get_c_color(&iColor, &nRed, &nGreen, &nBlue,
306 &sGrassColors)) {
307 sColor.c1 = nRed;
308 sColor.c2 = nGreen;
309 sColor.c3 = nBlue;
310 sColor.c4 = 255;
311
312 G_debug(3,
313 "Rast_get_c_color: Y, rcount %d, nRed %d, nGreen %d, nBlue %d",
314 rcount, nRed, nGreen, nBlue);
315 GDALSetColorEntry(hCT, iColor, &sColor);
316 }
317 else {
318 sColor.c1 = 0;
319 sColor.c2 = 0;
320 sColor.c3 = 0;
321 sColor.c4 = 0;
322
323 G_debug(3,
324 "Rast_get_c_color: N, rcount %d, nRed %d, nGreen %d, nBlue %d",
325 rcount, nRed, nGreen, nBlue);
326 GDALSetColorEntry(hCT, iColor, &sColor);
327 }
328 }
329
330 GDALSetRasterColorTable(hBand, hCT);
331 }
332
333 if (!no_metadata) {
334 if (rcount > 0) {
335 /* Create metadata entries for color table rules */
336 sprintf(value, "%d", rcount);
337 GDALSetMetadataItem(hBand, "COLOR_TABLE_RULES_COUNT", value,
338 NULL);
339 }
340
341 /* Add the rules in reverse order */
342 /* This can cause a GDAL warning with many rules, something like
343 * Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in tag. */
344 for (i = rcount - 1; i >= 0; i--) {
345 DCELL val1, val2;
346 unsigned char r1, g1, b1, r2, g2, b2;
347
348 Rast_get_fp_color_rule(&val1, &r1, &g1, &b1, &val2, &r2, &g2, &b2,
349 &sGrassColors, i);
350
351
352 sprintf(key, "COLOR_TABLE_RULE_RGB_%d", rcount - i - 1);
353 sprintf(value, "%e %e %d %d %d %d %d %d", val1, val2, r1, g1, b1,
354 r2, g2, b2);
355 GDALSetMetadataItem(hBand, key, value, NULL);
356 }
357 }
358 }
359
360 /* Create GRASS raster buffer */
361 void *bufer = Rast_allocate_buf(maptype);
362
363 if (bufer == NULL) {
364 G_warning(_("Unable to allocate buffer for reading raster map"));
365 return -1;
366 }
367
368 /* the following routine must be kept identical to exact_checks */
369
370 /* Copy data form GRASS raster to GDAL raster */
371 int row, col;
372 int n_nulls = 0;
373
374 /* Better use selected GDAL datatype instead of
375 * the best match with GRASS raster map types ? */
376
377 if (maptype == FCELL_TYPE) {
378
379 /* Source datatype understandable by GDAL */
380 GDALDataType datatype = GDT_Float32;
381 FCELL fnullval = (FCELL) nodataval;
382
383 G_debug(1, "FCELL nodata val: %f", fnullval);
384
385 for (row = 0; row < rows; row++) {
386
387 Rast_get_row(fd, bufer, row, maptype);
388 for (col = 0; col < cols; col++) {
389 if (Rast_is_f_null_value(&((FCELL *) bufer)[col])) {
390 ((FCELL *) bufer)[col] = fnullval;
391 if (n_nulls == 0) {
392 GDALSetRasterNoDataValue(hBand, nodataval);
393 }
394 n_nulls++;
395 }
396 }
397
398 if (GDALRasterIO
399 (hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
400 0, 0) >= CE_Failure) {
401 G_warning(_("Unable to write GDAL raster file"));
402 return -1;
403 }
404 G_percent(row + 1, rows, 2);
405 }
406 }
407 else if (maptype == DCELL_TYPE) {
408
409 GDALDataType datatype = GDT_Float64;
410 DCELL dnullval = (DCELL) nodataval;
411
412 G_debug(1, "DCELL nodata val: %f", dnullval);
413
414 for (row = 0; row < rows; row++) {
415
416 Rast_get_row(fd, bufer, row, maptype);
417 for (col = 0; col < cols; col++) {
418 if (Rast_is_d_null_value(&((DCELL *) bufer)[col])) {
419 ((DCELL *) bufer)[col] = dnullval;
420 if (n_nulls == 0) {
421 GDALSetRasterNoDataValue(hBand, nodataval);
422 }
423 n_nulls++;
424 }
425 }
426
427 if (GDALRasterIO
428 (hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
429 0, 0) >= CE_Failure) {
430 G_warning(_("Unable to write GDAL raster file"));
431 return -1;
432 }
433 G_percent(row + 1, rows, 2);
434 }
435 }
436 else {
437
438 GDALDataType datatype = GDT_Int32;
439 CELL inullval = (CELL) nodataval;
440
441 G_debug(1, "CELL nodata val: %d", inullval);
442
443 for (row = 0; row < rows; row++) {
444
445 Rast_get_row(fd, bufer, row, maptype);
446 for (col = 0; col < cols; col++) {
447 if (Rast_is_c_null_value(&((CELL *) bufer)[col])) {
448 ((CELL *) bufer)[col] = inullval;
449 if (n_nulls == 0) {
450 GDALSetRasterNoDataValue(hBand, nodataval);
451 }
452 n_nulls++;
453 }
454 }
455
456 if (GDALRasterIO
457 (hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
458 0, 0) >= CE_Failure) {
459 G_warning(_("Unable to write GDAL raster file"));
460 return -1;
461 }
462 G_percent(row + 1, rows, 2);
463 }
464 }
465 if (writenodata && n_nulls == 0)
466 GDALSetRasterNoDataValue(hBand, nodataval);
467
468 Rast_close(fd);
469
470 G_free(bufer);
471
472 return ret;
473 }
474
exact_range_check(double min,double max,GDALDataType datatype,const char * name)475 int exact_range_check(double min, double max, GDALDataType datatype,
476 const char *name)
477 {
478
479 switch (datatype) {
480 case GDT_Byte:
481 if (min < TYPE_BYTE_MIN || max > TYPE_BYTE_MAX) {
482 G_warning(_("Selected GDAL datatype does not cover data range."));
483 G_warning(_("GDAL datatype: %s, range: %d - %d"),
484 GDALGetDataTypeName(datatype), TYPE_BYTE_MIN,
485 TYPE_BYTE_MAX);
486 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
487 return 1;
488 }
489 else
490 return 0;
491
492 case GDT_UInt16:
493 if (min < TYPE_UINT16_MIN || max > TYPE_UINT16_MAX) {
494 G_warning(_("Selected GDAL datatype does not cover data range."));
495 G_warning(_("GDAL datatype: %s, range: %d - %d"),
496 GDALGetDataTypeName(datatype), TYPE_UINT16_MIN,
497 TYPE_UINT16_MAX);
498 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
499 return 1;
500 }
501 else
502 return 0;
503
504 case GDT_Int16:
505 case GDT_CInt16:
506 if (min < TYPE_INT16_MIN || max > TYPE_INT16_MAX) {
507 G_warning(_("Selected GDAL datatype does not cover data range."));
508 G_warning(_("GDAL datatype: %s, range: %d - %d"),
509 GDALGetDataTypeName(datatype), TYPE_INT16_MIN,
510 TYPE_INT16_MAX);
511 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
512 return 1;
513 }
514 else
515 return 0;
516
517 case GDT_Int32:
518 case GDT_CInt32:
519 if (min < TYPE_INT32_MIN || max > TYPE_INT32_MAX) {
520 G_warning(_("Selected GDAL datatype does not cover data range."));
521 G_warning(_("GDAL datatype: %s, range: %d - %d"),
522 GDALGetDataTypeName(datatype), TYPE_INT32_MIN,
523 TYPE_INT32_MAX);
524 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
525 return 1;
526 }
527 else
528 return 0;
529
530 case GDT_UInt32:
531 if (min < TYPE_UINT32_MIN || max > TYPE_UINT32_MAX) {
532 G_warning(_("Selected GDAL datatype does not cover data range."));
533 G_warning(_("GDAL datatype: %s, range: %u - %u"),
534 GDALGetDataTypeName(datatype), TYPE_UINT32_MIN,
535 TYPE_UINT32_MAX);
536 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
537 return 1;
538 }
539 else
540 return 0;
541
542 case GDT_Float32:
543 case GDT_CFloat32:
544 /* support export of inf / -inf ? */
545 if ((float)min !=TYPE_FLOAT32_MIN && (float)max != TYPE_FLOAT32_MAX &&
546 (min < TYPE_FLOAT32_MIN || max > TYPE_FLOAT32_MAX)) {
547 G_warning(_("Selected GDAL datatype does not cover data range."));
548 G_warning(_("GDAL datatype: %s, range: %.7g - %.7g"),
549 GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
550 TYPE_FLOAT32_MAX);
551 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
552 return 1;
553 }
554 else
555 return 0;
556
557 case GDT_Float64:
558 case GDT_CFloat64:
559 /* not possible because DCELL is FLOAT64, not 128bit floating point, but anyway... */
560 /* support export of inf / -inf ? */
561 if (min < TYPE_FLOAT64_MIN || max > TYPE_FLOAT64_MAX) {
562 G_warning(_("Selected GDAL datatype does not cover data range."));
563 G_warning(_("GDAL datatype: %s, range: %.16g - %.16g"),
564 GDALGetDataTypeName(datatype), TYPE_FLOAT64_MIN,
565 TYPE_FLOAT64_MAX);
566 G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
567 return 1;
568 }
569 else
570 return 0;
571
572 default:
573 return 0;
574 }
575 }
576