# The examples make use of the FARSITE LCP format specification at:
# https://gdal.org/drivers/raster/lcp.html
# An LCP file is a raw format with a 7,316-byte header. The format
# specification gives byte offets and data types for fields in the header.
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
# identify a FARSITE v.4 LCP file
# function to check if the first three fields have valid data
# input is the first twelve raw bytes in the file
is_lcp <- function(bytes) {
values <- readBin(bytes, "integer", n = 3)
if ((values[1] == 20 || values[1] == 21) &&
(values[2] == 20 || values[2] == 21) &&
(values[3] >= -90 && values[3] <= 90)) {
return(TRUE)
} else {
return(FALSE)
}
}
vf <- new(VSIFile, lcp_file)
vf$read(12) |> is_lcp()
vf$tell()
# read the whole file into memory
bytes <- vf$ingest(-1)
vf$close()
# write to a VSI in-memory file
mem_file <- "/vsimem/storml_copy.lcp"
vf <- new(VSIFile, mem_file, "w")
vf$write(bytes)
vf$tell()
vf$rewind()
vf$tell()
vf$seek(0, SEEK_END)
(vf$tell() == vsi_stat(lcp_file, "size")) # TRUE
vf$rewind()
vf$read(12) |> is_lcp()
# read/write an integer field
# write invalid data for the Latitude field and then set back
# save the original first
vf$seek(8, SEEK_SET)
lat_orig <- vf$read(4)
readBin(lat_orig, "integer") # 46
# latitude -99 out of range
vf$seek(8, SEEK_SET)
writeBin(-99L, raw()) |> vf$write()
vf$rewind()
vf$read(12) |> is_lcp() # FALSE
vf$seek(8, SEEK_SET)
vf$read(4) |> readBin("integer") # -99
# set back to original
vf$seek(8, SEEK_SET)
vf$write(lat_orig)
vf$rewind()
vf$read(12) |> is_lcp() # TRUE
# read a vector of doubles - xmax, xmin, ymax, ymin
# 327766.1, 323476.1, 5105082.0, 5101872.0
vf$seek(4172, SEEK_SET)
vf$read(32) |> readBin("double", n = 4)
# read a short int, the canopy cover units
vf$seek(4232, SEEK_SET)
vf$read(2) |> readBin("integer", size = 2) # 1 = "percent"
# read the Description field
vf$seek(6804, SEEK_SET)
bytes <- vf$read(512)
rawToChar(bytes)
# edit the Description
desc <- paste(rawToChar(bytes),
"Storm Lake AOI,",
"Beaverhead-Deerlodge National Forest, Montana.")
vf$seek(6804, SEEK_SET)
charToRaw(desc) |> vf$write()
vf$close()
# verify the file as a raster dataset
ds <- new(GDALRaster, mem_file)
ds$info()
# retrieve Description from the metadata
# band = 0 for dataset-level metadata, domain = "" for default domain
ds$getMetadata(band = 0, domain = "")
ds$getMetadataItem(band = 0, mdi_name = "DESCRIPTION", domain = "")
ds$close()
vsi_unlink(mem_file)
Run the code above in your browser using DataLab