GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-535c39c9fc
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 || G_strcasecmp(p, "nan") == 0 ||
111  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  const char *dir;
178  const char *ext;
179  const char *format;
180  char **options;
181 };
182 
183 static struct state {
184  int initialized;
185  struct GDAL_Options opts;
186  struct Key_Value *projinfo, *projunits, *projepsg;
187  char *srswkt;
188 } state;
189 
190 static struct state *st = &state;
191 
192 static void read_gdal_options(void)
193 {
194  FILE *fp;
195  struct Key_Value *key_val;
196  const char *p;
197 
198  fp = G_fopen_old("", "GDAL", G_mapset());
199  if (!fp)
200  G_fatal_error(_("Unable to open GDAL file"));
201  key_val = G_fread_key_value(fp);
202  fclose(fp);
203 
204  p = G_find_key_value("directory", key_val);
205  if (!p)
206  p = "gdal";
207  if (*p == '/') {
208  st->opts.dir = G_store(p);
209  }
210  else {
211  char path[GPATH_MAX];
212 
213  G_file_name(path, p, "", G_mapset());
214  st->opts.dir = G_store(path);
215  if (access(path, 0) != 0)
217  }
218 
219  p = G_find_key_value("extension", key_val);
220  st->opts.ext = G_store(p ? p : "");
221 
222  p = G_find_key_value("format", key_val);
223  st->opts.format = G_store(p ? p : "GTiff");
224 
225  p = G_find_key_value("options", key_val);
226  st->opts.options = p ? G_tokenize(p, ",") : NULL;
227 
228  G_free_key_value(key_val);
229 }
230 
231 /*!
232  \brief Create GDAL settings for given raster map
233 
234  \param name map name
235  \param map_type map type (CELL, FCELL, DCELL)
236 
237  \return pointer to allocated GDAL_link structure
238  \return NULL on error
239  */
241  RASTER_MAP_TYPE map_type)
242 {
243 #ifdef GDAL_LINK
244  char path[GPATH_MAX];
245  GDALDriverH driver;
246  double transform[6];
247  struct GDAL_link *gdal;
248  FILE *fp;
249  struct Key_Value *key_val;
250  char buf[32];
251 
253 
254  Rast_init_gdal();
255 
256  if (!G_is_initialized(&st->initialized)) {
257  read_gdal_options();
258  st->projinfo = G_get_projinfo();
259  st->projunits = G_get_projunits();
260  st->projepsg = G_get_projepsg();
261  if (st->projinfo && st->projunits)
262  st->srswkt = GPJ_grass_to_wkt2(st->projinfo, st->projunits,
263  st->projepsg, 0, 0);
264  G_initialize_done(&st->initialized);
265  }
266 
267  gdal = G_calloc(1, sizeof(struct GDAL_link));
268 
269  sprintf(path, "%s/%s%s", st->opts.dir, name, st->opts.ext);
270  gdal->filename = G_store(path);
271  gdal->band_num = 1;
272  gdal->hflip = 0;
273  gdal->vflip = 0;
274 
275  switch (map_type) {
276  case CELL_TYPE:
277  switch (R__.nbytes) {
278  case 1:
279  gdal->type = GDT_Byte;
280  gdal->null_val = (DCELL)0xFF;
281  break;
282  case 2:
283  gdal->type = GDT_UInt16;
284  gdal->null_val = (DCELL)0xFFFF;
285  break;
286  case 3:
287  case 4:
288  gdal->type = GDT_Int32;
289  gdal->null_val = (DCELL)0x80000000U;
290  break;
291  }
292  break;
293  case FCELL_TYPE:
294  gdal->type = GDT_Float32;
295  Rast_set_d_null_value(&gdal->null_val, 1);
296  break;
297  case DCELL_TYPE:
298  gdal->type = GDT_Float64;
299  Rast_set_d_null_value(&gdal->null_val, 1);
300  break;
301  default:
302  G_fatal_error(_("Invalid map type <%d>"), map_type);
303  break;
304  }
305 
306  driver = GDALGetDriverByName(st->opts.format);
307  if (!driver)
308  G_fatal_error(_("Unable to get <%s> driver"), st->opts.format);
309 
310  /* Does driver support GDALCreate ? */
311  if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, NULL)) {
312  gdal->data =
313  GDALCreate(driver, gdal->filename, R__.wr_window.cols,
314  R__.wr_window.rows, 1, gdal->type, st->opts.options);
315  if (!gdal->data)
316  G_fatal_error(_("Unable to create <%s> dataset using <%s> driver"),
317  name, st->opts.format);
318  }
319  /* If not - create MEM driver for intermediate dataset.
320  * Check if raster can be created at all (with GDALCreateCopy) */
321  else if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATECOPY, NULL)) {
322  GDALDriverH mem_driver;
323 
324  G_message(_("Driver <%s> does not support direct writing. "
325  "Using MEM driver for intermediate dataset."),
326  st->opts.format);
327 
328  mem_driver = GDALGetDriverByName("MEM");
329  if (!mem_driver)
330  G_fatal_error(_("Unable to get in-memory raster driver"));
331 
332  gdal->data =
333  GDALCreate(mem_driver, "", R__.wr_window.cols, R__.wr_window.rows,
334  1, gdal->type, st->opts.options);
335  if (!gdal->data)
337  _("Unable to create <%s> dataset using memory driver"), name);
338  }
339  else
340  G_fatal_error(_("Driver <%s> does not support creating rasters"),
341  st->opts.format);
342 
343  gdal->band = GDALGetRasterBand(gdal->data, gdal->band_num);
344 
345  GDALSetRasterNoDataValue(gdal->band, gdal->null_val);
346 
347  /* Set Geo Transform */
350  transform[2] = 0.0;
352  transform[4] = 0.0;
354 
355  if (GDALSetGeoTransform(gdal->data, transform) >= CE_Failure)
356  G_warning(_("Unable to set geo transform"));
357 
358  if (st->srswkt)
359  if (GDALSetProjection(gdal->data, st->srswkt) == CE_Failure)
360  G_warning(_("Unable to set projection"));
361 
362  fp = G_fopen_new_misc("cell_misc", "gdal", name);
363  if (!fp)
364  G_fatal_error(_("Unable to create cell_misc/%s/gdal file"), name);
365 
366  key_val = G_create_key_value();
367 
368  G_set_key_value("file", gdal->filename, key_val);
369 
370  sprintf(buf, "%d", gdal->band_num);
371  G_set_key_value("band", buf, key_val);
372 
373  sprintf(buf, "%.22g", gdal->null_val);
374  G_set_key_value("null", buf, key_val);
375 
376  sprintf(buf, "%d", gdal->type);
377  G_set_key_value("type", buf, key_val);
378 
379  if (G_fwrite_key_value(fp, key_val) < 0)
380  G_fatal_error(_("Error writing cell_misc/%s/gdal file"), name);
381 
382  G_free_key_value(key_val);
383 
384  fclose(fp);
385 
386  return gdal;
387 #else
388  return NULL;
389 #endif
390 }
391 
392 /*!
393  \brief Close existing GDAL link
394 
395  \param gdal pointer to GDAL_link to be closed
396  */
397 void Rast_close_gdal_link(struct GDAL_link *gdal)
398 {
399 #ifdef GDAL_LINK
400  GDALClose(gdal->data);
401 #endif
402  G_free(gdal->filename);
403  G_free(gdal);
404 }
405 
406 /*!
407  \brief Close existing GDAL link and write out data
408 
409  \param gdal pointer to GDAL_link to be closed
410 
411  \return 1 on success
412  \return -1 on failure
413  */
415 {
416  int stat = 1;
417 
418 #ifdef GDAL_LINK
419  GDALDriverH src_drv = GDALGetDatasetDriver(gdal->data);
420 
421  if (G_strcasecmp(GDALGetDriverShortName(src_drv), "MEM") == 0) {
422  GDALDriverH dst_drv = GDALGetDriverByName(st->opts.format);
423  GDALDatasetH dst = GDALCreateCopy(dst_drv, gdal->filename, gdal->data,
424  FALSE, st->opts.options, NULL, NULL);
425 
426  if (!dst) {
427  G_warning(_("Unable to create output file <%s> using driver <%s>"),
428  gdal->filename, st->opts.format);
429  stat = -1;
430  }
431  GDALClose(dst);
432  }
433 
434  GDALClose(gdal->data);
435 
436 #endif
437  G_free(gdal->filename);
438  G_free(gdal);
439 
440  return stat;
441 }
442 
443 #ifdef GDAL_LINK
444 /*!
445  \brief Input/output function for GDAL links
446 
447  See GDAL's RasterIO for details.
448  */
449 CPLErr Rast_gdal_raster_IO(GDALRasterBandH band, GDALRWFlag rw_flag, int x_off,
450  int y_off, int x_size, int y_size, void *buffer,
451  int buf_x_size, int buf_y_size,
452  GDALDataType buf_type, int pixel_size, int line_size)
453 {
454  return GDALRasterIO(band, rw_flag, x_off, y_off, x_size, y_size, buffer,
455  buf_x_size, buf_y_size, buf_type, pixel_size,
456  line_size);
457 }
458 #endif
CPLErr Rast_gdal_raster_IO(GDALRasterBandH, GDALRWFlag, int, int, int, int, void *, int, int, GDALDataType, int, int)
#define NULL
Definition: ccmath.h:32
FILE * G_fopen_old(const char *, const char *, const char *)
Open a database file for reading.
Definition: gis/open.c:251
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
Definition: get_projinfo.c:61
FILE * G_fopen_old_misc(const char *, const char *, const char *, const char *)
open a database misc file for reading
Definition: open_misc.c:205
#define G_calloc(m, n)
Definition: defs/gis.h:95
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
struct Key_Value * G_get_projunits(void)
Gets units information for location.
Definition: get_projinfo.c:32
void G_warning(const char *,...) __attribute__((format(printf
char * G_file_name(char *, const char *, const char *, const char *)
Builds full path names to GIS data files.
Definition: file_name.c:61
const char * G_find_key_value(const char *, const struct Key_Value *)
Find given key (case sensitive)
Definition: key_value1.c:85
int G_fwrite_key_value(FILE *, const struct Key_Value *)
Write key/value pairs to file.
Definition: key_value2.c:25
int G_make_mapset_object_group(const char *)
Create directory for group of elements of a given type.
Definition: mapset_msc.c:74
const char * G_find_raster2(const char *, const char *)
Find a raster map (look but don't touch)
Definition: find_rast.c:76
void G_free_key_value(struct Key_Value *)
Free allocated Key_Value structure.
Definition: key_value1.c:104
FILE * G_fopen_new_misc(const char *, const char *, const char *)
open a new database misc file
Definition: open_misc.c:178
const char * G_mapset(void)
Get current mapset name.
Definition: gis/mapset.c:33
void G_set_key_value(const char *, const char *, struct Key_Value *)
Set value for given key.
Definition: key_value1.c:39
struct Key_Value * G_fread_key_value(FILE *)
Read key/values pairs from file.
Definition: key_value2.c:49
int G_is_initialized(int *)
Definition: counter.c:60
void G_initialize_done(int *)
Definition: counter.c:77
int int G_strcasecmp(const char *, const char *)
String compare ignoring case (upper or lower)
Definition: strings.c:47
struct Key_Value * G_get_projepsg(void)
Gets EPSG information for the current location.
Definition: get_projinfo.c:102
void G_message(const char *,...) __attribute__((format(printf
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition: key_value1.c:23
char * G_store(const char *)
Copy string to allocated memory.
Definition: strings.c:87
char ** G_tokenize(const char *, const char *)
Tokenize string.
Definition: gis/token.c:47
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:150
void Rast__init_window(void)
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:153
RASTER_MAP_TYPE Rast_map_type(const char *, const char *)
Determine raster data type.
Definition: raster/open.c:894
const struct driver * driver
Definition: driver/init.c:25
int Rast_close_gdal_write_link(struct GDAL_link *gdal)
Close existing GDAL link and write out data.
Definition: gdal.c:414
void Rast_close_gdal_link(struct GDAL_link *gdal)
Close existing GDAL link.
Definition: gdal.c:397
void Rast_init_gdal(void)
Initialization.
Definition: gdal.c:39
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:240
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
#define GPATH_MAX
Definition: gis.h:194
#define FALSE
Definition: gis.h:83
double DCELL
Definition: gis.h:626
#define _(str)
Definition: glocale.h:10
const char * name
Definition: named_colr.c:6
struct state state
Definition: parser.c:103
struct state * st
Definition: parser.c:104
#define FCELL_TYPE
Definition: raster.h:12
#define DCELL_TYPE
Definition: raster.h:13
#define CELL_TYPE
Definition: raster.h:11
int RASTER_MAP_TYPE
Definition: raster.h:25
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:473
double north
Extent coordinates (north)
Definition: gis.h:483
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:477
int rows
Number of rows for 2D data.
Definition: gis.h:452
int cols
Number of columns for 2D data.
Definition: gis.h:456
double west
Extent coordinates (west)
Definition: gis.h:489
Definition: gis.h:525
Definition: R.h:87
int nbytes
Definition: R.h:92
struct Cell_head wr_window
Definition: R.h:98
Definition: driver.h:21
Definition: path.h:15