Package 'gdalraster'

Description

gdalraster is an interface to the Geospatial Data Abstraction Library (GDAL) for low level raster I/O. Calling signatures resemble those of the native C, C++ and Python APIs provided by the GDAL project. See https://gdal.org/api/ for details of the GDAL Raster API.

Details

Core functionality is contained in class GDALRaster and several related stand-alone functions:

• GDALRaster-class is an exposed C++ class that allows opening a raster dataset and calling methods on the GDALDataset, GDALDriver and GDALRasterBand objects in the underlying API (e.g., get/set parameters, read/write pixel data).

• raster creation: create(), createCopy(), rasterFromRaster(), translate(), getCreationOptions()

• virtual raster: buildVRT(), rasterToVRT()

• reproject/resample/crop/mosaic: warp()

• algorithms: dem_proc(), fillNodata(), footprint(), polygonize(), rasterize(), sieveFilter(), GDALRaster$getChecksum()

• raster attribute tables: buildRAT(), displayRAT(), GDALRaster$getDefaultRAT(), GDALRaster$setDefaultRAT()

• geotransform conversion: apply_geotransform(), get_pixel_line(), inv_geotransform()

• coordinate transformation: transform_xy(), inv_project()

• spatial reference convenience functions: epsg_to_wkt(), srs_to_wkt(),
  srs_is_geographic(), srs_is_projected(), srs_is_same()

• geometry convenience functions: bbox_from_wkt(), bbox_to_wkt(), bbox_intersect(), bbox_union(), bbox_transform(), g_area(), g_buffer(), g_centroid(), g_contains(), g_crosses(), g_difference(), g_disjoint(), g_distance(), g_equals(), g_intersection(), g_intersects(), g_is_empty(), g_is_valid(), g_length(), g_name(), g_overlaps(), g_sym_difference(), g_touches(), g_transform(), g_union(), g_within(), geos_version()

• data management: addFilesInZip(), copyDatasetFiles(), deleteDataset(),
  renameDataset(), bandCopyWholeRaster()

• OGR vector utilities: ogr2ogr(), ogrinfo(), ogr_manage, ogr_define

• virtual file systems: VSIFile, vsi_clear_path_options(), vsi_copy_file(), vsi_curl_clear_cache(), vsi_get_disk_free_space(), vsi_get_file_metadata(), vsi_get_fs_options(), vsi_get_fs_prefixes(), vsi_mkdir(), vsi_read_dir(), vsi_rename(), vsi_rmdir(), vsi_set_path_option(), vsi_stat(), vsi_supports_rnd_write(), vsi_supports_seq_write(), vsi_sync(), vsi_unlink(), vsi_unlink_batch()

• GDAL configuration: gdal_version(), gdal_formats(), get_cache_used(), get_config_option(), set_config_option(), get_num_cpus(), get_usable_physical_ram(), has_spatialite(), http_enabled(), push_error_handler(), pop_error_handler(), dump_open_datasets()

• PROJ configuration: proj_version(), proj_search_paths(), proj_networking()

Additional functionality includes:

• RunningStats-class calculates mean and variance in one pass. The min, max, sum, and count are also tracked (efficient summary statistics on data streams).

• CmbTable-class implements a hash table for counting unique combinations of integer values.

• combine() overlays multiple rasters so that a unique ID is assigned to each unique combination of input values. Pixel counts for each unique combination are obtained, and combination IDs are optionally written to an output raster.

• calc() evaluates an R expression for each pixel in a raster layer or stack of layers. Individual pixel coordinates are available as variables in the R expression, as either x/y in the raster projected coordinate system or inverse projected longitude/latitude.

• plot_raster() displays raster data using base R graphics. Supports single-band grayscale, RGB, color tables and color map functions (e.g., color ramp).

Note

Documentation for GDALRaster-class and several wrapper functions borrows from the GDAL API documentation, (c) 1998-2024, Frank Warmerdam, Even Rouault, and others, MIT license.

Sample datasets included with the package are used in examples throughout the documentation. The sample data include LANDFIRE raster layers describing terrain, vegetation and wildland fuels (LF 2020 version), Landsat C2 Analysis Ready Data from USGS Earth Explorer, and Monitoring Trends in Burn Severity (MTBS) fire perimeters from 1984-2022. Metadata for the sample datasets are in inst/extdata/metadata.zip.

system.file() is used in the examples to access the sample datasets. This enables the code to run regardless of where R is installed. Users will normally give file names as a regular full path or relative to the current working directory.

Author(s)

GDAL is by: Frank Warmerdam, Even Rouault and others
(see https://github.com/OSGeo/gdal/graphs/contributors)

R interface/additional functionality: Chris Toney

Maintainer: Chris Toney <chris.toney at usda.gov>

See Also

GDAL Raster Data Model:
https://gdal.org/user/raster_data_model.html

Raster format descriptions:
https://gdal.org/drivers/raster/index.html

Geotransform tutorial:
https://gdal.org/tutorials/geotransforms_tut.html

GDAL Virtual File Systems:
https://gdal.org/user/virtual_file_systems.html

Description

addFilesInZip() will create new or open existing ZIP file, and add one or more compressed files potentially using the seek optimization extension. This function is basically a wrapper for CPLAddFileInZip() in the GDAL Common Portability Library, but optionally creates a new ZIP file first (with CPLCreateZip()). It provides a subset of functionality in the GDAL sozip command-line utility (https://gdal.org/programs/sozip.html). Requires GDAL >= 3.7.

Usage

addFilesInZip(
  zip_file,
  add_files,
  overwrite = FALSE,
  full_paths = TRUE,
  sozip_enabled = NULL,
  sozip_chunk_size = NULL,
  sozip_min_file_size = NULL,
  num_threads = NULL,
  content_type = NULL,
  quiet = FALSE
)

Arguments

Details

A Seek-Optimized ZIP file (SOZip) contains one or more compressed files organized and annotated such that a SOZip-aware reader can perform very fast random access within the .zip file (see https://github.com/sozip/sozip-spec). Large compressed files can be accessed directly from SOZip without prior decompression. The .zip file is otherwise fully backward compatible.

If sozip_enabled="AUTO" (the default), a file is seek-optimized only if its size is above the values of sozip_min_file_size (default 1 MB) and sozip_chunk_size (default 32768). In "YES" mode, all input files will be seek-optimized. In "NO" mode, no input files will be seek-optimized. The default can be changed with the CPL_SOZIP_ENABLED configuration option.

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

Note

The GDAL_NUM_THREADS configuration option can be set to ALL_CPUS or an integer value to specify the number of threads to use for SOZip-compressed files (see set_config_option()).

SOZip can be validated with:

vsi_get_file_metadata(zip_file, domain="ZIP")

where zip_file uses the /vsizip/ prefix.

See Also

vsi_get_file_metadata()

Examples

lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
zip_file <- file.path(tempdir(), "storml_lcp.zip")

# Requires GDAL >= 3.7
if (as.integer(gdal_version()[2]) >= 3070000) {
  addFilesInZip(zip_file, lcp_file, full_paths=FALSE, sozip_enabled="YES",
                num_threads=1)

  print("Files in zip archive:")
  print(unzip(zip_file, list=TRUE))

  # Open with GDAL using Virtual File System handler '/vsizip/'
  # see: https://gdal.org/user/virtual_file_systems.html#vsizip-zip-archives
  lcp_in_zip <- file.path("/vsizip", zip_file, "storm_lake.lcp")
  print("SOZip metadata:")
  print(vsi_get_file_metadata(lcp_in_zip, domain="ZIP"))

  ds <- new(GDALRaster, lcp_in_zip)
  ds$info()
  ds$close()

  vsi_unlink(zip_file)
}

Description

apply_geotransform() applies geotransform coefficients to raster coordinates in pixel/line space (column/row), converting into georeferenced (x/y) coordinates. Wrapper of GDALApplyGeoTransform() in the GDAL API, operating on matrix input.

Usage

apply_geotransform(col_row, gt)

Arguments

Value

Numeric matrix of geospatial x/y coordinates.

Note

Bounds checking on the input coordinates is done if gt is obtained from an object of class GDALRaster. See Note for get_pixel_line().

See Also

GDALRaster$getGeoTransform(), get_pixel_line()

Examples

raster_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds <- new(GDALRaster, raster_file)

# compute some raster coordinates in column/row space
set.seed(42)
col_coords <- runif(10, min = 0, max = ds$getRasterXSize() - 0.00001)
row_coords <- runif(10, min = 0, max = ds$getRasterYSize() - 0.00001)
col_row <- cbind(col_coords, row_coords)

# convert to geospatial x/y coordinates
gt <- ds$getGeoTransform()
apply_geotransform(col_row, gt)

# or, using the class method
ds$apply_geotransform(col_row)

# bounds checking
col_row <- rbind(col_row, c(ds$getRasterXSize(), ds$getRasterYSize()))
ds$apply_geotransform(col_row)

ds$close()

Description

bandCopyWholeRaster() copies the complete raster contents of one band to another similarly configured band. The source and destination bands must have the same xsize and ysize. The bands do not have to have the same data type. It implements efficient copying, in particular "chunking" the copy in substantial blocks. This is a wrapper for GDALRasterBandCopyWholeRaster() in the GDAL API.

Usage

bandCopyWholeRaster(
  src_filename,
  src_band,
  dst_filename,
  dst_band,
  options = NULL,
  quiet = FALSE
)

Arguments

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

See Also

GDALRaster-class, create(), createCopy(), rasterFromRaster()

Examples

## copy Landsat data from a single-band file to a new multi-band image
b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
dst_file <- file.path(tempdir(), "sr_multi.tif")
rasterFromRaster(b5_file, dst_file, nbands=7, init=0)
opt <- c("COMPRESSED=YES", "SKIP_HOLES=YES")
bandCopyWholeRaster(b5_file, 1, dst_file, 5, options=opt)
ds <- new(GDALRaster, dst_file)
ds$getStatistics(band=5, approx_ok=FALSE, force=TRUE)
ds$close()

deleteDataset(dst_file)

Description

bbox_from_wkt() returns the bounding box of a WKT 2D geometry (e.g., LINE, POLYGON, MULTIPOLYGON).

Usage

bbox_from_wkt(wkt, extend_x = 0, extend_y = 0)

Arguments

Value

Numeric vector of length four containing the xmin, ymin, xmax, ymax of the geometry specified by wkt (possibly extended by values in extend_x, extend_y).

See Also

bbox_to_wkt()

Examples

bnd <- "POLYGON ((324467.3 5104814.2, 323909.4 5104365.4, 323794.2
5103455.8, 324970.7 5102885.8, 326420.0 5103595.3, 326389.6 5104747.5,
325298.1 5104929.4, 325298.1 5104929.4, 324467.3 5104814.2))"
bbox_from_wkt(bnd, 100, 100)

Description

bbox_intersect() returns the bounding box intersection, and bbox_union() returns the bounding box union, for input of either raster file names or list of bounding boxes. All of the inputs must be in the same projected coordinate system.

Usage

bbox_intersect(x, as_wkt = FALSE)

bbox_union(x, as_wkt = FALSE)

Arguments

Value

The intersection (bbox_intersect()) or union (bbox_union()) of inputs. If as_wkt = FALSE, a numeric vector of length four containing xmin, ymin, xmax, ymax. If as_wkt = TRUE, a character string containing OGC WKT for the bbox as POLYGON.

See Also

bbox_from_wkt(), bbox_to_wkt()

Examples

bbox_list <-list()

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)
bbox_list[[1]] <- ds$bbox()
ds$close()

b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
ds <- new(GDALRaster, b5_file)
bbox_list[[2]] <- ds$bbox()
ds$close()

bnd <- "POLYGON ((324467.3 5104814.2, 323909.4 5104365.4, 323794.2
5103455.8, 324970.7 5102885.8, 326420.0 5103595.3, 326389.6 5104747.5,
325298.1 5104929.4, 325298.1 5104929.4, 324467.3 5104814.2))"
bbox_list[[3]] <- bbox_from_wkt(bnd)

print(bbox_list)
bbox_intersect(bbox_list)
bbox_union(bbox_list)

Description

bbox_to_wkt() returns a WKT POLYGON string for the given bounding box. Requires GDAL built with the GEOS library.

Usage

bbox_to_wkt(bbox, extend_x = 0, extend_y = 0)

Arguments

Value

Character string for an OGC WKT polygon. NA is returned if GDAL was built without the GEOS library.

See Also

bbox_from_wkt(), g_buffer()

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file, read_only=TRUE)
bbox_to_wkt(ds$bbox())
ds$close()

Description

bbox_transform() is a convenience function for:

bbox_to_wkt(bbox) |>
  g_transform(srs_from, srs_to) |>
  bbox_from_wkt()

Usage

bbox_transform(bbox, srs_from, srs_to)

Arguments

Value

Numeric vector of length four containing a transformed bounding box (xmin, ymin, xmax, ymax).

See Also

g_transform(), bbox_from_wkt(), bbox_to_wkt()

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)
ds$bbox()
bbox_transform(ds$bbox(), ds$getProjection(), epsg_to_wkt(4326))
ds$close()

Description

buildRAT() reads all pixels of an input raster to obtain the set of unique values and their counts. The result is returned as a data frame suitable for use with the class method GDALRaster$setDefaultRAT(). The returned data frame might be further modified before setting as a Raster Attribute Table in a dataset, for example, by adding columns containing class names, color values, or other information (see Details). An optional input data frame containing such attributes may be given, in which case buildRAT() will attempt to join the additional columns and automatically assign the appropriate metadata on the output data frame (i.e., assign R attributes on the data frame and its columns that define usage in a GDAL Raster Attribute Table).

Usage

buildRAT(
  raster,
  band = 1L,
  col_names = c("VALUE", "COUNT"),
  table_type = "athematic",
  na_value = NULL,
  join_df = NULL,
  quiet = FALSE
)

Arguments

Details

A GDAL Raster Attribute Table (or RAT) provides attribute information about pixel values. Raster attribute tables can be used to represent histograms, color tables, and classification information. Each row in the table applies to either a single pixel value or a range of values, and might have attributes such as the histogram count for that value (or range), the color that pixels of that value (or range) should be displayed, names of classes, or various other information.

