library(soilDB)
library(terra)
#> terra 1.6.40
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
data(us_states, package = "spData")
# query soil data access for dominant condition (by mapunit) parent material information for NY
x <- get_SDA_pmgroupname(WHERE = "areasymbol LIKE 'NY%'",
method = "Dominant Condition",
simplify = FALSE)
# regular expression to find carbonate-containing parent materials
x$is_carbonate <- grepl("limestone|marble|carbonate", x$pmgroupname)
r <- mukey.wcs(sf::st_as_sf(subset(us_states, NAME == "New York")), res = 800)
levels(r) <- NULL
# replace MUKEY values with is_carbonate
r2 <- subst(r, x$mukey, x$is_carbonate, others = NA)
plot(r2)
# inspect cell coordinates
head(crds(r2))
#> x y
#> [1,] 1766123 2658356
#> [2,] 1762923 2657556
#> [3,] 1763723 2657556
#> [4,] 1764523 2657556
#> [5,] 1765323 2657556
#> [6,] 1766123 2657556
# write result to file
writeRaster(r2, filename = "ny_carbonate_parent_material.tif", overwrite = TRUE)
# you can read geotiff from file as follows
r3 <- rast("ny_carbonate_parent_material.tif")
# you can also pass the filename argument to many terra functions
r4 <- project(r3, "OGC:CRS84", filename="ny_carbonate_parent_material_wgs84.tif")
r4
#> class : SpatRaster
#> dimensions : 835, 1267, 1 (nrow, ncol, nlyr)
#> resolution : 0.008246337, 0.008246337 (x, y)
#> extent : -80.76919, -70.32108, 39.29279, 46.17848 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source : ny_carbonate_parent_material_wgs84.tif
#> name : mukey
#> min value : 0
#> max value : 1
# example extracting data at points
my_points <- data.frame(x = -75, y = 42)
my_point_vector <- vect(my_points, geom = c("x", "y"), crs = "OGC:CRS84")
my_point_vector
#> class : SpatVector
#> geometry : points
#> dimensions : 1, 0 (geometries, attributes)
#> extent : -75, -75, 42, 42 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
extract(r4, my_point_vector)
#> ID mukey
#> 1 1 0
