GRASS GIS 8 Programmer's Manual  8.2.1dev(2022)-5aae26b9b
gdal.c
Go to the documentation of this file.
1 /*!
2  \file lib/raster/gdal.c
3 
4  \brief Raster Library - Utilization of GDAL library.
5 
6  (C) 2010 by the GRASS Development Team
7 
8  This program is free software under the GNU General Public License
9  (>=v2). Read the file COPYING that comes with GRASS for details.
10 
11  \author Glynn Clements
12 */
13 
14 #include <stdlib.h>
15 #include <string.h>
16 #include <unistd.h>
17 
18 #include <grass/config.h>
19 #include <grass/gis.h>
20 #include <grass/raster.h>
21 #include <grass/gprojects.h>
22 #include <grass/glocale.h>
23 
24 #include "R.h"
25 
26 #ifndef HAVE_GDAL
27 #undef GDAL_LINK
28 #endif
29 
30 #ifdef GDAL_LINK
31 #include <gdal.h>
32 #endif
33 
34 /*!
35  \brief Initialization
36 
37  Register all GDAL drivers.
38 */
39 void Rast_init_gdal(void)
40 {
41 #ifdef GDAL_LINK
42  static int initialized;
43 
44  if (G_is_initialized(&initialized))
45  return;
46 
47  GDALAllRegister();
48  G_initialize_done(&initialized);
49 #endif
50 }
51 
52 /*!
53  \brief Get GDAL link settings for given raster map
54 
55  \param name map name
56  \param mapset name of mapset
57 
58  \return pointer to GDAL_link structure
59  \return NULL if link not found
60 */
61 struct GDAL_link *Rast_get_gdal_link(const char *name, const char *mapset)
62 {
63 #ifdef GDAL_LINK
64  GDALDatasetH data;
65  GDALRasterBandH band;
66  GDALDataType type;
67  RASTER_MAP_TYPE req_type;
68 #endif
69  const char *filename;
70  int band_num;
71  struct GDAL_link *gdal;
72  RASTER_MAP_TYPE map_type;
73  FILE *fp;
74  struct Key_Value *key_val;
75  const char *p;
76  DCELL null_val;
77  int hflip, vflip;
78 
79  if (!G_find_raster2(name, mapset))
80  return NULL;
81 
82  map_type = Rast_map_type(name, mapset);
83  if (map_type < 0)
84  return NULL;
85 
86  fp = G_fopen_old_misc("cell_misc", "gdal", name, mapset);
87  if (!fp)
88  return NULL;
89  key_val = G_fread_key_value(fp);
90  fclose(fp);
91 
92  if (!key_val)
93  return NULL;
94 
95  filename = G_find_key_value("file", key_val);
96  if (!filename)
97  return NULL;
98 
99  p = G_find_key_value("band", key_val);
100  if (!p)
101  return NULL;
102  band_num = atoi(p);
103  if (!band_num)
104  return NULL;
105 
106  p = G_find_key_value("null", key_val);
107  if (!p)
108  return NULL;
109  /* atof on windows can not read "nan" and returns 0 instead */
110  if (strcmp(p, "none") == 0 ||
111  G_strcasecmp(p, "nan") == 0 || G_strcasecmp(p, "-nan") == 0) {
112  Rast_set_d_null_value(&null_val, 1);
113  }
114  else
115  null_val = atof(p);
116 
117  hflip = G_find_key_value("hflip", key_val) ? 1 : 0;
118  vflip = G_find_key_value("vflip", key_val) ? 1 : 0;
119 
120 #ifdef GDAL_LINK
121  p = G_find_key_value("type", key_val);
122  if (!p)
123  return NULL;
124  type = atoi(p);
125 
126  switch (type) {
127  case GDT_Byte:
128  case GDT_Int16:
129  case GDT_UInt16:
130  case GDT_Int32:
131  case GDT_UInt32:
132  req_type = CELL_TYPE;
133  break;
134  case GDT_Float32:
135  req_type = FCELL_TYPE;
136  break;
137  case GDT_Float64:
138  req_type = DCELL_TYPE;
139  break;
140  default:
141  return NULL;
142  }
143 
144  if (req_type != map_type)
145  return NULL;
146 
147  Rast_init_gdal();
148 
149  data = GDALOpen(filename, GA_ReadOnly);
150  if (!data)
151  return NULL;
152 
153  band = GDALGetRasterBand(data, band_num);
154  if (!band) {
155  GDALClose(data);
156  return NULL;
157  }
158 #endif
159 
160  gdal = G_calloc(1, sizeof(struct GDAL_link));
161 
162  gdal->filename = G_store(filename);
163  gdal->band_num = band_num;
164  gdal->null_val = null_val;
165  gdal->hflip = hflip;
166  gdal->vflip = vflip;
167 #ifdef GDAL_LINK
168  gdal->data = data;
169  gdal->band = band;
170  gdal->type = type;
171 #endif
172 
173  return gdal;
174 }
175 
176 struct GDAL_Options
177 {
178  const char *dir;
179  const char *ext;
180  const char *format;
181  char **options;
182 };
183 
184 static struct state
185 {
186  int initialized;
187  struct GDAL_Options opts;
188  struct Key_Value *projinfo, *projunits, *projepsg;
189  char *srswkt;
190 } state;
191 
192 static struct state *st = &state;
193 
194 static void read_gdal_options(void)
195 {
196  FILE *fp;
197  struct Key_Value *key_val;
198  const char *p;
199 
200  fp = G_fopen_old("", "GDAL", G_mapset());
201  if (!fp)
202  G_fatal_error(_("Unable to open GDAL file"));
203  key_val = G_fread_key_value(fp);
204  fclose(fp);
205 
206  p = G_find_key_value("directory", key_val);
207  if (!p)
208  p = "gdal";
209  if (*p == '/') {
210  st->opts.dir = G_store(p);
211  }
212  else {
213  char path[GPATH_MAX];
214 
215  G_file_name(path, p, "", G_mapset());
216  st->opts.dir = G_store(path);
217  if (access(path, 0) != 0)
219  }
220 
221  p = G_find_key_value("extension", key_val);
222  st->opts.ext = G_store(p ? p : "");
223 
224  p = G_find_key_value("format", key_val);
225  st->opts.format = G_store(p ? p : "GTiff");
226 
227  p = G_find_key_value("options", key_val);
228  st->opts.options = p ? G_tokenize(p, ",") : NULL;
229 
230  G_free_key_value(key_val);
231 }
232 
233 /*!
234  \brief Create GDAL settings for given raster map
235 
236  \param name map name
237  \param map_type map type (CELL, FCELL, DCELL)
238 
239  \return pointer to allocated GDAL_link structure
240  \return NULL on error
241 */
243  RASTER_MAP_TYPE map_type)
244 {
245 #ifdef GDAL_LINK
246  char path[GPATH_MAX];
247  GDALDriverH driver;
248  double transform[6];
249  struct GDAL_link *gdal;
250  FILE *fp;
251  struct Key_Value *key_val;
252  char buf[32];
253 
255 
256  Rast_init_gdal();
257 
258  if (!G_is_initialized(&st->initialized)) {
259  read_gdal_options();
260  st->projinfo = G_get_projinfo();
261  st->projunits = G_get_projunits();
262  st->projepsg = G_get_projepsg();
263  if (st->projinfo && st->projunits)
264  st->srswkt = GPJ_grass_to_wkt2(st->projinfo, st->projunits,
265  st->projepsg, 0, 0);
266  G_initialize_done(&st->initialized);
267  }
268 
269  gdal = G_calloc(1, sizeof(struct GDAL_link));
270 
271  sprintf(path, "%s/%s%s", st->opts.dir, name, st->opts.ext);
272  gdal->filename = G_store(path);
273  gdal->band_num = 1;
274  gdal->hflip = 0;
275  gdal->vflip = 0;
276 
277  switch (map_type) {
278  case CELL_TYPE:
279  switch (R__.nbytes) {
280  case 1:
281  gdal->type = GDT_Byte;
282  gdal->null_val = (DCELL) 0xFF;
283  break;
284  case 2:
285  gdal->type = GDT_UInt16;
286  gdal->null_val = (DCELL) 0xFFFF;
287  break;
288  case 3:
289  case 4:
290  gdal->type = GDT_Int32;
291  gdal->null_val = (DCELL) 0x80000000U;
292  break;
293  }
294  break;
295  case FCELL_TYPE:
296  gdal->type = GDT_Float32;
297  Rast_set_d_null_value(&gdal->null_val, 1);
298  break;
299  case DCELL_TYPE:
300  gdal->type = GDT_Float64;
301  Rast_set_d_null_value(&gdal->null_val, 1);
302  break;
303  default:
304  G_fatal_error(_("Invalid map type <%d>"), map_type);
305  break;
306  }
307 
308  driver = GDALGetDriverByName(st->opts.format);
309  if (!driver)
310  G_fatal_error(_("Unable to get <%s> driver"), st->opts.format);
311 
312  /* Does driver support GDALCreate ? */
313  if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, NULL)) {
314  gdal->data =
315  GDALCreate(driver, gdal->filename,
317  1, gdal->type, st->opts.options);
318  if (!gdal->data)
319  G_fatal_error(_("Unable to create <%s> dataset using <%s> driver"),
320  name, st->opts.format);
321  }
322  /* If not - create MEM driver for intermediate dataset.
323  * Check if raster can be created at all (with GDALCreateCopy) */
324  else if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATECOPY, NULL)) {
325  GDALDriverH mem_driver;
326 
327  G_message(_("Driver <%s> does not support direct writing. "
328  "Using MEM driver for intermediate dataset."),
329  st->opts.format);
330 
331  mem_driver = GDALGetDriverByName("MEM");
332  if (!mem_driver)
333  G_fatal_error(_("Unable to get in-memory raster driver"));
334 
335  gdal->data =
336  GDALCreate(mem_driver, "",
338  1, gdal->type, st->opts.options);
339  if (!gdal->data)
340  G_fatal_error(_("Unable to create <%s> dataset using memory driver"),
341  name);
342  }
343  else
344  G_fatal_error(_("Driver <%s> does not support creating rasters"),
345  st->opts.format);
346 
347  gdal->band = GDALGetRasterBand(gdal->data, gdal->band_num);
348 
349  GDALSetRasterNoDataValue(gdal->band, gdal->null_val);
350 
351  /* Set Geo Transform */
352  transform[0] = R__.wr_window.west;
353  transform[1] = R__.wr_window.ew_res;
354  transform[2] = 0.0;
355  transform[3] = R__.wr_window.north;
356  transform[4] = 0.0;
357  transform[5] = -R__.wr_window.ns_res;
358 
359  if (GDALSetGeoTransform(gdal->data, transform) >= CE_Failure)
360  G_warning(_("Unable to set geo transform"));
361 
362  if (st->srswkt)
363  if (GDALSetProjection(gdal->data, st->srswkt) == CE_Failure)
364  G_warning(_("Unable to set projection"));
365 
366  fp = G_fopen_new_misc("cell_misc", "gdal", name);
367  if (!fp)
368  G_fatal_error(_("Unable to create cell_misc/%s/gdal file"), name);
369 
370  key_val = G_create_key_value();
371 
372  G_set_key_value("file", gdal->filename, key_val);
373 
374  sprintf(buf, "%d", gdal->band_num);
375  G_set_key_value("band", buf, key_val);
376 
377  sprintf(buf, "%.22g", gdal->null_val);
378  G_set_key_value("null", buf, key_val);
379 
380  sprintf(buf, "%d", gdal->type);
381  G_set_key_value("type", buf, key_val);
382 
383  if (G_fwrite_key_value(fp, key_val) < 0)
384  G_fatal_error(_("Error writing cell_misc/%s/gdal file"), name);
385 
386  G_free_key_value(key_val);
387 
388  fclose(fp);
389 
390  return gdal;
391 #else
392  return NULL;
393 #endif
394 }
395 
396 /*!
397  \brief Close existing GDAL link
398 
399  \param gdal pointer to GDAL_link to be closed
400 */
401 void Rast_close_gdal_link(struct GDAL_link *gdal)
402 {
403 #ifdef GDAL_LINK
404  GDALClose(gdal->data);
405 #endif
406  G_free(gdal->filename);
407  G_free(gdal);
408 }
409 
410 /*!
411  \brief Close existing GDAL link and write out data
412 
413  \param gdal pointer to GDAL_link to be closed
414 
415  \return 1 on success
416  \return -1 on failure
417 */
419 {
420  int stat = 1;
421 
422 #ifdef GDAL_LINK
423  GDALDriverH src_drv = GDALGetDatasetDriver(gdal->data);
424 
425  if (G_strcasecmp(GDALGetDriverShortName(src_drv), "MEM") == 0) {
426  GDALDriverH dst_drv = GDALGetDriverByName(st->opts.format);
427  GDALDatasetH dst =
428  GDALCreateCopy(dst_drv, gdal->filename, gdal->data, FALSE,
429  st->opts.options, NULL, NULL);
430 
431  if (!dst) {
432  G_warning(_("Unable to create output file <%s> using driver <%s>"),
433  gdal->filename, st->opts.format);
434  stat = -1;
435  }
436  GDALClose(dst);
437  }
438 
439  GDALClose(gdal->data);
440 
441 #endif
442  G_free(gdal->filename);
443  G_free(gdal);
444 
445  return stat;
446 }
447 
448 #ifdef GDAL_LINK
449 /*!
450  \brief Input/output function for GDAL links
451 
452  See GDAL's RasterIO for details.
453 */
454 CPLErr Rast_gdal_raster_IO(GDALRasterBandH band, GDALRWFlag rw_flag,
455  int x_off, int y_off, int x_size, int y_size,
456  void *buffer, int buf_x_size, int buf_y_size,
457  GDALDataType buf_type, int pixel_size,
458  int line_size)
459 {
460  return GDALRasterIO(band, rw_flag, x_off, y_off, x_size, y_size,
461  buffer, buf_x_size, buf_y_size, buf_type,
462  pixel_size, line_size);
463 }
464 #endif
#define CELL_TYPE
Definition: raster.h:11
char * G_file_name(char *, const char *, const char *, const char *)
Builds full path names to GIS data files.
Definition: file_name.c:61
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
Definition: get_projinfo.c:61
Definition: R.h:88
CPLErr Rast_gdal_raster_IO(GDALRasterBandH, GDALRWFlag, int, int, int, int, void *, int, int, GDALDataType, int, int)
int Rast_close_gdal_write_link(struct GDAL_link *gdal)
Close existing GDAL link and write out data.
Definition: gdal.c:418
double west
Extent coordinates (west)
Definition: gis.h:475
double DCELL
Definition: gis.h:614
struct Key_Value * G_fread_key_value(FILE *)
Read key/values pairs from file.
Definition: key_value2.c:49
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
char * dst
Definition: lz4.h:599
char ** G_tokenize(const char *, const char *)
Tokenize string.
Definition: gis/token.c:48
#define NULL
Definition: ccmath.h:32
void Rast__init_window(void)
void G_initialize_done(int *)
Definition: counter.c:76
#define G_calloc(m, n)
Definition: defs/gis.h:113
struct Cell_head wr_window
Definition: R.h:100
struct Key_Value * G_get_projunits(void)
Gets units information for location.
Definition: get_projinfo.c:32
void G_message(const char *,...) __attribute__((format(printf
struct state * st
Definition: parser.c:104
#define DCELL_TYPE
Definition: raster.h:13
double north
Extent coordinates (north)
Definition: gis.h:469
struct GDAL_link * Rast_create_gdal_link(const char *name, RASTER_MAP_TYPE map_type)
Create GDAL settings for given raster map.
Definition: gdal.c:242
int int G_strcasecmp(const char *, const char *)
String compare ignoring case (upper or lower)
Definition: strings.c:47
#define FALSE
Definition: gis.h:63
char * GPJ_grass_to_wkt2(const struct Key_Value *, const struct Key_Value *, const struct Key_Value *, int, int)
Converts a GRASS co-ordinate system representation to WKT style. EPSG code is preferred if available...
Definition: convert.c:152
FILE * G_fopen_old_misc(const char *, const char *, const char *, const char *)
open a database misc file for reading
Definition: open_misc.c:210
const struct driver * driver
Definition: driver/init.c:25
void G_set_key_value(const char *, const char *, struct Key_Value *)
Set value for given key.
Definition: key_value1.c:38
#define GPATH_MAX
Definition: gis.h:180
FILE * G_fopen_old(const char *, const char *, const char *)
Open a database file for reading.
Definition: gis/open.c:253
Definition: gis.h:512
int G_fwrite_key_value(FILE *, const struct Key_Value *)
Write key/value pairs to file.
Definition: key_value2.c:25
int cols
Number of columns for 2D data.
Definition: gis.h:442
int G_is_initialized(int *)
Definition: counter.c:59
FILE * G_fopen_new_misc(const char *, const char *, const char *)
open a new database misc file
Definition: open_misc.c:182
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:463
const char * G_mapset(void)
Get current mapset name.
Definition: gis/mapset.c:33
void G_warning(const char *,...) __attribute__((format(printf
void Rast_init_gdal(void)
Initialization.
Definition: gdal.c:39
Definition: path.h:16
#define _(str)
Definition: glocale.h:10
int RASTER_MAP_TYPE
Definition: raster.h:25
#define FCELL_TYPE
Definition: raster.h:12
const char * G_find_raster2(const char *, const char *)
Find a raster map (look but don&#39;t touch)
Definition: find_rast.c:76
char * G_store(const char *)
Copy string to allocated memory.
Definition: strings.c:87
void Rast_close_gdal_link(struct GDAL_link *gdal)
Close existing GDAL link.
Definition: gdal.c:401
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition: key_value1.c:23
int nbytes
Definition: R.h:94
const char * name
Definition: named_colr.c:7
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:459
int rows
Number of rows for 2D data.
Definition: gis.h:438
struct state state
Definition: parser.c:103
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:155
void G_free_key_value(struct Key_Value *)
Free allocated Key_Value structure.
Definition: key_value1.c:103
struct Key_Value * G_get_projepsg(void)
Gets EPSG information for the current location.
Definition: get_projinfo.c:102
int G_make_mapset_object_group(const char *)
Create directory for group of elements of a given type.
Definition: mapset_msc.c:74
struct GDAL_link * Rast_get_gdal_link(const char *name, const char *mapset)
Get GDAL link settings for given raster map.
Definition: gdal.c:61
const char * G_find_key_value(const char *, const struct Key_Value *)
Find given key (case sensitive)
Definition: key_value1.c:84
RASTER_MAP_TYPE Rast_map_type(const char *, const char *)
Determine raster data type.
Definition: raster/open.c:880