Each column in a raster attribute table has a name, a type (integer, double, or string), and a GDALRATFieldUsage. The usage distinguishes columns with particular understood purposes (such as color, histogram count, class name), and columns that have other purposes not understood by the library (long labels, ancillary attributes, etc).

In the general case, each row has a field indicating the minimum pixel value falling into that category, and a field indicating the maximum pixel value. In the GDAL API, these are indicated with usage values of GFU_Min and GFU_Max. In the common case where each row is a discrete pixel value, a single column with usage GFU_MinMax would be used instead. In R, the table is represented as a data frame with column attribute "GFU" containing the field usage as a string, e.g., "Max", "Min" or "MinMax" (case-sensitive). The full set of possible field usage descriptors is:

buildRAT() assigns GFU "MinMax" on the column of pixel values (named "VALUE" by default) and GFU "PixelCount" on the column of counts (named "COUNT" by default). If join_df is given, the additional columns that result from joining will have GFU assigned automatically based on the column names (ignoring case). First, the additional column names are checked for containing the string "name" (e.g., "classname", "TypeName", "EVT_NAME", etc). The first matching column (if any) will be assigned a GFU of "Name" (=GFU_Name, the field usage descriptor for class names). Next, columns named "R" or "Red" will be assigned GFU "Red", columns named "G" or "Green" will be assigned GFU "Green", columns named "B" or "Blue" will be assigned GFU "Blue", and columns named "A" or "Alpha" will be assigned GFU "Alpha". Finally, any remaining columns that have not been assigned a GFU will be assigned "Generic".

In a variation of RAT, all the categories are of equal size and regularly spaced, and the categorization can be determined by knowing the value at which the categories start and the size of a category. This is called "Linear Binning" and the information is kept specially on the raster attribute table as a whole. In R, a RAT that uses linear binning would have the following attributes set on the data frame: attribute "Row0Min" = the numeric lower bound (pixel value) of the first category, and attribute "BinSize" = the numeric width of each category (in pixel value units). buildRAT() does not create tables with linear binning, but one could be created manually based on the specifications above, and applied to a raster with the class method GDALRaster$setDefaultRAT().

A raster attribute table is thematic or athematic (continuous). In R, this is defined by an attribute on the data frame named "GDALRATTableType" with value of either "thematic" or "athematic".

Value

A data frame with at least two columns containing the set of unique pixel values and their counts. These columns have attribute "GFU" set to "MinMax" for the values, and "PixelCount" for the counts. If join_df is given, the returned data frame will have additional columns that result from merge(). The "GFU" attribute of the additional columns will be assigned automatically based on the column names (case-insensitive matching, see Details). The returned data frame has attribute "GDALRATTableType" set to table_type.

Note

The full raster will be scanned.

If na_value is not specified, then an NA pixel value (if present) will not be recoded in the output data frame. This may have implications if joining to other data (NA will not match), or when using the returned data frame to set a default RAT on a dataset (NA will be interpreted as the value that R uses internally to represent it for the type, e.g., -2147483648 for NA_integer_). In some cases, removing the row in the output data frame with value NA, rather than recoding, may be desirable (i.e., by removing manually or by side effect of joining via merge(), for example). Users should consider what is appropriate for a particular case.

See Also

GDALRaster$getDefaultRAT(), GDALRaster$setDefaultRAT(), displayRAT()

vignette("raster-attribute-tables")

Examples

evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
# make a copy to modify
f <- file.path(tempdir(), "storml_evt_tmp.tif")
file.copy(evt_file,  f)

ds <- new(GDALRaster, f, read_only=FALSE)
ds$getDefaultRAT(band=1) # NULL

# get the full attribute table for LANDFIRE EVT from the CSV file
evt_csv <- system.file("extdata/LF20_EVT_220.csv", package="gdalraster")
evt_df <- read.csv(evt_csv)
nrow(evt_df)
head(evt_df)
evt_df <- evt_df[,1:7]

tbl <- buildRAT(ds,
                table_type = "thematic",
                na_value = -9999,
                join_df = evt_df)

nrow(tbl)
head(tbl)

# attributes on the data frame and its columns define usage in a GDAL RAT
attributes(tbl)
attributes(tbl$VALUE)
attributes(tbl$COUNT)
attributes(tbl$EVT_NAME)
attributes(tbl$EVT_LF)
attributes(tbl$EVT_PHYS)
attributes(tbl$R)
attributes(tbl$G)
attributes(tbl$B)

ds$setDefaultRAT(band=1, tbl)
ds$flushCache()

tbl2 <- ds$getDefaultRAT(band=1)
nrow(tbl2)
head(tbl2)

ds$close()
deleteDataset(f)

# Display
evt_gt <- displayRAT(tbl2, title = "Storm Lake EVT Raster Attribute Table")
class(evt_gt)  # an object of class "gt_tbl" from package gt
# To show the table:
# evt_gt
# or simply call `displayRAT()` as above but without assignment
# `vignette("raster-attribute-tables")` has example output

Description

buildVRT() is a wrapper of the gdalbuildvrt command-line utility for building a VRT (Virtual Dataset) that is a mosaic of the list of input GDAL datasets (see https://gdal.org/programs/gdalbuildvrt.html).

Usage

buildVRT(vrt_filename, input_rasters, cl_arg = NULL, quiet = FALSE)

Arguments

Details

Several command-line options are described in the GDAL documentation at the URL above. By default, the input files are considered as tiles of a larger mosaic and the VRT file has as many bands as one of the input files. Alternatively, the -separate argument can be used to put each input raster into a separate band in the VRT dataset.

Some amount of checks are done to assure that all files that will be put in the resulting VRT have similar characteristics: number of bands, projection, color interpretation.... If not, files that do not match the common characteristics will be skipped. (This is true in the default mode for virtual mosaicing, and not when using the -separate option).

In a virtual mosaic, if there is spatial overlap between input rasters then the order of files appearing in the list of sources matter: files that are listed at the end are the ones from which the data will be fetched. Note that nodata will be taken into account to potentially fetch data from less priority datasets.

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

See Also

rasterToVRT()

Examples

# build a virtual 3-band RGB raster from individual Landsat band files
b4_file <- system.file("extdata/sr_b4_20200829.tif", package="gdalraster")
b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
b6_file <- system.file("extdata/sr_b6_20200829.tif", package="gdalraster")
band_files <- c(b6_file, b5_file, b4_file)
vrt_file <- file.path(tempdir(), "storml_b6_b5_b4.vrt")
buildVRT(vrt_file, band_files, cl_arg = "-separate")
ds <- new(GDALRaster, vrt_file)
ds$getRasterCount()
plot_raster(ds, nbands=3, main="Landsat 6-5-4 (vegetative analysis)")
ds$close()

vsi_unlink(vrt_file)

Description

calc() evaluates an R expression for each pixel in a raster layer or stack of layers. Each layer is defined by a raster filename, band number, and a variable name to use in the R expression. If not specified, band defaults to 1 for each input raster. Variable names default to LETTERS if not specified (A (layer 1), B (layer 2), ...). All of the input layers must have the same extent and cell size. The projection will be read from the first raster in the list of inputs. Individual pixel coordinates are also available as variables in the R expression, as either x/y in the raster projected coordinate system or inverse projected longitude/latitude. Multiband output is supported as of gdalraster 1.11.0.

Usage

calc(
  expr,
  rasterfiles,
  bands = NULL,
  var.names = NULL,
  dstfile = tempfile("rastcalc", fileext = ".tif"),
  fmt = NULL,
  dtName = "Int16",
  out_band = NULL,
  options = NULL,
  nodata_value = NULL,
  setRasterNodataValue = FALSE,
  usePixelLonLat = NULL,
  write_mode = "safe",
  quiet = FALSE
)

Arguments

Details

The variables in expr are vectors of length raster xsize (row vectors of the input raster layer(s)). The expression should return a vector also of length raster xsize (an output row). Four special variable names are available in expr: pixelX and pixelY provide pixel center coordinates in projection units. pixelLon and pixelLat can also be used, in which case the pixel x/y coordinates will be inverse projected to longitude/latitude (in the same geographic coordinate system used by the input projection, which is read from the first input raster). Note that inverse projection adds computation time.

To refer to specific bands in a multi-band input file, repeat the filename in rasterfiles and specify corresponding band numbers in bands, along with optional variable names in var.names, for example,

rasterfiles = c("multiband.tif", "multiband.tif")
bands = c(4, 5)
var.names = c("B4", "B5")

Output will be written to dstfile. To update a file that already exists, set write_mode = "update" and set out_band to an existing band number(s) in dstfile (new bands cannot be created in dstfile).

To write multiband output, expr must return a vector of values interleaved by band. This is equivalent to, and can also be returned as, a matrix m with nrow(m) equal to length() of an input vector, and ncol(m) equal to the number of output bands. In matrix form, each column contains a vector of output values for a band. length(m) must be equal to the length() of an input vector multiplied by length(out_band). The dimensions described above are assumed and not read from the return value of expr.

Value

Returns the output filename invisibly.

See Also

GDALRaster-class, combine(), rasterToVRT()

Examples

## Using pixel longitude/latitude

# Hopkins bioclimatic index (HI) as described in:
# Bechtold, 2004, West. J. Appl. For. 19(4):245-251.
# Integrates elevation, latitude and longitude into an index of the
# phenological occurrence of springtime. Here it is relativized to
# mean values for an eight-state region in the western US.
# Positive HI means spring is delayed by that number of days relative
# to the reference position, while negative values indicate spring is
# advanced. The original equation had elevation units as feet, so
# converting m to ft in `expr`.

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")

# expression to calculate HI
expr <- "round( ((ELEV_M * 3.281 - 5449) / 100) +
                ((pixelLat - 42.16) * 4) +
                ((-116.39 - pixelLon) * 1.25) )"

# calc() writes to a tempfile by default
hi_file <- calc(expr = expr,
                rasterfiles = elev_file,
                var.names = "ELEV_M",
                dtName = "Int16",
                nodata_value = -32767,
                setRasterNodataValue = TRUE)

ds <- new(GDALRaster, hi_file)
# min, max, mean, sd
ds$getStatistics(band=1, approx_ok=FALSE, force=TRUE)
ds$close()
deleteDataset(hi_file)


## Calculate normalized difference vegetation index (NDVI)

# Landast band 4 (red) and band 5 (near infrared):
b4_file <- system.file("extdata/sr_b4_20200829.tif", package="gdalraster")
b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")

expr <- "((B5 * 0.0000275 - 0.2) - (B4 * 0.0000275 - 0.2)) /
         ((B5 * 0.0000275 - 0.2) + (B4 * 0.0000275 - 0.2))"
ndvi_file <- calc(expr = expr,
                  rasterfiles = c(b4_file, b5_file),
                  var.names = c("B4", "B5"),
                  dtName = "Float32",
                  nodata_value = -32767,
                  setRasterNodataValue = TRUE)

ds <- new(GDALRaster, ndvi_file)
ds$getStatistics(band=1, approx_ok=FALSE, force=TRUE)
ds$close()
deleteDataset(ndvi_file)


## Reclassify a variable by rule set

# Combine two raster layers and look for specific combinations. Then
# recode to a new value by rule set.
#
# Based on example in:
#   Stratton, R.D. 2009. Guidebook on LANDFIRE fuels data acquisition,
#   critique, modification, maintenance, and model calibration.
#   Gen. Tech. Rep. RMRS-GTR-220. U.S. Department of Agriculture,
#   Forest Service, Rocky Mountain Research Station. 54 p.
# Context: Refine national-scale fuels data to improve fire simulation
#   results in localized applications.
# Issue: Areas with steep slopes (40+ degrees) were mapped as
#   GR1 (101; short, sparse dry climate grass) and
#   GR2 (102; low load, dry climate grass) but were not carrying fire.
# Resolution: After viewing these areas in Google Earth,
#   NB9 (99; bare ground) was selected as the replacement fuel model.

# look for combinations of slope >= 40 and FBFM 101 or 102
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
rasterfiles <- c(lcp_file, lcp_file)
var.names <- c("SLP", "FBFM")
bands <- c(2, 4)
tbl <- combine(rasterfiles, var.names, bands)
nrow(tbl)
tbl_subset <- subset(tbl, SLP >= 40 & FBFM %in% c(101,102))
print(tbl_subset)       # twelve combinations meet the criteria
sum(tbl_subset$count)   # 85 total pixels

# recode these pixels to 99 (bare ground)
# the LCP driver does not support in-place write so make a copy as GTiff
tif_file <- file.path(tempdir(), "storml_lndscp.tif")
createCopy("GTiff", tif_file, lcp_file)

expr <- "ifelse( SLP >= 40 & FBFM %in% c(101,102), 99, FBFM)"
calc(expr = expr,
     rasterfiles = c(lcp_file, lcp_file),
     bands = c(2, 4),
     var.names = c("SLP", "FBFM"),
     dstfile = tif_file,
     out_band = 4,
     write_mode = "update")

# verify the ouput
rasterfiles <- c(tif_file, tif_file)
tbl <- combine(rasterfiles, var.names, bands)
tbl_subset <- subset(tbl, SLP >= 40 & FBFM %in% c(101,102))
print(tbl_subset)
sum(tbl_subset$count)

# if LCP file format is needed:
# createCopy("LCP", "storml_edited.lcp", tif_file)

deleteDataset(tif_file)

Description

CmbTable implements a hash table having a vector of integers as the key, and the count of occurrences of each unique integer combination as the value. A unique ID is assigned to each unique combination of input values.

Arguments

Value

An object of class CmbTable. Contains a hash table having a vector of keyLen integers as the key and the count of occurrences of each unique integer combination as the value, along with methods that operate on the table as described in Details. CmbTable is a C++ class exposed directly to R (via RCPP_EXPOSED_CLASS). Methods of the class are accessed in R using the $ operator.

Usage

## Constructors
cmb <- new(CmbTable, keyLen)
# or, with variable names
cmb <- new(CmbTable, keyLen, varNames)

## Methods (see Details)
cmb$update(int_cmb, incr)
cmb$updateFromMatrix(int_cmbs, incr)
cmb$updateFromMatrixByRow(int_cmbs, incr)
cmb$asDataFrame()
cmb$asMatrix()

