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: autoCreateWarpedVRT(), 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_find_epsg(), srs_get_name(), 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_make_valid(), g_name(), g_overlaps(), g_summary(), g_sym_difference(), g_touches(), g_transform(), g_union(), g_within(), g_wk2wk(), geos_version()

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

• 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.

Temporary files are created in some examples which have cleanup code wrapped in dontshow{}. While the cleanup code is not shown in the documentation, note that this code runs by default if examples are run with example().

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()
  
}

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

autoCreateWarpedVRT() creates a warped virtual dataset representing the input raster warped into a target coordinate system. The output virtual dataset will be "north-up" in the target coordinate system. GDAL automatically determines the bounds and resolution of the output virtual raster which should be large enough to include all the input raster. Wrapper of GDALAutoCreateWarpedVRT() in the GDAL Warper API.

Usage

autoCreateWarpedVRT(
  src_ds,
  dst_wkt,
  resample_alg,
  src_wkt = "",
  max_err = 0,
  alpha_band = FALSE
)

Arguments

Value

An object of class GDALRaster for the new virtual dataset. An error is raised if the operation fails.

Note

The returned dataset will have no associated filename for itself. If you want to write the virtual dataset to a VRT file, use the ⁠$setFilename()⁠ method on the returned GDALRaster object to assign a filename before it is closed.

Examples

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

ds2 <- autoCreateWarpedVRT(ds, epsg_to_wkt(5070), "Bilinear")
ds2$info()

## set filename before close if a VRT file is needed for the virtual dataset
# ds2$setFilename("/path/to/file.vrt")

ds2$close()
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()

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()


# 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()

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()



## 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()



## 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)

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 (see Details)

## Constructors
cmb <- new(CmbTable, keyLen)
# or, giving the variable names:
cmb <- new(CmbTable, keyLen, varNames)

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

Details

Constructors

new(CmbTable, keyLen)
Default 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.

Methods

$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()

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,
  return_obj = FALSE
)

Arguments

Value

By default, returns a logical value indicating success (invisible TRUE, output written to dst_filename). An error is raised if the operation fails. An object of class GDALRaster open on the output dataset will be returned if return_obj = TRUE.

Note

dst_filename may be an empty string ("") with format = "MEM" and return_obj = TRUE to create an In-memory Raster (https://gdal.org/drivers/raster/mem.html).

See Also

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

Examples

new_file <- file.path(tempdir(), "newdata.tif")
ds <- create(format="GTiff",
             dst_filename = new_file,
             xsize = 143,
             ysize = 107,
             nbands = 1,
             dataType = "Int16",
             return_obj=TRUE)

# EPSG:26912 - NAD83 / UTM zone 12N
ds$setProjection(epsg_to_wkt(26912))

gt <- c(323476, 30, 0, 5105082, 0, -30)
ds$setGeoTransform(gt)

ds$setNoDataValue(band = 1, -9999)
ds$fillRaster(band = 1, -9999, 0)

# ...

# close the dataset when done
ds$close()

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()

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,
  return_obj = FALSE
)

Arguments

Value

By default, returns a logical value indicating success (invisible TRUE, output written to dst_filename). An error is raised if the operation fails. An object of class GDALRaster open on the output dataset will be returned if return_obj = TRUE.

Note

dst_filename may be an empty string ("") with format = "MEM" and return_obj = TRUE to create an In-memory Raster (https://gdal.org/drivers/raster/mem.html).

See Also

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

Examples

lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
tif_file <- file.path(tempdir(), "storml_lndscp.tif")
ds <- createCopy(format = "GTiff",
                 dst_filename = tif_file,
                 src_filename = lcp_file,
                 options = "COMPRESS=LZW",
                 return_obj = TRUE)

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()

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)

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),]

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)
  
}

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(
  geom,
  dist,
  quad_segs = 30L,
  as_wkb = TRUE,
  as_iso = FALSE,
  byte_order = "LSB",
  quiet = FALSE
)

Arguments

Value

A polygon as WKB raw vector or WKT string, or a list/character vector of polygons as WKB/WKT with length equal to length(geom). NA is returned with a warning if an input geometry cannot be converted into an OGR geometry object, or if an error occurs in the call to Buffer() in the underlying OGR API.

Examples

g_buffer("POINT (0 0)", dist = 10, as_wkb = FALSE)

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(geom, quiet = FALSE)

Arguments

Value

logical vector of the same length as the number of input geometries in geom, containing TRUE for the corresponding geometies that are empty or FALSE for non-empty geometries.

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(geom, quiet = FALSE)

Arguments

Value

logical vector of the same length as the number of input geometries in geom, containing TRUE for the corresponding geometies that are valid or FALSE for invalid geometries.

Examples

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

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_make_valid() attempts to make an invalid geometry valid without losing vertices. Already-valid geometries are cloned without further intervention. Wrapper of OGR_G_MakeValid()/OGR_G_MakeValidEx() in the GDAL Vector C API.

Usage

g_make_valid(
  geom,
  method = "LINEWORK",
  keep_collapsed = FALSE,
  as_wkb = TRUE,
  as_iso = FALSE,
  byte_order = "LSB",
  quiet = FALSE
)

Arguments

Details

LINEWORK is the default method, which combines all rings into a set of noded lines and then extracts valid polygons from that linework. The STRUCTURE method (requires GEOS >= 3.10 and GDAL >= 3.4) first makes all rings valid, then merges shells and subtracts holes from shells to generate a valid result. Assumes that holes and shells are correctly categorized.

