GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-498676c526
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 /* GDT_Int8 was introduced in GDAL 3.7 */
129 #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 7, 0)
130  case GDT_Int8:
131 #endif
132  case GDT_Int16:
133  case GDT_UInt16:
134  case GDT_Int32:
135  case GDT_UInt32:
136  req_type = CELL_TYPE;
137  break;
138  case GDT_Float32:
139  req_type = FCELL_TYPE;
140  break;
141  case GDT_Float64:
142  req_type = DCELL_TYPE;
143  break;
144  default:
145  return NULL;
146  }
147 
148  if (req_type != map_type)
149  return NULL;
150 
151  Rast_init_gdal();
152 
153  data = GDALOpen(filename, GA_ReadOnly);
154  if (!data)
155  return NULL;
156 
157  band = GDALGetRasterBand(data, band_num);
158  if (!band) {
159  GDALClose(data);
160  return NULL;
161  }
162 #endif
163 
164  gdal = G_calloc(1, sizeof(struct GDAL_link));
165 
166  gdal->filename = G_store(filename);
167  gdal->band_num = band_num;
168  gdal->null_val = null_val;
169  gdal->hflip = hflip;
170  gdal->vflip = vflip;
171 #ifdef GDAL_LINK
172  gdal->data = data;
173  gdal->band = band;
174  gdal->type = type;
175 #endif
176 
177  return gdal;
178 }
179 
180 struct GDAL_Options {
181  const char *dir;
182  const char *ext;
183  const char *format;
184  char **options;
185 };
186 
187 static struct state {
188  int initialized;
189  struct GDAL_Options opts;
190  struct Key_Value *projinfo, *projunits, *projepsg;
191  char *srswkt;
192 } state;
193 
194 static struct state *st = &state;
195 
196 static void read_gdal_options(void)
197 {
198  FILE *fp;
199  struct Key_Value *key_val;
200  const char *p;
201 
202  fp = G_fopen_old("", "GDAL", G_mapset());
203  if (!fp)
204  G_fatal_error(_("Unable to open GDAL file"));
205  key_val = G_fread_key_value(fp);
206  fclose(fp);
207 
208  p = G_find_key_value("directory", key_val);
209  if (!p)
210  p = "gdal";
211  if (*p == '/') {
212  st->opts.dir = G_store(p);
213  }
214  else {
215  char path[GPATH_MAX];
216 
217  G_file_name(path, p, "", G_mapset());
218  st->opts.dir = G_store(path);
219  if (access(path, 0) != 0)
221  }
222 
223  p = G_find_key_value("extension", key_val);
224  st->opts.ext = G_store(p ? p : "");
225 
226  p = G_find_key_value("format", key_val);
227  st->opts.format = G_store(p ? p : "GTiff");
228 
229  p = G_find_key_value("options", key_val);
230  st->opts.options = p ? G_tokenize(p, ",") : NULL;
231 
232  G_free_key_value(key_val);
233 }
234 
235 /*!
236  \brief Create GDAL settings for given raster map
237 
238  \param name map name
239  \param map_type map type (CELL, FCELL, DCELL)
240 
241  \return pointer to allocated GDAL_link structure
242  \return NULL on error
243  */
245  RASTER_MAP_TYPE map_type)
246 {
247 #ifdef GDAL_LINK
248  char path[GPATH_MAX];
249  GDALDriverH driver;
250  double transform[6];
251  struct GDAL_link *gdal;
252  FILE *fp;
253  struct Key_Value *key_val;
254  char buf[32];
255 
257 
258  Rast_init_gdal();
259 
260  if (!G_is_initialized(&st->initialized)) {
261  read_gdal_options();
262  st->projinfo = G_get_projinfo();
263  st->projunits = G_get_projunits();
264  st->projepsg = G_get_projepsg();
265  if (st->projinfo && st->projunits)
266  st->srswkt = GPJ_grass_to_wkt2(st->projinfo, st->projunits,
267  st->projepsg, 0, 0);
268  G_initialize_done(&st->initialized);
269  }
270 
271  gdal = G_calloc(1, sizeof(struct GDAL_link));
272 
273  sprintf(path, "%s/%s%s", st->opts.dir, name, st->opts.ext);
274  gdal->filename = G_store(path);
275  gdal->band_num = 1;
276  gdal->hflip = 0;
277  gdal->vflip = 0;
278 
279  switch (map_type) {
280  case CELL_TYPE:
281  switch (R__.nbytes) {
282  case 1:
283  gdal->type = GDT_Byte;
284  gdal->null_val = (DCELL)0xFF;
285  break;
286  case 2:
287  gdal->type = GDT_UInt16;
288  gdal->null_val = (DCELL)0xFFFF;
289  break;
290  case 3:
291  case 4:
292  gdal->type = GDT_Int32;
293  gdal->null_val = (DCELL)0x80000000U;
294  break;
295  }
296  break;
297  case FCELL_TYPE:
298  gdal->type = GDT_Float32;
299  Rast_set_d_null_value(&gdal->null_val, 1);
300  break;
301  case DCELL_TYPE:
302  gdal->type = GDT_Float64;
303  Rast_set_d_null_value(&gdal->null_val, 1);
304  break;
305  default:
306  G_fatal_error(_("Invalid map type <%d>"), map_type);
307  break;
308  }
309 
310  driver = GDALGetDriverByName(st->opts.format);
311  if (!driver)
312  G_fatal_error(_("Unable to get <%s> driver"), st->opts.format);
313 
314  /* Does driver support GDALCreate ? */
315  if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, NULL)) {
316  gdal->data =
317  GDALCreate(driver, gdal->filename, R__.wr_window.cols,
318  R__.wr_window.rows, 1, gdal->type, st->opts.options);
319  if (!gdal->data)
320  G_fatal_error(_("Unable to create <%s> dataset using <%s> driver"),
321  name, st->opts.format);
322  }
323  /* If not - create MEM driver for intermediate dataset.
324  * Check if raster can be created at all (with GDALCreateCopy) */
325  else if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATECOPY, NULL)) {
326  GDALDriverH mem_driver;
327 
328  G_message(_("Driver <%s> does not support direct writing. "
329  "Using MEM driver for intermediate dataset."),
330  st->opts.format);
331 
332  mem_driver = GDALGetDriverByName("MEM");
333  if (!mem_driver)
334  G_fatal_error(_("Unable to get in-memory raster driver"));
335 
336  gdal->data =
337  GDALCreate(mem_driver, "", R__.wr_window.cols, R__.wr_window.rows,
338  1, gdal->type, st->opts.options);
339  if (!gdal->data)
341  _("Unable to create <%s> dataset using memory driver"), 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 */
354  transform[2] = 0.0;
356  transform[4] = 0.0;
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 = GDALCreateCopy(dst_drv, gdal->filename, gdal->data,
428  FALSE, st->opts.options, NULL, NULL);
429 
430  if (!dst) {
431  G_warning(_("Unable to create output file <%s> using driver <%s>"),
432  gdal->filename, st->opts.format);
433  stat = -1;
434  }
435  GDALClose(dst);
436  }
437 
438  GDALClose(gdal->data);
439 
440 #endif
441  G_free(gdal->filename);
442  G_free(gdal);
443 
444  return stat;
445 }
446 
447 #ifdef GDAL_LINK
448 /*!
449  \brief Input/output function for GDAL links
450 
451  See GDAL's RasterIO for details.
452  */
453 CPLErr Rast_gdal_raster_IO(GDALRasterBandH band, GDALRWFlag rw_flag, int x_off,
454  int y_off, int x_size, int y_size, void *buffer,
455  int buf_x_size, int buf_y_size,
456  GDALDataType buf_type, int pixel_size, int line_size)
457 {
458  return GDALRasterIO(band, rw_flag, x_off, y_off, x_size, y_size, buffer,
459  buf_x_size, buf_y_size, buf_type, pixel_size,
460  line_size);
461 }
462 #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:418
void Rast_close_gdal_link(struct GDAL_link *gdal)
Close existing GDAL link.
Definition: gdal.c:401
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:244
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:628
#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:475
double north
Extent coordinates (north)
Definition: gis.h:485
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:479
int rows
Number of rows for 2D data.
Definition: gis.h:454
int cols
Number of columns for 2D data.
Definition: gis.h:458
double west
Extent coordinates (west)
Definition: gis.h:491
Definition: gis.h:527
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