Details

new(CmbTable, keyLen) Constructor. Variable names will be assigned as V1, V2, .... Returns an object of class CmbTable.

new(CmbTable, keyLen, varNames) Alternate constructor to specify variable names. Returns an object of class CmbTable.

$update(int_cmb, incr) Updates the hash table for the integer combination in the numeric vector int_cmb (coerced to integer by truncation). If this combination exists in the table, its count will be incremented by incr. If the combination is not found in the table, it will be inserted with count set to incr. Returns the unique ID assigned to this combination. Combination IDs are sequential integers starting at 1.

$updateFromMatrix(int_cmbs, incr) This method is the same as ⁠$update()⁠ but for a numeric matrix of integer combinations int_cmbs (coerced to integer by truncation). The matrix is arranged with each column vector forming an integer combination. For example, the rows of the matrix could be one row each from a set of keyLen rasters all read at the same extent and pixel resolution (i.e., row-by-row raster overlay). The method calls ⁠$update()⁠ on each combination (each column of int_cmbs), incrementing count by incr for existing combinations, or inserting new combinations with count set to incr. Returns a numeric vector of length ncol(int_cmbs) containing the IDs assigned to the combinations.

$updateFromMatrixByRow(int_cmbs, incr) This method is the same as ⁠$updateFromMatrix()⁠ above except the integer combinations are in rows of the matrix int_cmbs (columns are the variables). The method calls ⁠$update()⁠ on each combination (each row of int_cmbs), incrementing count by incr for existing combinations, or inserting new combinations with count set to incr. Returns a numeric vector of length nrow(int_cmbs) containing the IDs assigned to the combinations.

$asDataFrame() Returns the CmbTable as a data frame with column cmbid containing the unique combination IDs, column count containing the counts of occurrences, and keyLen columns (with names from varNames) containing the integer values comprising each unique combination.

$asMatrix() Returns the CmbTable as a matrix with column 1 (cmbid) containing the unique combination IDs, column 2 (count) containing the counts of occurrences, and columns 3:keyLen+2 (with names from varNames) containing the integer values comprising each unique combination.

Examples

m <- matrix(c(1,2,3,1,2,3,4,5,6,1,3,2,4,5,6,1,1,1), 3, 6, byrow=FALSE)
rownames(m) <- c("layer1", "layer2", "layer3")
print(m)
cmb <- new(CmbTable, 3, rownames(m))
cmb$updateFromMatrix(m, 1)
cmb$asDataFrame()
cmb$update(c(4,5,6), 1)
cmb$update(c(1,3,5), 1)
cmb$asDataFrame()

# same as above but matrix arranged with integer combinations in the rows
m <- matrix(c(1,2,3,1,2,3,4,5,6,1,3,2,4,5,6,1,1,1), 6, 3, byrow=TRUE)
colnames(m) <- c("V1", "V2", "V3")
print(m)
cmb <- new(CmbTable, 3)
cmb$updateFromMatrixByRow(m, 1)
cmb$asDataFrame()
cmb$update(c(4,5,6), 1)
cmb$update(c(1,3,5), 1)
cmb$asDataFrame()

Description

combine() overlays multiple rasters so that a unique ID is assigned to each unique combination of input values. The input raster layers typically have integer data types (floating point will be coerced to integer by truncation), and must have the same projection, extent and cell size. Pixel counts for each unique combination are obtained, and combination IDs are optionally written to an output raster.

Usage

combine(
  rasterfiles,
  var.names = NULL,
  bands = NULL,
  dstfile = NULL,
  fmt = NULL,
  dtName = "UInt32",
  options = NULL,
  quiet = FALSE
)

Arguments

Details

To specify input raster layers that are bands of a multi-band raster file, repeat the filename in rasterfiles and provide the corresponding band numbers in bands. For example:

rasterfiles <- c("multi-band.tif", "multi-band.tif", "other.tif")
bands <- c(4, 5, 1)
var.names <- c("multi_b4", "multi_b5", "other")

rasterToVRT() provides options for virtual clipping, resampling and pixel alignment, which may be helpful here if the input rasters are not already aligned on a common extent and cell size.

If an output raster of combination IDs is written, the user should verify that the number of combinations obtained did not exceed the range of the output data type. Combination IDs are sequential integers starting at 1. Typical output data types are the unsigned types: Byte (0 to 255), UInt16 (0 to 65,535) and UInt32 (the default, 0 to 4,294,967,295).

Value

A data frame with column cmbid containing the combination IDs, column count containing the pixel counts for each combination, and length(rasterfiles) columns named var.names containing the integer values comprising each unique combination.

See Also

CmbTable-class, GDALRaster-class, calc(), rasterToVRT()

buildRAT() to compute a table of the unique pixel values and their counts for a single raster layer

Examples

evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
evc_file <- system.file("extdata/storml_evc.tif", package="gdalraster")
evh_file <- system.file("extdata/storml_evh.tif", package="gdalraster")
rasterfiles <- c(evt_file, evc_file, evh_file)
var.names <- c("veg_type", "veg_cov", "veg_ht")
tbl <- combine(rasterfiles, var.names)
nrow(tbl)
tbl <- tbl[order(-tbl$count),]
head(tbl, n = 20)

# combine two bands from a multi-band file and write the combination IDs
# to an output raster
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
rasterfiles <- c(lcp_file, lcp_file)
bands <- c(4, 5)
var.names <- c("fbfm", "tree_cov")
cmb_file <- file.path(tempdir(), "fbfm_cov_cmbid.tif")
opt <- c("COMPRESS=LZW")
tbl <- combine(rasterfiles, var.names, bands, cmb_file, options = opt)
head(tbl)
ds <- new(GDALRaster, cmb_file)
ds$info()
ds$close()
deleteDataset(cmb_file)

Description

copyDatasetFiles() copies all the files associated with a dataset. Wrapper for GDALCopyDatasetFiles() in the GDAL API.

Usage

copyDatasetFiles(new_filename, old_filename, format = "")

Arguments

Value

Logical TRUE if no error or FALSE on failure.

Note

If format is set to an empty string "" (the default) then the function will try to identify the driver from old_filename. This is done internally in GDAL by invoking the Identify method of each registered GDALDriver in turn. The first driver that successful identifies the file name will be returned. An error is raised if a format cannot be determined from the passed file name.

See Also

GDALRaster-class, create(), createCopy(), deleteDataset(), renameDataset(), vsi_copy_file()

Examples

lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds <- new(GDALRaster, lcp_file)
ds$getFileList()
ds$close()

lcp_tmp <- file.path(tempdir(), "storm_lake_copy.lcp")
copyDatasetFiles(lcp_tmp, lcp_file)
ds_copy <- new(GDALRaster, lcp_tmp)
ds_copy$getFileList()
ds_copy$close()

deleteDataset(lcp_tmp)

Description

create() makes an empty raster in the specified format.

Usage

create(format, dst_filename, xsize, ysize, nbands, dataType, options = NULL)

Arguments

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

See Also

GDALRaster-class, createCopy(), rasterFromRaster(), getCreationOptions()

Examples

new_file <- file.path(tempdir(), "newdata.tif")
create(format="GTiff", dst_filename=new_file, xsize=143, ysize=107,
       nbands=1, dataType="Int16")
ds <- new(GDALRaster, new_file, read_only=FALSE)
## EPSG:26912 - NAD83 / UTM zone 12N
ds$setProjection(epsg_to_wkt(26912))
gt <- c(323476.1, 30, 0, 5105082.0, 0, -30)
ds$setGeoTransform(gt)
ds$setNoDataValue(band = 1, -9999)
ds$fillRaster(band = 1, -9999, 0)
## ...
## close the dataset when done
ds$close()

deleteDataset(new_file)

Description

createColorRamp() is a wrapper for GDALCreateColorRamp() in the GDAL API. It automatically creates a color ramp from one color entry to another. Output is an integer matrix in color table format for use with GDALRaster$setColorTable().

Usage

createColorRamp(
  start_index,
  start_color,
  end_index,
  end_color,
  palette_interp = "RGB"
)

Arguments

Value

Integer matrix with five columns containing the color ramp from start_index to end_index, with raster index values in column 1 and color entries in columns 2:5).

Note

createColorRamp() could be called several times, using rbind() to combine multiple ramps into the same color table. Possible duplicate rows in the resulting table are not a problem when used in GDALRaster$setColorTable() (i.e., when end_color of one ramp is the same as start_color of the next ramp).

See Also

GDALRaster$getColorTable(), GDALRaster$getPaletteInterp()

Examples

# create a color ramp for tree canopy cover percent
# band 5 of an LCP file contains canopy cover
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds <- new(GDALRaster, lcp_file)
ds$getDescription(band=5)
ds$getMetadata(band=5, domain="")
ds$close()

# create a GTiff file with Byte data type for the canopy cover band
# recode nodata -9999 to 255
tcc_file <- calc(expr = "ifelse(CANCOV == -9999, 255, CANCOV)",
                 rasterfiles = lcp_file,
                 bands = 5,
                 var.names = "CANCOV",
                 fmt = "GTiff",
                 dtName = "Byte",
                 nodata_value = 255,
                 setRasterNodataValue = TRUE)

ds_tcc <- new(GDALRaster, tcc_file, read_only=FALSE)

# create a color ramp from 0 to 100 and set as the color table
colors <- createColorRamp(start_index = 0,
                          start_color = c(211, 211, 211),
                          end_index = 100,
                          end_color = c(0, 100, 0))

print(colors)
ds_tcc$setColorTable(band=1, col_tbl=colors, palette_interp="RGB")
ds_tcc$setRasterColorInterp(band=1, col_interp="Palette")

# close and re-open the dataset in read_only mode
ds_tcc$open(read_only=TRUE)

plot_raster(ds_tcc, interpolate=FALSE, legend=TRUE,
            main="Storm Lake Tree Canopy Cover (%)")
ds_tcc$close()

deleteDataset(tcc_file)

Description

createCopy() copies a raster dataset, optionally changing the format. The extent, cell size, number of bands, data type, projection, and geotransform are all copied from the source raster.

Usage

createCopy(
  format,
  dst_filename,
  src_filename,
  strict = FALSE,
  options = NULL,
  quiet = FALSE
)

Arguments

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

See Also

GDALRaster-class, create(), rasterFromRaster(), getCreationOptions(), translate()

Examples

lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
tif_file <- file.path(tempdir(), "storml_lndscp.tif")
opt <- c("COMPRESS=LZW")
createCopy(format="GTiff", dst_filename=tif_file, src_filename=lcp_file,
           options=opt)
file.size(lcp_file)
file.size(tif_file)
ds <- new(GDALRaster, tif_file, read_only=FALSE)
ds$getMetadata(band=0, domain="IMAGE_STRUCTURE")
for (band in 1:ds$getRasterCount())
    ds$setNoDataValue(band, -9999)
ds$getStatistics(band=1, approx_ok=FALSE, force=TRUE)
ds$close()

deleteDataset(tif_file)

Description

These values are used in dem_proc() as the default processing options:

    list(
         hillshade =    c("-z", "1", "-s", "1", "-az", "315",
                          "-alt", "45", "-alg", "Horn",
                          "-combined", "-compute_edges"),
         slope =        c("-s", "1", "-alg", "Horn", "-compute_edges"),
         aspect =       c("-alg", "Horn", "-compute_edges"),
         color_relief = character(),
         TRI =          c("-alg", "Riley", "-compute_edges"),
         TPI =          c("-compute_edges"),
         roughness =    c("-compute_edges")
         )

Usage

DEFAULT_DEM_PROC

Format

An object of class list of length 7.

See Also

dem_proc()

https://gdal.org/programs/gdaldem.html for a description of all available command-line options for each processing mode

Description

These values are currently used in gdalraster when a nodata value is needed but has not been specified:

    list("Byte" = 255, "Int8" = -128,
         "UInt16" = 65535, "Int16" = -32767,
         "UInt32" = 4294967293, "Int32" = -2147483647,
         "Float32" = -99999.0, "Float64" = -99999.0)

Usage

DEFAULT_NODATA

Format

An object of class list of length 8.

Description

deleteDataset() will attempt to delete the named dataset in a format specific fashion. Full featured drivers will delete all associated files, database objects, or whatever is appropriate. The default behavior when no format specific behavior is provided is to attempt to delete all the files that would be returned by GDALRaster$getFileList() on the dataset. The named dataset should not be open in any existing GDALRaster objects when deleteDataset() is called. Wrapper for GDALDeleteDataset() in the GDAL API.

Usage

deleteDataset(filename, format = "")

Arguments

Value

Logical TRUE if no error or FALSE on failure.

Note

If format is set to an empty string "" (the default) then the function will try to identify the driver from filename. This is done internally in GDAL by invoking the Identify method of each registered GDALDriver in turn. The first driver that successful identifies the file name will be returned. An error is raised if a format cannot be determined from the passed file name.

See Also

GDALRaster-class, create(), createCopy(), copyDatasetFiles(), renameDataset()

Examples

b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
b5_tmp <- file.path(tempdir(), "b5_tmp.tif")
file.copy(b5_file,  b5_tmp)

ds <- new(GDALRaster, b5_tmp)
ds$buildOverviews("BILINEAR", levels = c(2, 4, 8), bands = c(1))
files <- ds$getFileList()
print(files)
ds$close()
file.exists(files)
deleteDataset(b5_tmp)
file.exists(files)

Description

dem_proc() generates DEM derivatives from an input elevation raster. This function is a wrapper for the gdaldem command-line utility. See https://gdal.org/programs/gdaldem.html for details.

Usage

dem_proc(
  mode,
  srcfile,
  dstfile,
  mode_options = DEFAULT_DEM_PROC[[mode]],
  color_file = NULL,
  quiet = FALSE
)

Arguments

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

Note

Band 1 of the source elevation raster is read by default, but this can be changed by including a -b command-line argument in mode_options. See the documentation for gdaldem for a description of all available options for each processing mode.

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
slp_file <- file.path(tempdir(), "storml_slp.tif")
dem_proc("slope", elev_file, slp_file)

deleteDataset(slp_file)

Description

displayRAT() generates a presentation table. Colors are shown if the Raster Attribute Table contains RGB columns. This function requires package gt.