KEEP_COLLAPSED only applies to the STRUCTURE method:

• FALSE (the default): collapses are converted to empty geometries

• TRUE: collapses are converted to a valid geometry of lower dimension

Value

A geometry as WKB raw vector or WKT string, or a list/character vector of geometries as WKB/WKT with length equal to length(geom). NA is returned with a warning if an input geometry cannot be converted into an OGR geometry object, or if an error occurs in the call to MakeValid() in the underlying OGR API.

Note

This function is built on the GEOS >= 3.8 library, check it for the definition of the geometry operation. If OGR is built without GEOS >= 3.8, this function will return a clone of the input geometry if it is valid, or NA if it is invalid.

Examples

# requires GEOS >= 3.8, otherwise is only a validity test (see Note)
geos_version()

# valid
wkt <- "POINT (0 0)"
g_make_valid(wkt, as_wkb = FALSE)

# invalid to valid
wkt <- "POLYGON ((0 0,10 10,0 10,10 0,0 0))"
g_make_valid(wkt, as_wkb = FALSE)

# invalid - error
wkt <- "LINESTRING (0 0)"
g_make_valid(wkt)  # NA

Description

g_name() returns geometry type names in well known text format.

Usage

g_name(geom, quiet = FALSE)

Arguments

Value

character vector of the same length as the number of input geometries in geom, containing the WKT names for the corresponding geometies.

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_summary() returns text summaries of WKB/WKT geometries. Requires GDAL >= 3.7.

Usage

g_summary(geom, quiet = FALSE)

Arguments

Value

character vector of the same length as the number of input geometries in geom, containing summaries for the corresponding geometies.

Examples

