Get DPIRD CVT Zones as an {sf} Object in Your R Session

dpird_cvt_sf.R

# get DPIRD CVT zones from the Public Services Slip WA as an {sf} object in R
# data are CC 4.0 and available from: <https://data.wa.gov.au/slip>
wfs_wa <-
  "https://public-services.slip.wa.gov.au/public/services/SLIP_Public_Services/Boundaries_WFS/MapServer/WFSServer"

url <- httr2::url_parse(wfs_wa)
url$query <- list(service = "wfs",
                  #version = "2.0.0", # facultative
                  request = "GetCapabilities")
request <- httr2::url_build(url)
wa_client <- ows4R::WFSClient$new(wfs_wa,
                                  serviceVersion = "2.0.0")

wa_features <- wa_client$getFeatureTypes(pretty = TRUE)
wa_features <-
  gsub("SLIP_Public_Services_Boundaries_WFS:", "", wa_features$name)
dpird_features <- grep("DPIRD", wa_features, value = TRUE)

#' Get a named data set from WA's Web Feature Service
#' @param x The name of the desired data set to fetch.
#'
get_wfs_data <- function(x) {
  url$query <- list(
    service = "wfs",
    version = "2.0.0",
    request = "GetFeature",
    typename = x,
    srsName = "EPSG:4326"
  )
  request <- httr2::url_build(url)
  sf::read_sf(request)
}

cvt <- get_wfs_data("esri:Crop_Variety_Testing_Zones_of_WA__DPIRD-028_")

library("ggplot2")
ggplot(cvt) + 
  geom_sf()

If you want to get all DPIRD boundary data sets, use this on line 35.

# get all the data sets for DPIRD
dpird_data <- lapply(X = dpird_features, FUN = get_wfs_data)

# tidy up the names
names(dpird_data) <- tolower(sub("\\__[^.]*$", "", dpird_features))