Usage

displayRAT(tbl, title = "Raster Attribute Table")

Arguments

Value

An object of class "gt_tbl" (i.e., a table created with gt::gt()).

See Also

buildRAT(), GDALRaster$getDefaultRAT()

vignette("raster-attribute-tables")

Examples

# see examples for `buildRAT()`

Description

dump_open_datasets() dumps a list of all open datasets (shared or not) to the console. This function is primarily intended to assist in debugging "dataset leaks" and reference counting issues. The information reported includes the dataset name, referenced count, shared status, driver name, size, and band count. This a wrapper for GDALDumpOpenDatasets() with output to the console.

Usage

dump_open_datasets()

Value

Number of open datasets.

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)
dump_open_datasets()
ds2 <- new(GDALRaster, elev_file)
dump_open_datasets()
# open without using shared mode
ds3 <- new(GDALRaster, elev_file, read_only = TRUE,
           open_options = NULL, shared = FALSE)
dump_open_datasets()
ds$close()
dump_open_datasets()
ds2$close()
dump_open_datasets()
ds3$close()
dump_open_datasets()

Description

epsg_to_wkt() exports the spatial reference for an EPSG code to WKT format.

Usage

epsg_to_wkt(epsg, pretty = FALSE)

Arguments

Details

As of GDAL 3.0, the default format for WKT export is OGC WKT 1. The WKT version can be overridden by using the OSR_WKT_FORMAT configuration option (see set_config_option()). Valid values are one of: SFSQL, WKT1_SIMPLE, WKT1, WKT1_GDAL, WKT1_ESRI, WKT2_2015, WKT2_2018, WKT2, DEFAULT. If SFSQL, a WKT1 string without AXIS, TOWGS84, AUTHORITY or EXTENSION node is returned. If WKT1_SIMPLE, a WKT1 string without AXIS, AUTHORITY or EXTENSION node is returned. WKT1 is an alias of WKT1_GDAL. WKT2 will default to the latest revision implemented (currently WKT2_2018). WKT2_2019 can be used as an alias of WKT2_2018 since GDAL 3.2

Value

Character string containing OGC WKT.

See Also

srs_to_wkt()

Examples

epsg_to_wkt(5070)
writeLines(epsg_to_wkt(5070, pretty=TRUE))
set_config_option("OSR_WKT_FORMAT", "WKT2")
writeLines(epsg_to_wkt(5070, pretty=TRUE))
set_config_option("OSR_WKT_FORMAT", "")

Description

fillNodata() is a wrapper for GDALFillNodata() in the GDAL Algorithms API. This algorithm will interpolate values for all designated nodata pixels (pixels having an intrinsic nodata value, or marked by zero-valued pixels in the optional raster specified in mask_file). For each nodata pixel, a four direction conic search is done to find values to interpolate from (using inverse distance weighting). Once all values are interpolated, zero or more smoothing iterations (3x3 average filters on interpolated pixels) are applied to smooth out artifacts.

Usage

fillNodata(
  filename,
  band,
  mask_file = "",
  max_dist = 100,
  smooth_iterations = 0L,
  quiet = FALSE
)

Arguments

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

Note

The input raster will be modified in place. It should not be open in a GDALRaster object while processing with fillNodata().

Examples

## fill nodata edge pixels in the elevation raster
elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")

## get count of nodata
tbl <- buildRAT(elev_file)
head(tbl)
tbl[is.na(tbl$VALUE),]

## make a copy that will be modified
mod_file <- file.path(tempdir(), "storml_elev_fill.tif")
file.copy(elev_file,  mod_file)

fillNodata(mod_file, band=1)

mod_tbl = buildRAT(mod_file)
head(mod_tbl)
mod_tbl[is.na(mod_tbl$VALUE),]

deleteDataset(mod_file)

Description

footprint() is a wrapper of the gdal_footprint command-line utility (see https://gdal.org/programs/gdal_footprint.html). The function can be used to compute the footprint of a raster file, taking into account nodata values (or more generally the mask band attached to the raster bands), and generating polygons/multipolygons corresponding to areas where pixels are valid, and write to an output vector file. Refer to the GDAL documentation at the URL above for a list of command-line arguments that can be passed in cl_arg. Requires GDAL >= 3.8.

Usage

footprint(src_filename, dst_filename, cl_arg = NULL)

Arguments

Details

Post-vectorization geometric operations are applied in the following order:

• optional splitting (-split_polys)

• optional densification (-densify)

• optional reprojection (-t_srs)

• optional filtering by minimum ring area (-min_ring_area)

• optional application of convex hull (-convex_hull)

• optional simplification (-simplify)

• limitation of number of points (-max_points)

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

See Also

polygonize()

Examples

evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
out_file <- file.path(tempdir(), "storml.geojson")

# Requires GDAL >= 3.8
if (as.integer(gdal_version()[2]) >= 3080000) {
  # command-line arguments for gdal_footprint
  args <- c("-t_srs", "EPSG:4326")
  footprint(evt_file, out_file, args)

  deleteDataset(out_file)
}

Description

g_area() computes the area for a LinearRing, Polygon or MultiPolygon. Undefined for all other geometry types (returns zero).

Usage

g_area(wkt)

Arguments

Value

Numeric scalar. Area of the geometry or 0.

Note

LinearRing is a non-standard geometry type, used in GDAL just for geometry creation.

Description

These functions implement operations on pairs of geometries in OGC WKT format.

Usage

g_intersection(this_geom, other_geom)

g_union(this_geom, other_geom)

g_difference(this_geom, other_geom)

g_sym_difference(this_geom, other_geom)

Arguments

Details

These functions use the GEOS library via GDAL headers.

g_intersection() returns a new geometry which is the region of intersection of the two geometries operated on. g_intersects() can be used to test if two geometries intersect.

g_union() returns a new geometry which is the region of union of the two geometries operated on.

g_difference() returns a new geometry which is the region of this geometry with the region of the other geometry removed.

g_sym_difference() returns a new geometry which is the symmetric difference of this geometry and the other geometry (union minus intersection).

Value

Character string. The resulting geometry as OGC WKT.

Note

Geometry validity is not checked. In case you are unsure of the validity of the input geometries, call g_is_valid() before, otherwise the result might be wrong.

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)
g1 <- ds$bbox() |> bbox_to_wkt()
ds$close()

g2 <- "POLYGON ((327381.9 5104541.2, 326824.0 5104092.5, 326708.8 5103182.9,
  327885.2 5102612.9, 329334.5 5103322.4, 329304.2 5104474.5,328212.7
  5104656.4, 328212.7 5104656.4, 327381.9 5104541.2))"

# see spatial predicate defintions at https://en.wikipedia.org/wiki/DE-9IM
g_intersects(g1, g2)  # TRUE
g_overlaps(g1, g2)  # TRUE
# therefore,
g_contains(g1, g2)  # FALSE

g_sym_difference(g1, g2) |> g_area()

g3 <- g_intersection(g1, g2)
g4 <- g_union(g1, g2)
g_difference(g4, g3) |> g_area()

Description

These functions implement tests for pairs of geometries in OGC WKT format.

Usage

g_intersects(this_geom, other_geom)

g_disjoint(this_geom, other_geom)

g_touches(this_geom, other_geom)

g_contains(this_geom, other_geom)

g_within(this_geom, other_geom)

g_crosses(this_geom, other_geom)

g_overlaps(this_geom, other_geom)

g_equals(this_geom, other_geom)

Arguments

Details

These functions use the GEOS library via GDAL headers.

g_intersects() tests whether two geometries intersect.

g_disjoint() tests if this geometry and the other geometry are disjoint.

g_touches() tests if this geometry and the other geometry are touching.

g_contains() tests if this geometry contains the other geometry.

g_within() tests if this geometry is within the other geometry.

g_crosses() tests if this geometry and the other geometry are crossing.

g_overlaps() tests if this geometry and the other geometry overlap, that is, their intersection has a non-zero area (they have some but not all points in common).

g_equals() tests whether two geometries are equivalent. The GDAL documentation says: "This operation implements the SQL/MM ST_OrderingEquals() operation. The comparison is done in a structural way, that is to say that the geometry types must be identical, as well as the number and ordering of sub-geometries and vertices. Or equivalently, two geometries are considered equal by this method if their WKT/WKB representation is equal. Note: this must be distinguished from equality in a spatial way."

Value

Logical scalar

Note

Geometry validity is not checked. In case you are unsure of the validity of the input geometries, call g_is_valid() before, otherwise the result might be wrong.

See Also

https://en.wikipedia.org/wiki/DE-9IM

Description

g_buffer() builds a new geometry containing the buffer region around the geometry on which it is invoked. The buffer is a polygon containing the region within the buffer distance of the original geometry.

Usage

g_buffer(wkt, dist, quad_segs = 30L)

Arguments

Value

Character string for an OGC WKT polygon.

See Also

bbox_from_wkt(), bbox_to_wkt()

Examples

g_buffer(wkt = "POINT (0 0)", dist = 10)

Description

g_centroid() returns a vector of point X, point Y.

Usage

g_centroid(wkt)

Arguments

Details

The GDAL documentation states "This method relates to the SFCOM ISurface::get_Centroid() method however the current implementation based on GEOS can operate on other geometry types such as multipoint, linestring, geometrycollection such as multipolygons. OGC SF SQL 1.1 defines the operation for surfaces (polygons). SQL/MM-Part 3 defines the operation for surfaces and multisurfaces (multipolygons)."

Value

Numeric vector of length 2 containing the centroid (X, Y).

Description

g_distance() returns the distance between two geometries or -1 if an error occurs. Returns the shortest distance between the two geometries. The distance is expressed into the same unit as the coordinates of the geometries.

Usage

g_distance(this_geom, other_geom)

Arguments

Value

Numeric. Distance or '-1' if an error occurs.

Note

Geometry validity is not checked. In case you are unsure of the validity of the input geometries, call g_is_valid() before, otherwise the result might be wrong.

Examples

g_distance("POINT (0 0)", "POINT (5 12)")

Description

g_is_empty() tests whether a geometry has no points.

Usage

g_is_empty(wkt)

Arguments

Value

logical scalar. TRUE if the geometry has no points, otherwise FALSE.

Examples

g1 <- "POLYGON ((0 0, 10 10, 10 0, 0 0))"
g2 <- "POLYGON ((5 1, 9 5, 9 1, 5 1))"
g_difference(g2, g1) |> g_is_empty()

Description

g_is_valid() tests whether a geometry is valid.

Usage

g_is_valid(wkt)

Arguments

Value

logical scalar. TRUE if the geometry is valid, otherwise FALSE.

Examples

g1 <- "POLYGON ((0 0, 10 10, 10 0, 0 0))"
g_is_valid(g1)

g2 <- "POLYGON ((0 0, 10 10, 10 0, 0 1))"
g_is_valid(g2)

Description

g_length() computes the length for LineString or MultiCurve objects. Undefined for all other geometry types (returns zero).

Usage

g_length(wkt)

Arguments

Value

Numeric scalar. Length of the geometry or 0.

Examples

g_length("LINESTRING (0 0, 3 4)")

Description

g_name() returns the name for this geometry type in well known text format.

Usage

g_name(wkt)

Arguments

Value

WKT name for this geometry type.

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)
bbox_to_wkt(ds$bbox()) |> g_name()
ds$close()

Description

g_transform() will transform the coordinates of a geometry from their current spatial reference system to a new target spatial reference system. Normally this means reprojecting the vectors, but it could include datum shifts, and changes of units.

Usage

g_transform(
  wkt,
  srs_from,
  srs_to,
  wrap_date_line = FALSE,
  date_line_offset = 10L
)

Arguments

Value

Character string for a transformed OGC WKT geometry.

Note

This function uses the OGR_GeomTransformer_Create() and OGR_GeomTransformer_Transform() functions in the GDAL API: "This is an enhanced version of OGR_G_Transform(). When reprojecting geometries from a Polar Stereographic projection or a projection naturally crossing the antimeridian (like UTM Zone 60) to a geographic CRS, it will cut geometries along the antimeridian. So a LineString might be returned as a MultiLineString."

The wrap_date_line = TRUE option might be specified for circumstances to correct geometries that incorrectly go from a longitude on a side of the antimeridian to the other side, e.g., ⁠LINESTRING (-179 0,179 0)⁠ will be transformed to ⁠MULTILINESTRING ((-179 0,-180 0),(180 0,179 0))⁠. For that use case, srs_to might be the same as srs_from.

See Also

bbox_transform()

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)

# the convenience function bbox_transform() does this:
bbox_to_wkt(ds$bbox()) |>
  g_transform(ds$getProjection(), epsg_to_wkt(4326)) |>
  bbox_from_wkt()

ds$close()

# correct geometries that incorrectly go from a longitude on a side of the
# antimeridian to the other side
geom <- "LINESTRING (-179 0,179 0)"
srs <- epsg_to_wkt(4326)
g_transform(geom, srs, srs, wrap_date_line = TRUE)

Description

gdal_formats() returns a table of the supported raster and vector formats, with information about the capabilities of each format driver.

Usage

gdal_formats(format = "")

Arguments

Value

A data frame containing the format short name, long name, raster (logical), vector (logical), read/write flag (ro is read-only, w supports CreateCopy, ⁠w+⁠ supports Create), virtual I/O supported (logical), and subdatasets (logical).

Note

Virtual I/O refers to operations on GDAL Virtual File Systems. See https://gdal.org/user/virtual_file_systems.html#virtual-file-systems.

Examples

nrow(gdal_formats())
head(gdal_formats())

gdal_formats("GPKG")

Description

gdal_version() returns runtime version information.

Usage

gdal_version()

Value

Character vector of length four containing:

• "–version" - one line version message, e.g., “GDAL 3.6.3, released 2023/03/12”

• "GDAL_VERSION_NUM" - formatted as a string, e.g., “3060300” for GDAL 3.6.3.0

• "GDAL_RELEASE_DATE" - formatted as a string, e.g., “20230312”

• "GDAL_RELEASE_NAME" - e.g., “3.6.3”

Examples

gdal_version()

Description