# Requires GDAL >= 3.7
if (as.integer(gdal_version()[2]) >= 3070000) {
  f <- system.file("extdata/ynp_fires_1984_2022.gpkg", package = "gdalraster")
  lyr <- new(GDALVector, f, "mtbs_perims")

  feat <- lyr$getNextFeature()
  g_summary(feat$geom)

  feat_set <- lyr$fetch(5)
  g_summary(feat_set$geom)

  lyr$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

g_wk2wk() converts geometries between Well Known Binary (WKB) and Well Known Text (WKT) formats. A geometry given as a raw vector of WKB will be converted to a WKT string, while a geometry given as a WKT string will be converted to a WKB raw vector. Input may also be a list of WKB raw vectors or a character vector of WKT strings.

Usage

g_wk2wk(geom, as_iso = FALSE, byte_order = "LSB")

Arguments

Value

For input of a WKB raw vector or list of raw vectors, returns a character vector of WKT strings, with length of the returned vector equal to the number of input raw vectors. For input of a single WKT string, returns a raw vector of WKB. For input of a character vector containing more than one WKT string, returns a list of WKB raw vectors, with length of the returned list equal to the number of input strings.

Note

With as_iso = FALSE (the default), geometries are exported as extended dimension (Z) WKB/WKT for types Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon and GeometryCollection. For other geometry types, it is equivalent to ISO.

When the return value is a list of WKB raw vectors, an element in the returned list will contain NA if the corresponding input string was NA or empty ("").

When input is a list of WKB raw vectors, a corresponding element in the returned character vector will be an empty string ("") if the input was a raw vector of length 0 (raw(0)). If an input list element is not a raw vector, then the corresponding element in the returned character vector will be NA.

See Also

GEOS reference for geometry formats:
https://libgeos.org/specifications/

Examples

wkt <- "POINT (-114 47)"
wkb <- g_wk2wk(wkt)
print(wkb)
g_wk2wk(wkb)

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 (see Details)

## 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
ds$infoOptions
ds$quiet
ds$readByteAsRaw

## Methods
ds$getFilename()
ds$setFilename(filename)
ds$open(read_only)
ds$isOpen()
ds$getFileList()

ds$info()
ds$infoAsJSON()

ds$getDriverShortName()
ds$getDriverLongName()

ds$getRasterXSize()
ds$getRasterYSize()
ds$getRasterCount()

ds$addBand(dataType, options)

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$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$getMaskFlags(band)
ds$getMaskBand(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$setMetadata(band, metadata, 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(band, value, ivalue)

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

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

ds$flushCache()

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

ds$close()

Details

Constructors

new(GDALRaster, filename, read_only)
Returns an object of class GDALRaster. The read_only argument 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 access, 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. The shared argument defaults to TRUE but can be set to FALSE with this constructor (see Note). All arguments are required with this form of the constructor, but open_options can be NULL. Returns an object of class GDALRaster.

Read/write fields

$infoOptions
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
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
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.

Methods

$getFilename()
Returns a character string containing the filename associated with this GDALRaster object (filename originally used to open the dataset). May be a regular filename, database connection string, URL, etc.

$setFilename(filename)
Sets the filename if the underlying dataset does not already have an associated filename. Explicitly setting the filename is an advanced setting that should only be used when the user has determined that it is needed. Writing certain virtual datasets to file is one potential use case (e.g., a dataset returned by autoCreateWarpedVRT()).

$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.

$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).

$addBand(dataType, options)
Adds a band to a dataset if the underlying format supports this action. Most formats do not, but "MEM" and "VRT" are notable exceptions that support adding bands. The added band will always be the last band. dataType is a character string containing the data type name (e.g., "Byte", "Int16", "UInt16", "Int32", "Float32", etc). The options argument is a character vector of NAME=VALUE option strings. Supported options are format specific. Note that the options argument is required but may be given as NULL. Returns logical TRUE on success or FALSE if a band could not be added.

$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.

$getDescription(band)
Returns a string containing the description for band. An empty string is returned if no description is set for the band. Passing 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. (Passing band = 0 can be used to set the dataset-level description. Note that the dataset description is generally the filename that was used to open the dataset. It usually should not be changed by calling this method on an existing dataset.)

$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.

$getMaskFlags(band)
Returns the status flags of the mask band associated with band. Masks are represented as Byte bands with a value of zero indicating nodata and non-zero values indicating valid data. Normally the value 255 will be used for valid data pixels. GDAL supports external (.msk) mask bands, and normal Byte alpha (transparency) band as mask (any value other than 0 to be treated as valid data). But masks may not be regular raster bands on the datasource, such as an implied mask from a band nodata value or the ALL_VALID mask. See RFC 15: Band Masks for more details (https://gdal.org/development/rfc/rfc15_nodatabitmask.html.

Returns a named list of GDAL mask flags and their logical values, with the following definitions:

• ALL_VALID: There are no invalid pixels, all mask values will be 255. When used this will normally be the only flag set.

• PER_DATASET: The mask band is shared between all bands on the dataset.

• ALPHA: The mask band is actually an alpha band and may have values other than 0 and 255.

• NODATA: Indicates the mask is actually being generated from nodata values (mutually exclusive of ALPHA).

$getMaskBand(band)
Returns the mask filename and band number associated with band. The return value is a named list with two elements. The MaskFile element gives the filename where the mask band is located, e.g., a file with the same name as the main dataset but suffixed with .msk in the case of a GDAL external mask file. MaskFile will be an empty string for the derived ALL_VALID and NODATA masks, which internally are freestanding bands not considered to be a part of a dataset. The MaskBand element gives the band number for a mask that is a regular alpha band in the main dataset or external mask file. BandNumber will be 0 for the ALL_VALID and NODATA masks.

$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.

$setMetadata(band, metadata, domain)
Sets metadata in the specified domain. The metadata argument is given as a character vector of NAME=VALUE pairs. Pass band = 0 to set dataset-level metadata, or pass an integer band number to set band-level metadata. Use domain = "" (empty string) to set an item in the default domain. Returns logical TRUE on success or FALSE if metadata could not be set.

$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. Pass band = 0 to set dataset-level metadata, or pass an integer band number to set band-level metadata. Use domain = "" (empty string) to set an item in the default domain. Returns logical TRUE on success or FALSE if metadata could not be set.

$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).

$fillRaster(band, value, ivalue)
Fills band with a constant value. GDAL makes no guarantees about what values the pixels in newly created files are set to, so this method can be used to clear a band to a specified "default" value. The fill value is passed as numeric, but this will be converted to the underlying raster data type before writing to the file. The ivalue argument allows setting the imaginary component of a complex value. Note that ivalue is a required argument but can be set to 0 for real data types. No return value. An error is raised if the operation fails.

$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:

• Gray: column 2 contains grayscale values (columns 3:5 unused)

• RGB: columns 2:5 contain red, green, blue, alpha

• CMYK: columns 2:5 contain cyan, magenta, yellow, black

• HLS: columns 2:4 contain hue, lightness, saturation (column 5 unused)

$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.

$clearColorTable(band)
Clears the raster color table for band. Returns logical TRUE on success or FALSE if the color table could not be cleared, e.g., if this action is not supported by the driver.

$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() is a 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()


## 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

GDALVector provides an interface for accessing a vector layer in a GDAL dataset and calling methods on the underlying OGRLayer object. An object of class GDALVector persists an open connection to the dataset, and exposes methods for retrieving layer information, setting attribute and spatial filters, and reading/writing feature data. See https://gdal.org/api/index.html for details of the GDAL Vector API.

Class GDALVector is currently under development. An initial implementation supporting read access was added in gdalraster 1.11.1.9100. A working document with draft specifications is available at:
https://usdaforestservice.github.io/gdalraster/articles/gdalvector-draft.html
and discussion thread/status updates at:
https://github.com/USDAForestService/gdalraster/issues/241.

Arguments

Value

An object of class GDALVector which contains pointers to the opened layer and the dataset that owns it, and methods that operate on the layer as described in Details. GDALVector is a C++ class exposed directly to R (via RCPP_EXPOSED_CLASS). Fields and methods of the class are accessed using the $ operator. Note that all arguments to exposed class methods are required (but do not have to be named). The read/write fields are per-object settings which can be changed as needed during the lifetime of the object.

Usage (see Details)

## Constructors
# for single-layer file formats such as shapefile
lyr <- new(GDALVector, dsn)
# specifying the layer name, or SQL statement defining the layer
lyr <- new(GDALVector, dsn, layer)
# for update access
lyr <- new(GDALVector, dsn, layer, read_only = FALSE)
# using dataset open options
lyr <- new(GDALVector, dsn, layer, read_only, open_options)
# setting a spatial filter and/or specifying the SQL dialect
lyr <- new(GDALVector, dsn, layer, read_only, open_options, spatial_filter, dialect)

## Read/write fields (per-object settings)
lyr$defaultGeomFldName
lyr$promoteToMulti
lyr$quiet
lyr$returnGeomAs
lyr$wkbByteOrder

## Methods
lyr$open(read_only)
lyr$isOpen()
lyr$getDsn()
lyr$getFileList()
lyr$getDriverShortName()
lyr$getDriverLongName()

lyr$getName()
lyr$getFieldNames()
lyr$testCapability()
lyr$getFIDColumn()
lyr$getGeomType()
lyr$getGeometryColumn()
lyr$getSpatialRef()
lyr$bbox()
lyr$getLayerDefn()

lyr$setAttributeFilter(query)
lyr$getAttributeFilter()
lyr$setIgnoredFields(fields)
lyr$setSelectedFields(fields)

lyr$setSpatialFilter(wkt)
lyr$setSpatialFilterRect(bbox)
lyr$getSpatialFilter()
lyr$clearSpatialFilter()

lyr$getFeatureCount()
lyr$getNextFeature()
lyr$setNextByIndex(i)
lyr$getFeature(fid)
lyr$resetReading()
lyr$fetch(n)

lyr$setFeature(feature)
lyr$createFeature(feature)
lyr$upsertFeature(feature)
lyr$getLastWriteFID()
lyr$deleteFeature(fid)
lyr$syncToDisk()

lyr$startTransaction(force)
lyr$commitTransaction()
lyr$rollbackTransaction()

ds$getMetadata()
ds$setMetadata(metadata)
ds$getMetadataItem(mdi_name)

lyr$close()

Details

Constructors

new(GDALVector, dsn)
The first layer by index is assumed if the layer argument is omitted, so this form of the constructor might be used for single-layer formats like shapefile.

new(GDALVector, dsn, layer)
Constructor specifying the name of a layer to open. The layer argument may also be given as an SQL SELECT statement to define a layer as the result set.

new(GDALVector, dsn, layer, read_only)
Constructor specifying read/write access (⁠read_only = {TRUE|FALSE})⁠. The layer argument is required in this form of the constructor, but may be given as empty string (""), in which case the first layer by index will be assumed.

new(GDALVector, dsn, layer, read_only, open_options)
Constructor specifying dataset open options as a character vector of NAME=VALUE pairs.

new(GDALVector, dsn, layer, read_only, open_options, spatial_filter, dialect))
Constructor to specify a spatial filter and/or SQL dialect. All arguments are required in this form of the constructor, but open_options may be NULL, and spatial_filter or dialect may be an empty string ("").

Read/write fields

$defaultGeomFldName
Character string specifying a name to use for returned columns when the geometry column name in the source layer is empty, like with shapefiles etc. Defaults to "geometry".

$promoteToMulti
A logical value specifying whether to automatically promote geometries from Polygon to MultiPolygon, Point to MultiPoint, or LineString to MultiLineString during read operations (i.e., with methods ⁠$getFeature()⁠, ⁠$getNextFeature()⁠, ⁠$fetch()⁠). Defaults to FALSE. Setting to TRUE may be useful when reading from layers such as shapefiles that mix, e.g., Polygons and MultiPolygons.

$quiet
A logical value, FALSE by default. Set to TRUE to suppress various messages and warnings.

$returnGeomAs
Character string specifying the return format of feature geometries. Must be one of WKB (the default), WKB_ISO, WKT, WKT_ISO, BBOX, or NONE. Using WKB/WKT exports as 99-402 extended dimension (Z) types for Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon and GeometryCollection. For other geometry types, it is equivalent to using WKB_ISO/WKT_ISO (see https://libgeos.org/specifications/wkb/). Using BBOX exports as a list of numeric vectors, each of length 4 with values ⁠xmin, ymin, xmax, ymax⁠. If an empty geometry is encountered these values will be NA_real_ in the corresponding location. Using NONE will result in no geometry value being present in the feature returned.

$wkbByteOrder
Character string specifying the byte order for WKB geometries. Must be either LSB (Least Significant Byte first, the default) or MSB (Most Significant Byte first).

Methods

$open(read_only)
(Re-)opens the vector layer on the existing DSN. Use this method to open a layer that has been closed using $close(). May be used to re-open a layer 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 a logical scalar indicating whether the vector dataset is open.

$getDsn()
Returns a character string containing the dsn associated with this GDALVector object (dsn originally used to open the layer).

$getFileList()
Returns a character vector of files believed to be part of the data source. 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.

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

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

$getName()
Returns the layer name.

$getFieldNames()
Returns a character vector of the layer's field names.

$testCapability()
Tests whether the layer supports named capabilities based on the current read/write access. Returns a list of capabilities with values TRUE or FALSE. 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().)

$getFIDColumn()
Returns the name of the underlying database column being used as the FID column, or empty string ("") if not supported.

$getGeomType()
Returns the well known name of the layer geometry type as character string. For layers with multiple geometry fields, this method only returns the geometry type of the first geometry column. For other columns, use ⁠$getLayerDefn()⁠. For layers without any geometry field, this method returns "NONE".

$getGeometryColumn()
Returns he name of the underlying database column being used as the geometry column, or an empty string ("") if not supported. For layers with multiple geometry fields, this method only returns the name of the first geometry column. For other columns, use ⁠$getLayerDefn()⁠.

$getSpatialRef()
Returns a WKT string containing the spatial reference system for this layer, or empty string ("") if no spatial reference exists.

$bbox()
Returns a numeric vector of length four containing the bounding box for this layer (xmin, ymin, xmax, ymax). Note that bForce = true is set in the underlying API call to OGR_L_GetExtent(), so the entire layer may be scanned to compute a minimum bounding rectangle (see FastGetExtent in the list returned by ⁠$testCapability()⁠). Depending on the format driver, a spatial filter may or may not be taken into account, so it is safer to call ⁠$bbox()⁠ without setting a spatial filter.

$getLayerDefn()
Returns a list containing the OGR feature class definition for this layer (a.k.a. layer definition). The list contains zero or more attribute field definitions, along with one or more geometry field definitions. See ogr_define for details of the field and feature class definitions.

$setAttributeFilter(query)
Sets an attribute query string to be used when fetching features via the ⁠$getNextFeature()⁠ or ⁠$fetch()⁠ methods. Only features for which query evaluates as true will be returned. The query string should be in the format of an SQL WHERE clause, described in the "WHERE" section of the OGR SQL dialect documentation (e.g., "population > 1000000 and population < 5000000", where population is an attribute in the layer). In some cases (RDBMS backed drivers, SQLite, GeoPackage) the native capabilities of the database may be used to to interpret the WHERE clause, in which case the capabilities will be broader than those of OGR SQL. Note that installing a query string will generally result in resetting the current reading position (as with ⁠$resetReading()⁠ described below). The query parameter may be set to empty string ("") to clear the current attribute filter.

$getAttributeFilter()
Returns the attribute query string currently in use, or empty string ("") if an attribute filter is not set.

$setIgnoredFields(fields)
Set which fields can be omitted when retrieving features from the layer. The fields argument is a character vector of field names. Passing an empty string ("") for fields will reset to no ignored fields. If the format driver supports this functionality (testable using ⁠$testCapability()$IgnoreFields⁠), it will not fetch the specified fields in subsequent calls to ⁠$getFeature()⁠ / ⁠$getNextFeature()⁠ / ⁠$fetch()⁠, and thus save some processing time and/or bandwidth. Besides field names of the layer, the following special fields can be passed: "OGR_GEOMETRY" to ignore geometry and "OGR_STYLE" to ignore layer style. By default, no fields are ignored. Note that fields that are used in an attribute filter should generally not be set as ignored fields, as most drivers (such as those relying on the OGR SQL engine) will be unable to correctly evaluate the attribute filter. No return value, called for side effects.

$setSelectedFields(fields)
Set which fields will be included when retrieving features from the layer. The fields argument is a character vector of field names. Passing an empty string ("") for fields will reset to no ignored fields. See the ⁠$setIgnoredFields()⁠ method above for more information. The data source must provide IgnoreFields capability in order to set selected fields. Note that geometry fields, if desired, must be specified when setting selected fields, either by including named geometry field(s) or the special field "OGR_GEOMETRY" in the fields argument. No return value, called for side effects.

$setSpatialFilter(wkt)
Sets a new spatial filter from a geometry in WKT format. This method sets the geometry to be used as a spatial filter when fetching features via the ⁠$getNextFeature()⁠ or ⁠$fetch()⁠ methods. Only features that geometrically intersect the filter geometry will be returned. Currently this test may be inaccurately implemented (depending on the vector format driver), but it is guaranteed that all features whose envelope overlaps the envelope of the spatial filter will be returned. This can result in more shapes being returned that should strictly be the case. wkt is a character string containing a WKT geometry in the same coordinate system as the layer. An empty string ("") may be passed indicating that the current spatial filter should be cleared, but no new one instituted.

$setSpatialFilterRect(bbox)
Sets a new rectangular spatial filter. This method sets a rectangle to be used as a spatial filter when fetching features via the ⁠$getNextFeature()⁠ or ⁠$fetch()⁠ methods. Only features that geometrically intersect the given rectangle will be returned. bbox is a numeric vector of length four containing xmin, ymin, xmax, ymax in the same coordinate system as the layer as a whole (as returned by ⁠$getSpatialRef()⁠).

$getSpatialFilter()
Returns the current spatial filter geometry as a WKT string, or empty string ("") if a spatial filter is not set.

$clearSpatialFilter()
Clears a spatial filter that was set with ⁠$setSpatialFilterRect()⁠. No return value, called for that side effect.

$getFeatureCount()
Returns the number of features in the layer. For dynamic databases the count may not be exact. This method forces a count in the underlying API call (i.e., bForce = TRUE in the call to OGR_L_GetFeatureCount()). Note that some vector drivers will actually scan the entire layer once to count features. The FastFeatureCount element in the list returned by the ⁠$testCapability()⁠ method can be checked if this might be a concern. The number of features returned takes into account the spatial and/or attribute filters. Some driver implementations of this method may alter the read cursor of the layer.

$getNextFeature()
Fetch the next available feature from this layer. Only features matching the current spatial and/or attribute filter (if defined) will be returned. This method implements sequential access to the features of a layer. The ⁠$resetReading()⁠ method can be used to start at the beginning again. Returns a list with the unique feature identifier (FID), the attribute and geometry field names, and their values. The returned list carries the OGRFeature class attribute with S3 methods for for print() and plot(). NULL is returned if no more features are available.

$setNextByIndex(i)
Moves the read cursor to the ith feature in the current result set (with 0-based indexing). This method allows positioning of a layer such that a call to ⁠$getNextFeature()⁠ or fetch() will read the requested feature(s), where i is an absolute index into the current result set. So, setting i = 3 would mean the next feature read with ⁠$getNextFeature()⁠ would have been the 4th feature read if sequential reading took place from the beginning of the layer, including accounting for spatial and attribute filters. This method is not implemented efficiently by all vector format drivers. The default implementation simply resets reading to the beginning and then calls GetNextFeature() i times. To determine if fast seeking is available on the current layer, check the FastSetNextByIndex element in the list returned by the ⁠$testCapability()⁠ method. No return value, called for side effect.

$getFeature(fid)
Returns a feature by its identifier. The value of fid must be a numeric scalar, optionally carrying the bit64::integer64 class attribute. Success or failure of this operation is unaffected by any spatial or attribute filters that may be in effect. The RandomRead element in the list returned by ⁠$testCapability()⁠ can be checked to establish if this layer supports efficient random access reading; however, the call should always work if the feature exists since a fallback implementation just scans all the features in the layer looking for the desired feature. Returns a list with the unique feature identifier (FID), the attribute and geometry field names, and their values, or NULL on failure. Note that sequential reads (with ⁠$getNextFeature()⁠) are generally considered interrupted by a call to ⁠$getFeature()⁠.

$resetReading()
Reset feature reading to start on the first feature. No return value, called for that side effect.

$fetch(n)
Fetches the next n features from the layer and returns them as a data frame. This allows retrieving the entire set of features, one page of features at a time, or the remaining features (from the current cursor position). Returns a data frame with as many rows as features were fetched, and as many columns as attribute plus geometry fields in the result set, even if the result is a single value or has one or zero rows. The returned data frame carries the OGRFeatureSet class attribute with S3 methods for for print() and plot().

This method is an analog of DBI::dbFetch().

The n argument is the maximum number of features to retrieve per fetch given as integer or numeric but assumed to be a whole number (will be truncated). Use n = -1 or n = Inf to retrieve all pending features (resets reading to the first feature). Otherwise, ⁠$fetch()⁠ can be called multiple times to perform forward paging from the current cursor position. Passing n = NA is also supported and returns the remaining features. Fetching zero features is possible to retrieve the structure of the feature set as a data frame (columns fully typed).

OGR field types are returned as the following R types (NA for OGR NULL values):

• OFTInteger: integer

• OFTInteger subtype OFSTBoolean: logical

• OFTIntegerList: vector of integer (list column)

• OFTInteger64: numeric carrying "integer64" class attribute {bit64}

• OFTInteger64 subtype OFSTBoolean: logical

• OFTInteger64List: vector of bit64::integer64 (list column)

• OFTReal: numeric

• OFTRealList: vector of numeric (list column)

• OFTString: character string

• OFTStringList: vector of character strings (list column)

• OFTDate: class "Date" (numeric)

• OFTDateTime: class "POSIXct" (numeric, millisecond accuracy)

• OFTTime: character string ("HH:MM:SS")

• OFTBinary: raw vector (list column, NULL entries for OGR NULL values)

Geometries are not returned if the field returnGeomAs is set to NONE. Omitting the geometries may be beneficial for performance and memory usage when access only to feature attributes is needed. Geometries are returned as raw vectors in a data frame list column when returnGeomAs is set to WKB (the default) or WKB_ISO, or as character strings when returnGeomAs is set to one of WKT or WKT_ISO.

Note that ⁠$getFeatureCount()⁠ is called internally when fetching the full feature set or all remaining features (but not for a page of features).

$setFeature(feature)
Rewrites/replaces an existing feature. This method writes a feature based on the feature id within the input feature. The feature argument is a named list of fields and their values, and must include a ⁠$FID⁠ element referencing the existing feature to rewrite. The RandomWrite element in the list returned by ⁠$testCapability()⁠ can be checked to establish if this layer supports random access writing via ⁠$setFeature()⁠. The way omitted fields in the passed feature are processed is driver dependent:

• SQL-based drivers which implement set feature through SQL UPDATE will skip unset fields, and thus the content of the existing feature will be preserved.

• The shapefile driver will write a NULL value in the DBF file.

• The GeoJSON driver will take into account unset fields to remove the corresponding JSON member.

Returns logical TRUE upon successful completion, or FALSE if setting the feature did not succeed. The FID of the last feature written to the layer may be obtained with the method ⁠$getLastWriteFID()⁠ (see below). To set a feature, but create it if it doesn't exist see the ⁠$upsertFeature()⁠ method.

$createFeature(feature)
Creates and writes a new feature within the layer. The feature argument is a named list of fields and their values. The passed feature is written to the layer as a new feature, rather than overwriting an existing one. If the feature has a ⁠$FID⁠ element other than NA, then the vector format driver may use that as the feature id of the new feature, but not necessarily. The FID of the last feature written to the layer may be obtained with the method ⁠$getLastWriteFID()⁠ (see below). Returns logical TRUE upon successful completion, or FALSE if creating the feature did not succeed. To create a feature, but set it if it already exists see the ⁠$upsertFeature()⁠ method.

$upsertFeature(feature)
Rewrites/replaces an existing feature or creates a new feature within the layer. This method will write a feature to the layer, based on the feature id within the input feature. The feature argument is a named list of fields and their values, potentially including a ⁠$FID⁠ element referencing an existing feature to rewrite. If the feature id doesn't exist a new feature will be written. Otherwise, the existing feature will be rewritten. The UpsertFeature element in the list returned by ⁠$testCapability()⁠ can be checked to determine if this layer supports upsert writing. See ⁠$setFeature()⁠ above for a description of how omitted fields in the passed feature are processed. Returns logical TRUE upon successful completion, or FALSE if upserting the feature did not succeed. Requires GDAL >= 3.6.

$getLastWriteFID()
Returns the FID of the last feature written (either newly created or updated existing). NULL is returned if no features have been written in the layer. Note that OGRNullFID (-1) may be returned after writing a feature in some formats. This is the case if a FID has not been assigned yet, and generally does not indicate an error (e.g., formats that do not store a persistent FID and assign FIDs upon a sequential read operation). The returned FID is a numeric scalar carrying the bit64::integer64 class attribute.

$deleteFeature(fid)
Deletes a feature from the layer. The feature with the indicated feature ID is deleted from the layer if supported by the format driver. The value of fid must be a numeric scalar, optionally carrying the bit64::integer64 class attribute (should be a whole number, will be truncated). The DeleteFeature element in the list returned by ⁠$testCapability()⁠ can be checked to establish if this layer has delete feature capability. Returns logical TRUE if the operation succeeds, or FALSE on failure.

$syncToDisk()
Flushes pending changes to disk. This call is intended to force the layer to flush any pending writes to disk, and leave the disk file in a consistent state. It would not normally have any effect on read-only datasources. Some formats do not implement this method, and will still return no error. An error is only returned if an error occurs while attempting to flush to disk. In any event, you should always close any opened datasource with ⁠$close()⁠ which will ensure all data is correctly flushed. Returns logical TRUE if no error occurs (even if nothing is done) or FALSE on error.

$startTransaction(force)
Creates a transaction if supported by the vector data source. The force argument is a logical value. If force = FALSE, only "efficient" transactions will be attempted. Some drivers may offer an emulation of transactions, but sometimes with significant overhead, in which case the user must explicitly allow for such an emulation by setting force =TRUE. The function ogr_ds_test_cap() can be used to determine whether a vector data source supports efficient or emulated transactions.

All changes done after the start of the transaction are definitely applied in the data source if ⁠$commitTransaction()⁠ is called. They can be canceled by calling rollbackTransaction() instead. Nested transactions are not supported. Transactions are implemented at the dataset level, so multiple GDALVector objects using the same data source should not have transactions active at the same time.

In case ⁠$startTransaction()⁠ fails, neither ⁠$commitTransaction()⁠ nor ⁠$rollbackTransaction()⁠ should be called. If an error occurs after a successful ⁠$startTransaction()⁠, the whole transaction may or may not be implicitly canceled, depending on the format driver (e.g., the PostGIS driver will cancel it, SQLite/GPKG will not). In any case, in the event of an error, an explicit call to rollbackTransaction() should be done to keep things balanced.

Returns logical TRUE if the transaction is created, or FALSE on failure.

$commitTransaction()
Commits a transaction if supported by the vector data source. Returns a logical value, TRUE if the transaction is successfully committed. Returns FALSE if no transaction is active, or the rollback fails, or if the data source does not support transactions. Depending on the format driver, this may or may not abort layer sequential reading that may be active.

$rollbackTransaction()
Rolls back a data source to its state before the start of the current transaction, if transactions are supported by the data source. Returns a logical value, TRUE if the transaction is successfully rolled back. Returns FALSE if no transaction is active, or the rollback fails, or if the data source does not support transactions.

$getMetadata()
Returns a character vector of all metadata NAME=VALUE pairs for the layer or empty string ("") if there are no metadata items.

$setMetadata(metadata)
Sets metadata on the layer if the format supports it. The metadata argument is given as a character vector of NAME=VALUE pairs. Returns logical TRUE on success or FALSE if metadata could not be set.

$getMetadataItem(mdi_name)
Returns the value of a specific metadata item named mdi_name, or empty string ("") if no matching item is found.

$close()
Closes the vector dataset (no return value, called for side effects). Calling $close() results in proper cleanup, and flushing of any pending writes. The GDALVector object is still available after calling $close(). The layer can be re-opened on the existing dsn with $open(read_only = {TRUE|FALSE}).

See Also

ogr_define, ogr_manage, ogr2ogr(), ogrinfo()

GDAL vector format descriptions:
https://gdal.org/drivers/vector/index.html

GDAL-supported SQL dialects:
https://gdal.org/user/ogr_sql_sqlite_dialect.html)

Examples

## MTBS fire perimeters in Yellowstone National Park 1984-2022
f <- system.file("extdata/ynp_fires_1984_2022.gpkg", package = "gdalraster")

## copy to a temporary file that is writeable
dsn <- file.path(tempdir(), basename(f))
file.copy(f, dsn)

lyr <- new(GDALVector, dsn, "mtbs_perims")

## object of class GDALVector
lyr
str(lyr)

## dataset info
lyr$getDriverShortName()
lyr$getDriverLongName()
lyr$getFileList()

## layer info
lyr$getName()
lyr$getGeomType()
lyr$getGeometryColumn()
lyr$getFIDColumn()
lyr$getSpatialRef()
lyr$bbox()

## layer capabilities
lyr$testCapability()

## re-open with write access
lyr$open(read_only = FALSE)
lyr$testCapability()$SequentialWrite
lyr$testCapability()$RandomWrite

## feature class definition - a list of field names and their definitions
defn <- lyr$getLayerDefn()
names(defn)
str(defn)

## default value of the read/write field 'returnGeomAs'
lyr$returnGeomAs

lyr$getFeatureCount()

## sequential read cursor
feat <- lyr$getNextFeature()
# a list of field names and their values, with class attribute `OGRFeature`
feat

## set an attribute filter
lyr$setAttributeFilter("ig_year = 2020")
lyr$getFeatureCount()

feat <- lyr$getNextFeature()
plot(feat)

## NULL when no more features are available
lyr$getNextFeature()

## reset reading to the start
lyr$resetReading()
lyr$getNextFeature()

## clear the attribute filter
lyr$setAttributeFilter("")
lyr$getFeatureCount()

## set a spatial filter
## get the bounding box of the largest 1988 fire and use as spatial filter
## first set a temporary attribute filter to do the lookup
lyr$setAttributeFilter("ig_year = 1988 ORDER BY burn_bnd_ac DESC")
feat <- lyr$getNextFeature()
feat

bbox <- g_wk2wk(feat$geom) |> bbox_from_wkt()

## set spatial filter on the full layer
lyr$setAttributeFilter("")  # clears
lyr$setSpatialFilterRect(bbox)
lyr$getFeatureCount()

## fetch in chunks and return as data frame (class `OGRFeatureSet`)
feat_set <- lyr$fetch(20)
head(feat_set)
plot(feat_set)

## the next chunk
feat_set <- lyr$fetch(20)
nrow(feat_set)

## no features remaining
feat_set <- lyr$fetch(20)
nrow(feat_set)
str(feat_set)  # 0-row data frame with columns typed

## fetch all pending features
feat_set <- lyr$fetch(-1)  # resets reading to the first feature
nrow(feat_set)
plot(feat_set)

lyr$clearSpatialFilter()
lyr$getFeatureCount()

lyr$close()


## simple example for feature write methods showing use of various data types
## create and write to a new layer in a GeoPackage data source
dsn2 <- tempfile(fileext = ".gpkg")

## define a feature class
defn <- ogr_def_layer("POINT", srs = epsg_to_wkt(4326))

## add field definitions
defn$unique_int <- ogr_def_field("OFTInteger", is_nullable = FALSE,
                                 is_unique = TRUE)
defn$bool_data <- ogr_def_field("OFTInteger", fld_subtype = "OFSTBoolean")
defn$large_ints <- ogr_def_field("OFTInteger64")
defn$doubles <- ogr_def_field("OFTReal")
defn$strings <- ogr_def_field("OFTString", fld_width = 50)
defn$dates <- ogr_def_field("OFTDate")
defn$dt_modified <- ogr_def_field("OFTDateTime",
                                  default_value = "CURRENT_TIMESTAMP")
defn$blobs <- ogr_def_field("OFTBinary")

ogr_ds_create("GPKG", dsn2, "test_layer", layer_defn = defn)

lyr <- new(GDALVector, dsn2, "test_layer", read_only = FALSE)
# lyr$getLayerDefn() |> str()

## define a feature to write
feat1 <- list()
## $FID is omitted since it is assigned when written (could also be NA)
## $dt_modified is omitted since the datasource sets a default timestamp
feat1$unique_int <- 1001
feat1$bool_data <- TRUE
## passing a string to as.integer64()
## this value is too large to be represented exactly as R numeric (double)
feat1$large_ints <- bit64::as.integer64("90071992547409910")
feat1$doubles <- 1.234
feat1$strings <- "A test string"
feat1$dates <- as.Date("2024-01-01")
feat1$blobs <- charToRaw("A binary object")
feat1$geom <- "POINT (1 1)"  # can be a WKT string or raw vector of WKB

## create as a new feature in the layer
lyr$createFeature(feat1)

## the assigned FID
lyr$getLastWriteFID()

## this fails due to the unique constraint
lyr$createFeature(feat1)

feat2 <- list()
feat2$unique_int <- 1002
feat2$bool_data <- FALSE
feat2$large_ints <- bit64::as.integer64("90071992547409920")
feat2$doubles <- 2.345
feat2$strings <- "A test string 2"
feat2$dates <- as.Date("2024-01-02")
feat2$blobs <- charToRaw("A binary object 2")
feat2$geom <- "POINT (2 2)"

lyr$createFeature(feat2)
lyr$getLastWriteFID()

## close and re-open as a read-only layer
lyr$open(read_only = TRUE)

lyr$getFeatureCount()
feat_set <- lyr$fetch(-1)  # -1 for all features reading from start
str(feat_set)

## edit an existing feature, e.g., feat <- lyr$getFeature(2)
## here we copy a row of the data frame returned by lyr$fetch() above
feat <- feat_set[2,]
str(feat)

Sys.sleep(1)  # only to ensure a timestamp difference

feat$bool_data <- TRUE
feat$strings <- paste(feat$strings, "- edited")
feat$dt_modified <- Sys.time()
feat$geom <- "POINT (2.001 2.001)"

lyr$open(read_only = FALSE)

## lyr$setFeature() re-writes the feature identified by the $FID element
## N.B., all fields are re-written:
##   any fields omitted from the input feature, or set to NA, will be
##   re-written as OGR NULL
lyr$setFeature(feat)

lyr$open(read_only = TRUE)
lyr$getFeatureCount()

lyr$returnGeomAs <- "WKT"
feat_set <- lyr$fetch(-1)
str(feat_set)

lyr$close()

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()