GDALRaster provides an interface for accessing a raster dataset via GDAL and calling methods on the underlying GDALDataset, GDALDriver and GDALRasterBand objects. See https://gdal.org/api/index.html for details of the GDAL Raster API.

Arguments

Value

An object of class GDALRaster which contains a pointer to the opened dataset, and methods that operate on the dataset as described in Details. GDALRaster is a C++ class exposed directly to R (via RCPP_EXPOSED_CLASS). Fields and methods of the class are accessed using the $ operator. The read/write fields can be used for per-object settings.

Usage

## Constructors
# read-only by default:
ds <- new(GDALRaster, filename)
# for update access:
ds <- new(GDALRaster, filename, read_only = FALSE)
# to use dataset open options:
ds <- new(GDALRaster, filename, read_only = TRUE|FALSE, open_options)
# to open without shared mode:
new(GDALRaster, filename, read_only, open_options, shared = FALSE)

## Read/write fields (see Details)
ds$infoOptions
ds$quiet
ds$readByteAsRaw

## Methods (see Details)
ds$getFilename()
ds$open(read_only)
ds$isOpen()
ds$getFileList()

ds$info()
ds$infoAsJSON()

ds$getDriverShortName()
ds$getDriverLongName()

ds$getRasterXSize()
ds$getRasterYSize()
ds$getGeoTransform()
ds$setGeoTransform(transform)
ds$getProjection()
ds$getProjectionRef()
ds$setProjection(projection)
ds$bbox()
ds$res()
ds$dim()
ds$apply_geotransform(col_row)
ds$get_pixel_line(xy)

ds$getRasterCount()
ds$getDescription(band)
ds$setDescription(band, desc)
ds$getBlockSize(band)
ds$getActualBlockSize(band, xblockoff, yblockoff)
ds$getOverviewCount(band)
ds$buildOverviews(resampling, levels, bands)
ds$getDataTypeName(band)
ds$getNoDataValue(band)
ds$setNoDataValue(band, nodata_value)
ds$deleteNoDataValue(band)
ds$getUnitType(band)
ds$setUnitType(band, unit_type)
ds$getScale(band)
ds$setScale(band, scale)
ds$getOffset(band)
ds$setOffset(band, offset)
ds$getRasterColorInterp(band)
ds$setRasterColorInterp(band, col_interp)

ds$getMinMax(band, approx_ok)
ds$getStatistics(band, approx_ok, force)
ds$clearStatistics()
ds$getHistogram(band, min, max, num_buckets, incl_out_of_range, approx_ok)
ds$getDefaultHistogram(band, force)

ds$getMetadata(band, domain)
ds$getMetadataItem(band, mdi_name, domain)
ds$setMetadataItem(band, mdi_name, mdi_value, domain)
ds$getMetadataDomainList(band)

ds$read(band, xoff, yoff, xsize, ysize, out_xsize, out_ysize)
ds$write(band, xoff, yoff, xsize, ysize, rasterData)
ds$fillRaster(value, ivalue)

ds$getColorTable(band)
ds$getPaletteInterp(band)
ds$setColorTable(band, col_tbl, palette_interp)

ds$getDefaultRAT(band)
ds$setDefaultRAT(band, df)

ds$flushCache()

ds$getChecksum(band, xoff, yoff, xsize, ysize)

ds$close()

Details

new(GDALRaster, filename, read_only) Constructor. Returns an object of class GDALRaster. read_only defaults to TRUE if not specified.

new(GDALRaster, filename, read_only, open_options) Alternate constructor for passing dataset open_options, a character vector of NAME=VALUE pairs. read_only is required for this form of the constructor, TRUE for read-only, or FALSE to open with write access. Returns an object of class GDALRaster.

new(GDALRaster, filename, read_only, open_options, shared) Alternate constructor for specifying the shared mode for dataset opening. shared defaults to TRUE but can be set to FALSE with this constructor (see Note). All parameters are required with this form of the constructor, but open_options can be NULL. Returns an object of class GDALRaster.

$infoOptions Read/write field. A character vector of command-line arguments to control the output of ⁠$info()⁠ and ⁠$infoAsJSON()⁠ (see below). Defaults to character(0). Can be set to a vector of strings specifying arguments to the gdalinfo command-line utility, e.g., c("-nomd", "-norat", "-noct"). Restore the default by setting to empty string ("") or character(0).

$quiet Read/write field. A logical value, FALSE by default. This field can be set to TRUE which will suppress various messages as well as progress reporting for potentially long-running processes such as building overviews and computation of statistics and histograms.

$readByteAsRaw Read/write field. A logical value, FALSE by default. This field can be set to TRUE which will affect the data type returned by ⁠$read()⁠ and read_ds(). When the underlying band data type is 'Byte' and readByteAsRaw is TRUE the output type will be raw rather than integer. See also the as_raw argument to read_ds() to control this in a non-persistent setting. If the underlying band data type is not Byte this setting has no effect.

$getFilename() Returns a character string containing the filename associated with this GDALRaster object (filename originally used to open the dataset).

$open(read_only) (Re-)opens the raster dataset on the existing filename. Use this method to open a dataset that has been closed using $close(). May be used to re-open a dataset with a different read/write access (read_only set to TRUE or FALSE). The method will first close an open dataset, so it is not required to call $close() explicitly in this case. No return value, called for side effects.

$isOpen() Returns logical indicating whether the associated raster dataset is open.

$getFileList() Returns a character vector of files believed to be part of this dataset. If it returns an empty string ("") it means there is believed to be no local file system files associated with the dataset (e.g., a virtual file system). The returned filenames will normally be relative or absolute paths depending on the path used to originally open the dataset.

$info() Prints various information about the raster dataset to the console (no return value, called for that side effect only). Equivalent to the output of the gdalinfo command-line utility (gdalinfo filename, if using the default infoOptions). See the field ⁠$infoOptions⁠ above for setting the arguments to gdalinfo.

$infoAsJSON() Returns information about the raster dataset as a JSON-formatted string. Equivalent to the output of the gdalinfo command-line utility (gdalinfo -json filename, if using the default infoOptions). See the field ⁠$infoOptions⁠ above for setting the arguments to gdalinfo.

$getDriverShortName() Returns the short name of the raster format driver.

$getDriverLongName() Returns the long name of the raster format driver.

$getRasterXSize() Returns the number of pixels along the x dimension.

$getRasterYSize() Returns the number of pixels along the y dimension.

$getGeoTransform() Returns the affine transformation coefficients for transforming between pixel/line raster space (column/row) and projection coordinate space (geospatial x/y). The return value is a numeric vector of length six. See https://gdal.org/tutorials/geotransforms_tut.html for details of the affine transformation. With 1-based indexing in R, the geotransform vector contains (in map units of the raster spatial reference system):

$setGeoTransform(transform) Sets the affine transformation coefficients on this dataset. transform is a numeric vector of length six. Returns logical TRUE on success or FALSE if the geotransform could not be set.

$getProjection() Returns the coordinate reference system of the raster as an OGC WKT format string. Equivalent to ds$getProjectionRef().

$getProjectionRef() Returns the coordinate reference system of the raster as an OGC WKT format string. An empty string is returned when a projection definition is not available.

$setProjection(projection) Sets the projection reference for this dataset. projection is a string in OGC WKT format. Returns logical TRUE on success or FALSE if the projection could not be set.

$bbox() Returns a numeric vector of length four containing the bounding box (xmin, ymin, xmax, ymax) assuming this is a north-up raster.

$res() Returns a numeric vector of length two containing the resolution (pixel width, pixel height as positive values) assuming this is a north-up raster.

$dim() Returns an integer vector of length three containing the raster dimensions. Equivalent to: c(ds$getRasterXSize(), ds$getRasterYSize(), ds$getRasterCount())

$apply_geotransform(col_row) Applies geotransform coefficients to raster coordinates in pixel/line space (column/row), converting into georeferenced (x/y) coordinates. col_row is a numeric matrix of raster col/row coordinates (or two-column data frame that will be coerced to numeric matrix). Returns a numeric matrix of geospatial x/y coordinates. See the stand-alone function of the same name (apply_geotransform()) for more info and examples.

$get_pixel_line(xy) Converts geospatial coordinates to pixel/line (raster column/row numbers). xy is a numeric matrix of geospatial x,y coordinates in the same spatial reference system as the raster (or two-column data frame that will be coerced to numeric matrix). Returns an integer matrix of raster pixel/line. See the stand-alone function of the same name (get_pixel_line()) for more info and examples.

$getRasterCount() Returns the number of raster bands on this dataset. For the methods described below that operate on individual bands, the band argument is the integer band number (1-based).

$getDescription(band) Returns a string containing the description for band. An empty string is returned if no description is set for the band. (Setting band = 0 will return the dataset-level description.)

$setDescription(band, desc) Sets a description for band. desc is the character string to set. No return value.

$getBlockSize(band) Returns an integer vector of length two (xsize, ysize) containing the "natural" block size of band. GDAL has a concept of the natural block size of rasters so that applications can organize data access efficiently for some file formats. The natural block size is the block size that is most efficient for accessing the format. For many formats this is simply a whole row in which case block xsize is the same as $getRasterXSize() and block ysize is 1. However, for tiled images block size will typically be the tile size. Note that the X and Y block sizes don't have to divide the image size evenly, meaning that right and bottom edge blocks may be incomplete.

$getActualBlockSize(band, xblockoff, yblockoff) Returns an integer vector of length two (xvalid, yvalid) containing the actual block size for a given block offset in band. Handles partial blocks at the edges of the raster and returns the true number of pixels. xblockoff is an integer scalar, the horizontal block offset for which to calculate the number of valid pixels, with zero indicating the left most block, 1 the next block, etc. yblockoff is likewise the vertical block offset, with zero indicating the top most block, 1 the next block, etc.

$getOverviewCount(band) Returns the number of overview layers (a.k.a. pyramids) available for band.

$buildOverviews(resampling, levels, bands) Build one or more raster overview images using the specified downsampling algorithm. resampling is a character string, one of AVERAGE, AVERAGE_MAGPHASE, RMS, BILINEAR, CUBIC, CUBICSPLINE, GAUSS, LANCZOS, MODE, NEAREST or NONE. levels is an integer vector giving the list of overview decimation factors to build (e.g., c(2, 4, 8)), or 0 to delete all overviews (at least for external overviews (.ovr) and GTiff internal overviews). bands is an integer vector giving a list of band numbers to build overviews for, or 0 to build for all bands. Note that for GTiff, overviews will be created internally if the dataset is open in update mode, while external overviews (.ovr) will be created if the dataset is open read-only. External overviews created in GTiff format may be compressed using the COMPRESS_OVERVIEW configuration option. All compression methods supported by the GTiff driver are available (e.g., set_config_option("COMPRESS_OVERVIEW", "LZW")). Since GDAL 3.6, COMPRESS_OVERVIEW is honoured when creating internal overviews of GTiff files. The GDAL documentation for gdaladdo command-line utility describes additional configuration for overview building. See also set_config_option(). No return value, called for side effects.

$getDataTypeName(band) Returns the name of the pixel data type for band. The possible data types are:

Some raster formats including GeoTIFF ("GTiff") and Erdas Imagine .img ("HFA") support sub-byte data types. Rasters can be created with these data types by specifying the "NBITS=n" creation option where n=1...7 for GTiff or n=1/2/4 for HFA. In these cases, $getDataTypeName() reports the apparent type "Byte". GTiff also supports n=9...15 (UInt16 type) and n=17...31 (UInt32 type), and n=16 is accepted for Float32 to generate half-precision floating point values.

$getNoDataValue(band) Returns the nodata value for band if one exists. This is generally a special value defined to mark pixels that are not valid data. NA is returned if a nodata value is not defined for band. Not all raster formats support a designated nodata value.

$setNoDataValue(band, nodata_value) Sets the nodata value for band. nodata_value is a numeric value to be defined as the nodata marker. Depending on the format, changing the nodata value may or may not have an effect on the pixel values of a raster that has just been created (often not). It is thus advised to call $fillRaster() explicitly if the intent is to initialize the raster to the nodata value. In any case, changing an existing nodata value, when one already exists on an initialized dataset, has no effect on the pixels whose values matched the previous nodata value. Returns logical TRUE on success or FALSE if the nodata value could not be set.

$deleteNoDataValue(band) Removes the nodata value for band. This affects only the definition of the nodata value for raster formats that support one (does not modify pixel values). No return value. An error is raised if the nodata value cannot be removed.

$getUnitType(band) Returns the name of the unit type of the pixel values for band (e.g., "m" or "ft"). An empty string "" is returned if no units are available.

$setUnitType(band, unit_type) Sets the name of the unit type of the pixel values for band. unit_type should be one of empty string "" (the default indicating it is unknown), "m" indicating meters, or "ft" indicating feet, though other nonstandard values are allowed. Returns logical TRUE on success or FALSE if the unit type could not be set.

$getScale(band) Returns the pixel value scale (units value = (raw value * scale) + offset) for band. This value (in combination with the $getOffset() value) can be used to transform raw pixel values into the units returned by $getUnitType(). Returns NA if a scale value is not defined for this band.

$setScale(band, scale) Sets the pixel value scale (units value = (raw value * scale) + offset) for band. Many raster formats do not implement this method. Returns logical TRUE on success or FALSE if the scale could not be set.

$getOffset(band) Returns the pixel value offset (units value = (raw value * scale) + offset) for band. This value (in combination with the $getScale() value) can be used to transform raw pixel values into the units returned by $getUnitType(). Returns NA if an offset value is not defined for this band.

$setOffset(band, offset) Sets the pixel value offset (units value = (raw value * scale) + offset) for band. Many raster formats do not implement this method. Returns logical TRUE on success or FALSE if the offset could not be set.

$getRasterColorInterp(band) Returns a string describing the color interpretation for band. The color interpretation values and their meanings are:

$setRasterColorInterp(band, col_interp) Sets the color interpretation for band. See above for the list of valid values for col_interp (passed as a string).

$getMinMax(band, approx_ok) Returns a numeric vector of length two containing the min/max values for band. If approx_ok is TRUE and the raster format knows these values intrinsically then those values will be returned. If that doesn't work, a subsample of blocks will be read to get an approximate min/max. If the band has a nodata value it will be excluded from the minimum and maximum. If approx_ok is FALSE, then all pixels will be read and used to compute an exact range.

$getStatistics(band, approx_ok, force) Returns a numeric vector of length four containing the minimum, maximum, mean and standard deviation of pixel values in band (excluding nodata pixels). Some raster formats will cache statistics allowing fast retrieval after the first request.

approx_ok:

• TRUE: Approximate statistics are sufficient, in which case overviews or a subset of raster tiles may be used in computing the statistics.

• FALSE: All pixels will be read and used to compute statistics (if computation is forced).

force:

• TRUE: The raster will be scanned to compute statistics. Once computed, statistics will generally be “set” back on the raster band if the format supports caching statistics. (Note: ComputeStatistics() in the GDAL API is called automatically here. This is a change in the behavior of GetStatistics() in the API, to a definitive force.)

• FALSE: Results will only be returned if it can be done quickly (i.e., without scanning the raster, typically by using pre-existing STATISTICS_xxx metadata items). NAs will be returned if statistics cannot be obtained quickly.

$clearStatistics() Clear statistics. Only implemented for now in PAM supported datasets (Persistable Auxiliary Metadata via .aux.xml file). GDAL >= 3.2.

$getHistogram(band, min, max, num_buckets, incl_out_of_range, approx_ok)
Computes raster histogram for band. min is the lower bound of the histogram. max is the upper bound of the histogram. num_buckets is the number of buckets to use (bucket size is (max - min) / num_buckets). incl_out_of_range is a logical scalar: if TRUE values below the histogram range will be mapped into the first bucket and values above will be mapped into the last bucket, if FALSE out of range values are discarded. approx_ok is a logical scalar: TRUE if an approximate histogram is OK (generally faster), or FALSE for an exactly computed histogram. Returns the histogram as a numeric vector of length num_buckets.

$getDefaultHistogram(band, force) Returns a default raster histogram for band. In the GDAL API, this method is overridden by derived classes (such as GDALPamRasterBand, VRTDataset, HFADataset...) that may be able to fetch efficiently an already stored histogram. force is a logical scalar: TRUE to force the computation of a default histogram; or if FALSE and no default histogram is available, a warning is emitted and the returned list has a 0-length histogram vector. Returns a list of length four containing named elements ⁠$min⁠ (lower bound), ⁠$max⁠ (upper bound), ⁠$num_buckets⁠ (number of buckets), and ⁠$histogram⁠ (a numeric vector of length num_buckets).

$getMetadata(band, domain) Returns a character vector of all metadata name=value pairs that exist in the specified domain, or "" (empty string) if there are no metadata items in domain (metadata in the context of the GDAL Raster Data Model: https://gdal.org/user/raster_data_model.html). Set band = 0 to retrieve dataset-level metadata, or to an integer band number to retrieve band-level metadata. Set domain = "" (empty string) to retrieve metadata in the default domain.

$getMetadataItem(band, mdi_name, domain) Returns the value of a specific metadata item named mdi_name in the specified domain, or "" (empty string) if no matching item is found. Set band = 0 to retrieve dataset-level metadata, or to an integer band number to retrieve band-level metadata. Set domain = "" (empty string) to retrieve an item in the default domain.

$setMetadataItem(band, mdi_name, mdi_value, domain) Sets the value (mdi_value) of a specific metadata item named mdi_name in the specified domain. Set band = 0 to set dataset-level metadata, or to an integer band number to set band-level metadata. Set domain = "" (empty string) to set an item in the default domain.

$getMetadataDomainList(band) Returns a character vector of metadata domains or "" (empty string). Set band = 0 to retrieve dataset-level domains, or to an integer band number to retrieve band-level domains.

$read(band, xoff, yoff, xsize, ysize, out_xsize, out_ysize) Reads a region of raster data from band. The method takes care of pixel decimation / replication if the output size (out_xsize * out_ysize) is different than the size of the region being accessed (xsize * ysize). xoff is the pixel (column) offset to the top left corner of the region of the band to be accessed (zero to start from the left side). yoff is the line (row) offset to the top left corner of the region of the band to be accessed (zero to start from the top). Note that raster row/column offsets use 0-based indexing. xsize is the width in pixels of the region to be accessed. ysize is the height in pixels of the region to be accessed. out_xsize is the width of the output array into which the desired region will be read (typically the same value as xsize). out_ysize is the height of the output array into which the desired region will be read (typically the same value as ysize). Returns a numeric or complex vector containing the values that were read. It is organized in left to right, top to bottom pixel order. NA will be returned in place of the nodata value if the raster dataset has a nodata value defined for this band. Data are read as R integer type when possible for the raster data type (Byte, Int8, Int16, UInt16, Int32), otherwise as type double (UInt32, Float32, Float64). No rescaling of the data is performed (see $getScale() and $getOffset() above). An error is raised if the read operation fails. See also the setting ⁠$readByteAsRaw⁠ above.

$write(band, xoff, yoff, xsize, ysize, rasterData) Writes a region of raster data to band. xoff is the pixel (column) offset to the top left corner of the region of the band to be accessed (zero to start from the left side). yoff is the line (row) offset to the top left corner of the region of the band to be accessed (zero to start from the top). Note that raster row/column offsets use 0-based indexing. xsize is the width in pixels of the region to write. ysize is the height in pixels of the region to write. rasterData is a numeric or complex vector containing values to write. It is organized in left to right, top to bottom pixel order. NA in rasterData should be replaced with a suitable nodata value prior to writing (see $getNoDataValue() and $setNoDataValue() above). An error is raised if the operation fails (no return value).

$getColorTable(band) Returns the color table associated with band, or NULL if there is no associated color table. The color table is returned as an integer matrix with five columns. To associate a color with a raster pixel, the pixel value is used as a subscript into the color table. This means that the colors are always applied starting at zero and ascending (see GDAL Color Table). Column 1 contains the pixel values. Interpretation of columns 2:5 depends on the value of ⁠$getPaletteInterp()⁠ (see below). For "RGB", columns 2:5 contain red, green, blue, alpha as 0-255 integer values.

$getPaletteInterp(band) If band has an associated color table, this method returns a character string with the palette interpretation for columns 2:5 of the table. An empty string ("") is returned if band does not have an associated color table. The palette interpretation values and their meanings are:

$setColorTable(band, col_tbl, palette_interp) Sets the raster color table for band (see GDAL Color Table). col_tbl is an integer matrix or data frame with either four or five columns (see $getColorTable() above). Column 1 contains the pixel values. Valid values are integers 0 and larger (note that GTiff format supports color tables only for Byte and UInt16 bands). Negative values will be skipped with a warning emitted. Interpretation of columns 2:5 depends on the value of ⁠$getPaletteInterp()⁠ (see above). For RGB, columns 2:4 contain red, green, blue as 0-255 integer values, and an optional column 5 contains alpha transparency values (defaults to 255 opaque). palette_interp is a string, one of Gray, RGB, CMYK or HLS (see $getPaletteInterp() above). Returns logical TRUE on success or FALSE if the color table could not be set.

$getDefaultRAT(band) Returns the Raster Attribute Table for band as a data frame, or NULL if there is no associated Raster Attribute Table. See the stand-alone function buildRAT() for details of the Raster Attribute Table format.

$setDefaultRAT(band, df) Sets a default Raster Attribute Table for band from data frame df. The input data frame will be checked for attribute "GDALRATTableType" which can have values of "thematic" or "athematic" (for continuous data). Columns of the data frame will be checked for attribute "GFU" (for "GDAL field usage"). If the "GFU" attribute is missing, a value of "Generic" will be used (corresponding to GFU_Generic in the GDAL API, for general purpose field). Columns with other, specific field usage values should generally be present in df, such as fields containing the set of unique (discrete) pixel values (GFU "MinMax"), pixel counts (GFU "PixelCount"), class names (GFU "Name"), color values (GFUs "Red", "⁠Green"⁠, "Blue"), etc. The data frame will also be checked for attributes "Row0Min" and "BinSize" which can have numeric values that describe linear binning. See the stand-alone function buildRAT() for details of the GDAL Raster Attribute Table format and its representation as data frame.

$flushCache() Flush all write cached data to disk. Any raster data written via GDAL calls, but buffered internally will be written to disk. Using this method does not preclude calling ⁠$close()⁠ to properly close the dataset and ensure that important data not addressed by ⁠$flushCache()⁠ is written in the file (see also ⁠$open()⁠ above). No return value, called for side effect.

$getChecksum(band, xoff, yoff, xsize, ysize) Returns a 16-bit integer (0-65535) checksum from a region of raster data on band. Floating point data are converted to 32-bit integer so decimal portions of such raster data will not affect the checksum. Real and imaginary components of complex bands influence the result. xoff is the pixel (column) offset of the window to read. yoff is the line (row) offset of the window to read. Raster row/column offsets use 0-based indexing. xsize is the width in pixels of the window to read. ysize is the height in pixels of the window to read.

$close() Closes the GDAL dataset (no return value, called for side effects). Calling $close() results in proper cleanup, and flushing of any pending writes. Forgetting to close a dataset opened in update mode on some formats such as GTiff could result in being unable to open it afterwards. The GDALRaster object is still available after calling $close(). The dataset can be re-opened on the existing filename with $open(read_only=TRUE) or $open(read_only=FALSE).

Note

If a dataset object is opened with update access (read_only = FALSE), it is not recommended to open a new dataset on the same underlying filename.

Datasets are opened in shared mode by default. This allows the sharing of GDALDataset handles for a dataset with other callers that open shared on the same filename, if the dataset is opened from the same thread. Functions in gdalraster that do processing will open input datasets in shared mode. This provides potential efficiency for cases when an object of class GDALRaster is already open in read-only mode on the same filename (avoids overhead associated with initial dataset opening by using the existing handle, and potentially makes use of existing data in the GDAL block cache). Opening in shared mode can be disabled by specifying the optional shared parameter in the class constructor.

The ⁠$read()⁠ method will perform automatic resampling if the specified output size (out_xsize * out_ysize) is different than the size of the region being read (xsize * ysize). In that case, the GDAL_RASTERIO_RESAMPLING configuration option could also be set to override the default resampling to one of BILINEAR, CUBIC, CUBICSPLINE, LANCZOS, AVERAGE or MODE (see set_config_option()).

See Also

Package overview in help("gdalraster-package")

vignette("raster-api-tutorial")

read_ds() convenience wrapper for GDALRaster$read()

Examples

lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds <- new(GDALRaster, lcp_file)

## print information about the dataset to the console
ds$info()

## retrieve the raster format name
ds$getDriverShortName()
ds$getDriverLongName()

## retrieve a list of files composing the dataset
ds$getFileList()

## retrieve dataset parameters
ds$getRasterXSize()
ds$getRasterYSize()
ds$getGeoTransform()
ds$getProjection()
ds$getRasterCount()
ds$bbox()
ds$res()
ds$dim()

## retrieve some band-level parameters
ds$getDescription(band = 1)
ds$getBlockSize(band = 1)
ds$getOverviewCount(band = 1)
ds$getDataTypeName(band = 1)
# LCP format does not support an intrinsic nodata value so this returns NA:
ds$getNoDataValue(band = 1)

## LCP driver reports several dataset- and band-level metadata
## see the format description at https://gdal.org/drivers/raster/lcp.html
## set band = 0 to retrieve dataset-level metadata
## set domain = "" (empty string) for the default metadata domain
ds$getMetadata(band = 0, domain = "")

## retrieve metadata for a band as a vector of name=value pairs
ds$getMetadata(band = 4, domain = "")

## retrieve the value of a specific metadata item
ds$getMetadataItem(band = 2, mdi_name = "SLOPE_UNIT_NAME", domain = "")

## read one row of pixel values from band 1 (elevation)
## raster row/column index are 0-based
## the upper left corner is the origin
## read the tenth row:
ncols <- ds$getRasterXSize()
rowdata <- ds$read(band = 1, xoff = 0, yoff = 9,
                   xsize = ncols, ysize = 1,
                   out_xsize = ncols, out_ysize = 1)
head(rowdata)

ds$close()

## create a new raster using lcp_file as a template
new_file <- file.path(tempdir(), "storml_newdata.tif")
rasterFromRaster(srcfile = lcp_file,
                 dstfile = new_file,
                 nbands = 1,
                 dtName = "Byte",
                 init = -9999)

ds_new <- new(GDALRaster, new_file, read_only = FALSE)

## write random values to all pixels
set.seed(42)
ncols <- ds_new$getRasterXSize()
nrows <- ds_new$getRasterYSize()
for (row in 0:(nrows - 1)) {
    rowdata <- round(runif(ncols, 0, 100))
    ds_new$write(band = 1,
                 xoff = 0,
                 yoff = row,
                 xsize = ncols,
                 ysize = 1,
                 rowdata)
}

## re-open in read-only mode when done writing
## this will ensure flushing of any pending writes (implicit $close)
ds_new$open(read_only = TRUE)

## getStatistics returns min, max, mean, sd, and sets stats in the metadata
ds_new$getStatistics(band = 1, approx_ok = FALSE, force = TRUE)
ds_new$getMetadataItem(band = 1, "STATISTICS_MEAN", "")

## close the dataset for proper cleanup
ds_new$close()
deleteDataset(new_file)


## using a GDAL Virtual File System handler '/vsicurl/'
## see: https://gdal.org/user/virtual_file_systems.html
url <- "/vsicurl/https://raw.githubusercontent.com/"
url <- paste0(url, "usdaforestservice/gdalraster/main/sample-data/")
url <- paste0(url, "lf_elev_220_mt_hood_utm.tif")

set_config_option("GDAL_HTTP_CONNECTTIMEOUT", "20")
set_config_option("GDAL_HTTP_TIMEOUT", "20")
if (http_enabled() && vsi_stat(url)) {
  ds <- new(GDALRaster, url)
  plot_raster(ds, legend = TRUE, main = "Mount Hood elevation (m)")
  ds$close()
}
set_config_option("GDAL_HTTP_CONNECTTIMEOUT", "")
set_config_option("GDAL_HTTP_TIMEOUT", "")

Description

geos_version() returns version information for the GEOS library in use by GDAL. Requires GDAL >= 3.4.

Usage

geos_version()

Value

A list of length four containing:

• name - a string formatted as "major.minor.patch"

• major - major version as integer

• minor - minor version as integer

• patch - patch version as integer

List elements will be NA if GDAL < 3.4.

See Also

gdal_version(), proj_version()

Examples

geos_version()

Description

get_cache_used() returns the amount of memory currently in use for GDAL block caching. This a wrapper for GDALGetCacheUsed64() with return value as MB.

Usage

get_cache_used()

Value

Integer. Amount of cache memory in use in MB.

See Also

GDAL Block Cache

Examples

get_cache_used()

Description

get_config_option() gets the value of GDAL runtime configuration option. Configuration options are essentially global variables the user can set. They are used to alter the default behavior of certain raster format drivers, and in some cases the GDAL core. For a full description and listing of available options see https://gdal.org/user/configoptions.html.

Usage

get_config_option(key)

Arguments

Value

Character. The value of a (key, value) option previously set with set_config_option(). An empty string ("") is returned if key is not found.

See Also

set_config_option()

vignette("gdal-config-quick-ref")

Examples

## this option is set during initialization of the gdalraster package
get_config_option("OGR_CT_FORCE_TRADITIONAL_GIS_ORDER")

Description

get_num_cpus() returns the number of processors detected by GDAL. Wrapper of CPLGetNumCPUs() in the GDAL Common Portability Library.

Usage

get_num_cpus()

Value

Integer scalar, number of CPUs.

Examples

get_num_cpus()

Description

get_pixel_line() converts geospatial coordinates to pixel/line (raster column, row numbers). The upper left corner pixel is the raster origin (0,0) with column, row increasing left to right, top to bottom.

Usage

get_pixel_line(xy, gt)

Arguments

Value

Integer matrix of raster pixel/line.

Note

This function applies the inverse geotransform to the input points. If gt is given as the numeric vector, no bounds checking is done (i.e., min pixel/line could be less than zero and max pixel/line could be greater than the raster x/y size). If gt is obtained from an object of class GDALRaster, then NA is returned for points that fall outside the raster extent and a warning emitted giving the number points that were outside. This latter case is equivalent to calling the ⁠$get_pixel_line()⁠ class method on the GDALRaster object (see Examples).

See Also

GDALRaster$getGeoTransform(), inv_geotransform()

Examples

pt_file <- system.file("extdata/storml_pts.csv", package="gdalraster")
# id, x, y in NAD83 / UTM zone 12N
pts <- read.csv(pt_file)
print(pts)

raster_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds <- new(GDALRaster, raster_file)
gt <- ds$getGeoTransform()
get_pixel_line(pts[, -1], gt)

# or, using the class method
ds$get_pixel_line(pts[, -1])

# add a point outside the raster extent
pts[11, ] <- c(11, 323318, 5105104)
get_pixel_line(pts[, -1], gt)

# with bounds checking on the raster extent
ds$get_pixel_line(pts[, -1])

ds$close()

Description

get_usable_physical_ram() returns the total physical RAM, usable by a process, in bytes. It will limit to 2 GB for 32 bit processes. Starting with GDAL 2.4.0, it will also take into account resource limits (virtual memory) on Posix systems. Starting with GDAL 3.6.1, it will also take into account RLIMIT_RSS on Linux. Wrapper of CPLGetUsablePhysicalRAM() in the GDAL Common Portability Library.

Usage

get_usable_physical_ram()

Value

Numeric scalar, number of bytes as bit64::integer64 type (or 0 in case of failure).

Note

This memory may already be partly used by other processes.

Examples

get_usable_physical_ram()

Description

getCreationOptions() returns the list of creation options supported by a GDAL format driver as an XML string (invisibly). Wrapper for GDALGetDriverCreationOptionList() in the GDAL API. Information about the available creation options is also printed to the console by default.

Usage

getCreationOptions(format, filter = NULL)

Arguments

Value

Invisibly, an XML string that describes the full list of creation options or empty string "" (full output of GDALGetDriverCreationOptionList() in the GDAL API).

See Also

GDALRaster-class, create(), createCopy()

Examples

getCreationOptions("GTiff", filter="COMPRESS")

Description

has_geos() returns a logical value indicating whether GDAL was built against the GEOS library. GDAL built with GEOS is a system requirement as of gdalraster 1.10.0, so this function will always return TRUE (may be removed in a future version).

Usage

has_geos()

Value

Logical. TRUE if GEOS is available, otherwise FALSE.

Examples

has_geos()

Description

has_spatialite() returns a logical value indicating whether GDAL was built with support for the SpatiaLite library. SpatiaLite extends the SQLite core to support full Spatial SQL capabilities.

Usage

has_spatialite()

Details

GDAL supports executing SQL statements against a datasource. For most file formats (e.g. Shapefiles, GeoJSON, FlatGeobuf files), the built-in OGR SQL dialect will be used by default. It is also possible to request the alternate "SQLite" dialect, which will use the SQLite engine to evaluate commands on GDAL datasets. This assumes that GDAL is built with support for SQLite, and preferably with Spatialite support too to benefit from spatial functions.

Value

Logical scalar. TRUE if SpatiaLite is available to GDAL.

Note

All GDAL/OGR drivers for database systems, e.g., PostgreSQL / PostGIS, Oracle Spatial, SQLite / Spatialite RDBMS, GeoPackage, etc., override the GDALDataset::ExecuteSQL() function with a dedicated implementation and, by default, pass the SQL statements directly to the underlying RDBMS. In these cases the SQL syntax varies in some particulars from OGR SQL. Also, anything possible in SQL can then be accomplished for these particular databases. For those drivers, it is also possible to explicitly request the OGRSQL or SQLite dialects, although performance will generally be much less than the native SQL engine of those database systems.

See Also

ogrinfo(), ogr_execute_sql()

OGR SQL dialect and SQLITE SQL dialect:
https://gdal.org/user/ogr_sql_sqlite_dialect.html

Examples

has_spatialite()

Description

http_enabled() returns TRUE if libcurl support is enabled. Wrapper of CPLHTTPEnabled() in the GDAL Common Portability Library.

Usage

http_enabled()

Value

Logical scalar, TRUE if GDAL was built with libcurl support.

Examples

http_enabled()

Description

inv_geotransform() inverts a vector of geotransform coefficients. This converts the equation from being:
raster pixel/line (column/row) -> geospatial x/y coordinate
to:
geospatial x/y coordinate -> raster pixel/line (column/row)

Usage

inv_geotransform(gt)

Arguments

Value

Numeric vector of length six containing the inverted geotransform. The output vector will contain NAs if the input geotransform is uninvertable.

See Also

GDALRaster$getGeoTransform(), get_pixel_line()

Examples

elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
ds <- new(GDALRaster, elev_file)
invgt <- ds$getGeoTransform() |> inv_geotransform()
ds$close()

ptX = 324181.7
ptY = 5103901.4

## for a point x, y in the spatial reference system of elev_file
## raster pixel (column number):
pixel <- floor(invgt[1] +
               invgt[2] * ptX +
               invgt[3] * ptY)

## raster line (row number):
line <- floor(invgt[4] +
              invgt[5] * ptX +
              invgt[6] * ptY)

## get_pixel_line() applies this conversion

Description

inv_project() transforms geospatial x/y coordinates to longitude/latitude in the same geographic coordinate system used by the given projected spatial reference system. The output long/lat can optionally be set to a specific geographic coordinate system by specifying a well known name (see Details).

Usage

inv_project(pts, srs, well_known_gcs = "")

Arguments

Details

By default, the geographic coordinate system of the projection specified by srs will be used. If a specific geographic coordinate system is desired, then well_known_gcs can be set to one of the values below:

The returned array will always be in longitude, latitude order (traditional GIS order) regardless of the axis order defined for the names above.

Value

Numeric array of longitude, latitude. An error is raised if the transformation cannot be performed.

See Also

transform_xy()

Examples

pt_file <- system.file("extdata/storml_pts.csv", package="gdalraster")
## id, x, y in NAD83 / UTM zone 12N
pts <- read.csv(pt_file)
print(pts)
inv_project(pts[,-1], epsg_to_wkt(26912))
inv_project(pts[,-1], epsg_to_wkt(26912), "NAD27")

Description

This topic contains documentation and helper functions for defining an OGR feature class. ogr_def_field() creates an attribute field definition, a list containing the field data type and potentially other optional field properties. ogr_def_geom_field() similarly creates a geometry field definition. A list containing zero or more attribute field definitions, along with one or more geometry field definitions, comprise an OGR feature class definition (a.k.a. layer definition). ogr_def_layer() initializes such a list with a geometry field. Attribute fields can be added to a feature class definition with calls to ogr_def_field() as in the examples.

Usage

ogr_def_field(
  fld_type,
  fld_subtype = NULL,
  fld_width = NULL,
  fld_precision = NULL,
  is_nullable = NULL,
  is_unique = NULL,
  is_ignored = NULL,
  default_value = NULL
)

ogr_def_geom_field(
  geom_type,
  srs = NULL,
  is_nullable = NULL,
  is_ignored = NULL
)

ogr_def_layer(geom_type, geom_fld_name = "geom", srs = NULL)

Arguments

Details

All features in an OGR Layer share a common schema (feature class), modeled in GDAL as OGR Feature Definition. The feature class definition includes the set of attribute fields and their data types and the geometry field(s). In R, a feature class definition is represented as a list, having as names the attribute/geometry field names, with each list element holding a field definition.

An attribute field definition is a list with named elements:

$type       : OGR Field Type ("OFTReal", "OFTString" etc.)
$subtype    : optional ("OFSTBoolean", ...)
$width      : optional max number of characters
$precision  : optional number of digits after the decimal point
$is_nullable: optional NOT NULL constraint (logical scalar)
$is_unique  : optional UNIQUE constraint (logical scalar)
$default    : optional default value as character string
$is_ignored : optionally ignored when retrieving features (logical scalar)
$is_geom    : FALSE (the default) for attribute fields

An OGR field type is specified as a character string with possible values: OFTInteger, OFTIntegerList, OFTReal, OFTRealList, OFTString, OFTStringList, OFTBinary, OFTDate, OFTTime, OFTDateTime, OFTInteger64, OFTInteger64List.

An optional field subtype is specified as a character string with possible values: OFSTNone, OFSTBoolean, OFSTInt16, OFSTFloat32, OFSTJSON, OFSTUUID.

By default, fields are nullable, have no unique constraint, and are not ignored (i.e., not omitted when fetching features). Not-null and unique constraints are not supported by all format drivers.

A default field value is taken into account by format drivers (generally those with a SQL interface) that support it at field creation time. If given in the field definition, ⁠$default⁠ must be a character string. The accepted values are "NULL", a numeric value (e.g., "0"), a literal value enclosed between single quote characters (e.g., "'a default value'", with any inner single quote characters escaped by repetition of the single quote character), "CURRENT_TIMESTAMP", "CURRENT_TIME", "CURRENT_DATE" or a driver-specific expression (that might be ignored by other drivers). For a datetime literal value, format should be "'YYYY/MM/DD HH:MM:SS[.sss]'" (considered as UTC time).

A geometry field definition is a list with named elements:

$type       : geom type ("Point", "Polygon", etc.)
$srs        : optional spatial reference as WKT string
$is_nullable: optional NOT NULL constraint (logical scalar)
$is_ignored : optionally ignored when retrieving features (logical scalar)
$is_geom    : TRUE (required) for geometry fields

Typically, there is one geometry field on a layer, but some formats support more than one geometry column per table (e.g., PostGIS).

Geometry types are specified as a character string containing OGC WKT. Common types include: Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon. See the GDAL documentation for a list of all supported geometry types:
https://gdal.org/api/vector_c_api.html#_CPPv418OGRwkbGeometryType

Format drivers may or may not support not-null constraints on attribute and geometry fields. If they support creating fields with not-null constraints, this is generally before creating any features to the layer. In some cases, a not-null constraint may be available as a layer creation option. For example, GeoPackage format has a layer creation option ⁠GEOMETRY_NULLABLE=[YES/NO]⁠.

Note

The feature id (FID) is a special property of a feature and not treated as an attribute of the feature. Additional information is given in the GDAL documentation for the OGR SQL and SQLite SQL dialects. Implications for SQL statements and result sets may depend on the dialect used.

Some vector formats do not support schema definition prior to creating features. For example, with GeoJSON only the Feature object has a member with name properties. The specification does not require all Feature objects in a collection to have the same schema of properties, nor does it require all Feature objects in a collection to have geometry of the same type (https://geojson.org/).

See Also

ogr_ds_create(), ogr_layer_create(), ogr_field_create(), ogrinfo()

WKT representation of geometry:
https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry

Examples

dsn <- file.path(tempdir(), "test.sqlite")
opt <- NULL
if (has_spatialite())
  opt <- "SPATIALITE=YES"
ogr_ds_create("SQLite", dsn, dsco = opt)

# define a layer
defn <- ogr_def_layer("Point", srs = epsg_to_wkt(4326))
defn$fld1_int64 <- ogr_def_field("OFTInteger64")
defn$fld2_string <- ogr_def_field("OFTString")

if (ogr_ds_test_cap(dsn)$CreateLayer)
  ogr_layer_create(dsn, "layer1", layer_defn = defn)

ogr_ds_layer_names(dsn)
ogr_layer_field_names(dsn, "layer1")

deleteDataset(dsn)

Description

This set of functions can be used to create new vector datasets, test existence of dataset/layer/field, test dataset and layer capabilities, create new layers in an existing dataset, delete layers, create new attribute and geometry fields on an existing layer, rename and delete fields, and edit data with SQL statements.

Usage

ogr_ds_exists(dsn, with_update = FALSE)

ogr_ds_format(dsn)

ogr_ds_test_cap(dsn, with_update = TRUE)

ogr_ds_create(
  format,
  dsn,
  layer = NULL,
  layer_defn = NULL,
  geom_type = NULL,
  srs = NULL,
  fld_name = NULL,
  fld_type = NULL,
  dsco = NULL,
  lco = NULL,
  overwrite = FALSE
)

ogr_ds_layer_count(dsn)

ogr_ds_layer_names(dsn)

ogr_layer_exists(dsn, layer)

ogr_layer_test_cap(dsn, layer, with_update = TRUE)

ogr_layer_create(
  dsn,
  layer,
  layer_defn = NULL,
  geom_type = NULL,
  srs = NULL,
  lco = NULL
)

ogr_layer_field_names(dsn, layer)

ogr_layer_delete(dsn, layer)

ogr_field_index(dsn, layer, fld_name)

ogr_field_create(
  dsn,
  layer,
  fld_name,
  fld_defn = NULL,
  fld_type = "OFTInteger",
  fld_subtype = "OFSTNone",
  fld_width = 0L,
  fld_precision = 0L,
  is_nullable = TRUE,
  is_ignored = FALSE,
  is_unique = FALSE,
  default_value = ""
)

ogr_geom_field_create(
  dsn,
  layer,
  fld_name,
  geom_fld_defn = NULL,
  geom_type = NULL,
  srs = NULL,
  is_nullable = TRUE,
  is_ignored = FALSE
)

ogr_field_rename(dsn, layer, fld_name, new_name)

ogr_field_delete(dsn, layer, fld_name)

ogr_execute_sql(dsn, sql, spatial_filter = NULL, dialect = NULL)

Arguments

Details

These functions are complementary to ogrinfo() and ogr2ogr() for vector data management. They are also intended to support vector I/O in a future release of gdalraster. Bindings to OGR wrap portions of the GDAL Vector API (ogr_core.h and ogr_api.h, https://gdal.org/api/vector_c_api.html).

ogr_ds_exists() tests whether a vector dataset can be opened from the given data source name (DSN), potentially testing for update access. Returns a logical scalar.

ogr_ds_format() returns a character string containing the short name of the format driver for a given DSN, or NULL if the dataset cannot be opened as a vector source.

ogr_ds_test_cap() tests the capabilities of a vector data source, attempting to open it with update access by default. Returns a list of capabilities with values TRUE or FALSE, or NULL is returned if dsn cannot be opened with the requested access. Wrapper of GDALDatasetTestCapability() in the GDAL API. The returned list contains the following named elements:

• CreateLayer: TRUE if this datasource can create new layers

• DeleteLayer: TRUE if this datasource can delete existing layers

• CreateGeomFieldAfterCreateLayer: TRUE if the layers of this datasource support geometry field creation just after layer creation

• CurveGeometries: TRUE if this datasource supports curve geometries

• Transactions: TRUE if this datasource supports (efficient) transactions

• EmulatedTransactions: TRUE if this datasource supports transactions through emulation

• RandomLayerRead: TRUE if this datasource has a dedicated GetNextFeature() implementation, potentially returning features from layers in a non-sequential way

• RandomLayerWrite: TRUE if this datasource supports calling CreateFeature() on layers in a non-sequential way

ogr_ds_create() creates a new vector datasource, optionally also creating a layer, and optionally creating one or more fields on the layer. The attribute fields and geometry field(s) to create can be specified as a feature class definition (layer_defn as list, see ogr_define), or alternatively, by giving the geom_type and srs, optionally along with one fld_name and fld_type to be created in the layer. Returns a logical scalar, TRUE indicating success.

ogr_ds_layer_count() returns the number of layers in a vector dataset.

ogr_ds_layer_names() returns a character vector of layer names in a vector dataset, or NULL if no layers are found.

ogr_layer_exists() tests whether a layer can be accessed by name in a given vector dataset. Returns a logical scalar.

ogr_layer_test_cap() tests whether a layer supports named capabilities, attempting to open the dataset with update access by default. Returns a list of capabilities with values TRUE or FALSE. NULL is returned if dsn cannot be opened with the requested access, or layer cannot be found. The returned list contains the following named elements: RandomRead, SequentialWrite, RandomWrite, UpsertFeature, FastSpatialFilter, FastFeatureCount, FastGetExtent, FastSetNextByIndex, CreateField, CreateGeomField, DeleteField, ReorderFields, AlterFieldDefn, AlterGeomFieldDefn, DeleteFeature, StringsAsUTF8, Transactions, CurveGeometries. See the GDAL documentation for OGR_L_TestCapability().

ogr_layer_create() creates a new layer in an existing vector data source, with a specified geometry type and spatial reference definition. This function also accepts a feature class definition given as a list of field names and their definitions (see ogr_define). (Note: use ogr_ds_create() to create single-layer formats such as "ESRI Shapefile", "FlatGeobuf", "GeoJSON", etc.) Returns a logical scalar, TRUE indicating success.

ogr_layer_field_names() returns a character vector of field names on a layer, or NULL if no fields are found.

ogr_layer_delete() deletes an existing layer in a vector dataset. Returns a logical scalar, TRUE indicating success.

ogr_field_index() tests for existence of an attribute field by name. Returns the field index on the layer (0-based), or -1 if the field does not exist.

ogr_field_create() creates a new attribute field of specified data type in a given DSN/layer. Several optional field properties can be specified in addition to the type. Returns a logical scalar, TRUE indicating success.

ogr_geom_field_create() creates a new geometry field of specified type in a given DSN/layer. Returns a logical scalar, TRUE indicating success.

ogr_field_rename() renames an existing field on a vector layer. Not all format drivers support this function. Some drivers may only support renaming a field while there are still no features in the layer. AlterFieldDefn is the relevant layer capability to check. Returns a logical scalar, TRUE indicating success.

ogr_field_delete() deletes an existing field on a vector layer. Not all format drivers support this function. Some drivers may only support deleting a field while there are still no features in the layer. Returns a logical scalar, TRUE indicating success.

ogr_execute_sql() executes an SQL statement against the data store. This function can be used to modify the schema or edit data using SQL (e.g., ⁠ALTER TABLE⁠, ⁠DROP TABLE⁠, ⁠CREATE INDEX⁠, ⁠DROP INDEX⁠, INSERT, UPDATE, DELETE). Currently, this function does not return a result set for a SELECT statement. Returns NULL invisibly. Wrapper of GDALDatasetExecuteSQL() in the GDAL C API.

Note

The OGR SQL document linked under See Also contains information on the SQL dialect supported internally by GDAL/OGR. Some format drivers (e.g., PostGIS) pass the SQL directly through to the underlying RDBMS (unless OGRSQL is explicitly passed as the dialect). The SQLite dialect can also be requested with the SQLite string passed as the dialect argument of ogr_execute_sql(). This assumes that GDAL/OGR is built with support for SQLite, and preferably also with Spatialite support to benefit from spatial functions. The GDAL document for SQLite dialect has detailed information.

Other SQL dialects may also be present for some vector formats. For example, the "INDIRECT_SQLITE" dialect might potentially be used with GeoPackage format (https://gdal.org/drivers/vector/gpkg.html#sql).

The function ogrinfo() can also be used to edit data with SQL statements (GDAL >= 3.7).

The name of the geometry column of a layer is empty ("") with some formats such as ESRI Shapefile and FlatGeobuf. Implications for SQL may depend on the dialect used. See the GDAL documentation for the "OGR SQL" and "SQLite" dialects for details.

See Also

gdal_formats(), has_spatialite(), ogr_def_field(), ogr_def_layer(), ogrinfo(), ogr2ogr()

OGR SQL dialect and SQLite SQL dialect:
https://gdal.org/user/ogr_sql_sqlite_dialect.html

Examples

dsn <- file.path(tempdir(), "test1.gpkg")
ogr_ds_create("GPKG", dsn)
ogr_ds_exists(dsn, with_update = TRUE)
ogr_ds_layer_count(dsn)
ogr_ds_test_cap(dsn)
ogr_layer_exists(dsn, "layer1")
if (ogr_ds_test_cap(dsn)$CreateLayer) {
  opt <- c("GEOMETRY_NULLABLE=NO", "DESCRIPTION=test layer")
  ogr_layer_create(dsn, "layer1", geom_type = "Polygon", srs = "EPSG:5070",
                   lco = opt)
}
ogr_ds_layer_count(dsn)
ogr_layer_exists(dsn, "layer1")
ogr_ds_layer_names(dsn)

ogr_layer_field_names(dsn, "layer1")
ogr_field_index(dsn, "layer1", "field1")
if (ogr_layer_test_cap(dsn, "layer1")$CreateField) {
  ogr_field_create(dsn, "layer1", "field1",
                   fld_type = "OFTInteger64",
                   is_nullable = FALSE)
  ogr_field_create(dsn, "layer1", "field2",
                   fld_type = "OFTString")
}
ogr_field_index(dsn, "layer1", "field1")
ogr_layer_field_names(dsn, "layer1")

# delete a field
if (ogr_layer_test_cap(dsn, "layer1")$DeleteField) {
  ogr_field_delete(dsn, "layer1", "field2")
}
ogr_layer_field_names(dsn, "layer1")

# define a feature class (layer definition)
defn <- ogr_def_layer("Point", srs = epsg_to_wkt(4326))
# add the attribute fields
defn$fld1_name <- ogr_def_field("OFTInteger64",
                                is_nullable = FALSE,
                                is_unique = TRUE)
defn$fld2_name <- ogr_def_field("OFTString",
                                fld_width = 25,
                                is_nullable = FALSE,
                                default_value = "'a default string'")
defn$third_field <- ogr_def_field("OFTReal",
                                  default_value = "0.0")

ogr_layer_create(dsn, "layer2", layer_defn = defn)
ogr_ds_layer_names(dsn)
ogr_layer_field_names(dsn, "layer2")

# add a field using SQL instead
sql <- "ALTER TABLE layer2 ADD field4 float"
ogr_execute_sql(dsn, sql)
ogr_layer_field_names(dsn, "layer2")

# rename a field
if (ogr_layer_test_cap(dsn, "layer1")$AlterFieldDefn) {
  ogr_field_rename(dsn, "layer2", "field4", "renamed_field")
}
ogr_layer_field_names(dsn, "layer2")

# GDAL >= 3.7
if (as.integer(gdal_version()[2]) >= 3070000)
  ogrinfo(dsn, "layer2")

deleteDataset(dsn)

# edit data using SQL
src <- system.file("extdata/ynp_fires_1984_2022.gpkg", package="gdalraster")
perims_shp <- file.path(tempdir(), "mtbs_perims.shp")
ogr2ogr(src, perims_shp, src_layers = "mtbs_perims")
ogr_ds_format(perims_shp)
ogr_ds_layer_names(perims_shp)
ogr_layer_field_names(perims_shp, "mtbs_perims")

if (ogr_layer_test_cap(perims_shp, "mtbs_perims")$CreateField) {
  sql <- "ALTER TABLE mtbs_perims ADD burn_bnd_ha float"
  ogr_execute_sql(perims_shp, sql)
  # with GDAL >= 3.7, equivalent to:
  # ogrinfo(perims_shp, cl_arg = c("-sql", sql), read_only = FALSE)
}

sql <- "UPDATE mtbs_perims SET burn_bnd_ha = (burn_bnd_ac / 2.471)"
ogr_execute_sql(perims_shp, sql, dialect = "SQLite")
ogr_layer_field_names(perims_shp, "mtbs_perims")

# if GDAL >= 3.7:
#   ogrinfo(perims_shp, "mtbs_perims")
# or, for output incl. the feature data (omit the default "-so" arg):
#   ogrinfo(perims_shp, "mtbs_perims", cl_arg = "-nomd")

deleteDataset(perims_shp)

Description

ogr2ogr() is a wrapper of the ogr2ogr command-line utility (see https://gdal.org/programs/ogr2ogr.html). This function can be used to convert simple features data between file formats. It can also perform various operations during the process, such as spatial or attribute selection, reducing the set of attributes, setting the output coordinate system or even reprojecting the features during translation. Refer to the GDAL documentation at the URL above for a description of command-line arguments that can be passed in cl_arg.

Usage

ogr2ogr(
  src_dsn,
  dst_dsn,
  src_layers = NULL,
  cl_arg = NULL,
  open_options = NULL
)

Arguments

Value

Logical indicating success (invisible TRUE). An error is raised if the operation fails.

Note

For progress reporting, see command-line argument -progress: Display progress on terminal. Only works if input layers have the "fast feature count" capability.

See Also

ogrinfo(), the ogr_manage utilities

translate() for raster data

Examples

src <- system.file("extdata/ynp_fires_1984_2022.gpkg", package="gdalraster")

# Convert GeoPackage to Shapefile
shp_file <- file.path(tempdir(), "ynp_fires.shp")
ogr2ogr(src, shp_file, src_layers = "mtbs_perims")

# Reproject to WGS84
ynp_wgs84 <- file.path(tempdir(), "ynp_fires_wgs84.gpkg")
args <- c("-t_srs", "EPSG:4326")
ogr2ogr(src, ynp_wgs84, cl_arg = args)

# Clip to a bounding box (xmin, ymin, xmax, ymax in the source SRS)
# This will select features whose geometry intersects the bounding box.
# The geometries themselves will not be clipped unless "-clipsrc" is
# specified.
# The source SRS can be overridden with "-spat_srs" "<srs_def>"
ynp_clip <- file.path(tempdir(), "ynp_fires_aoi_clip.gpkg")
bb <- c(469685.97, 11442.45, 544069.63, 85508.15)
args <- c("-spat", bb)
ogr2ogr(src, ynp_clip, cl_arg = args)

# Filter features by a -where clause
ynp_filtered <- file.path(tempdir(), "ynp_fires_2000_2022.gpkg")
sql <- "ig_year >= 2000 ORDER BY ig_year"
args <- c("-where", sql)
ogr2ogr(src, ynp_filtered, src_layers = "mtbs_perims", cl_arg = args)

deleteDataset(shp_file)
deleteDataset(ynp_wgs84)
deleteDataset(ynp_clip)
deleteDataset(ynp_filtered)

Description

ogrinfo() is a wrapper of the ogrinfo command-line utility (see https://gdal.org/programs/ogrinfo.html). This function lists information about an OGR-supported data source. It is also possible to edit data with SQL statements. Refer to the GDAL documentation at the URL above for a description of command-line arguments that can be passed in cl_arg. Requires GDAL >= 3.7.

Usage

ogrinfo(
  dsn,
  layers = NULL,
  cl_arg = as.character(c("-so", "-nomd")),
  open_options = NULL,
  read_only = TRUE,
  cout = TRUE
)