Getting started with readwfs

Source: vignettes/readwfs.Rmd

readwfs reads vector features from OGC Web Feature Services (WFS), OGC API Features (OAPIF), and ArcGIS REST endpoints. It returns plain tibbles with wk geometry columns and uses gdalraster for all I/O.

The package is designed to make navigating web feature services less painful – discovering what layers exist, what fields they have, and pulling just the features you need into R.

Finding services

readwfs ships with a small catalogue of known public services:

library(readwfs)
wfs_services()
#> # A tibble: 2 × 7
#>   name                   url               driver region description srs   notes
#>   <chr>                  <chr>             <chr>  <chr>  <chr>       <chr> <chr>
#> 1 LIST Tasmania          https://services… WFS    Tasma… Cadastral … EPSG… Rich…
#> 2 Esri SampleWorldCities https://samplese… WFS    Global Continents… EPSG… Alwa…

This is deliberately minimal. The real world has thousands of WFS endpoints run by government agencies, utilities, and research organisations. Most are poorly documented. readwfs helps you explore them.

Discovering layers

Point wfs_layers() at any WFS endpoint to see what’s available:

url <- wfs_example_url("list_tasmania")
layers <- wfs_layers(url, version = "2.0.0", srs = "EPSG:28355")
length(layers)
#> [1] 56
head(layers)
#> [1] "Public_OpenDataWFS:Climate_Futures_Tasmania_-_Scenario_A2"     
#> [2] "Public_OpenDataWFS:Climate_Futures_Tasmania_-_Scenario_B1"     
#> [3] "Public_OpenDataWFS:Interim_Planning_Scheme_-_Overlays"         
#> [4] "Public_OpenDataWFS:Interim_Planning_Scheme_-_Zoning_Boundaries"
#> [5] "Public_OpenDataWFS:Interim_Planning_Scheme_-__Zoning"          
#> [6] "Public_OpenDataWFS:LIDAR_Climate_Futures_Index"

Tasmania’s LIST service has over 100 layers. wfs_find_layers() narrows things down:

wfs_find_layers(url, "CADASTRAL|TASVEG|LOCAL_GOV",
                version = "2.0.0", srs = "EPSG:28355")
#> [1] "Public_OpenDataWFS:LIST_Cadastral_Parcels"        
#> [2] "Public_OpenDataWFS:LIST_Local_Government_Areas"   
#> [3] "Public_OpenDataWFS:LIST_Local_Government_Reserves"
#> [4] "Public_OpenDataWFS:TASVEG_3.0"                    
#> [5] "Public_OpenDataWFS:TASVEG_4.0"

Inspecting layers before downloading

wfs_layer_info() tells you the geometry type, feature count, and extent without downloading any features:

wfs_layer_info(
  url,
  layers = c(
    "Public_OpenDataWFS:LIST_CADASTRAL_PARCELS",
    "Public_OpenDataWFS:TASVEG_4.0"
  ),
  version = "2.0.0", srs = "EPSG:28355"
)
#> # A tibble: 2 × 9
#>   name   geom_column geom_type feature_count   xmin   ymin   xmax   ymax srs_wkt
#>   <chr>  <chr>       <chr>             <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  
#> 1 Publi… SHAPE       MULTISUR…        437755 2.25e5 5.14e6 6.30e5 5.66e6 "PROJC…
#> 2 Publi… SHAPE       MULTISUR…        463915 2.25e5 5.16e6 6.30e5 5.66e6 "PROJC…

wfs_fields() shows the attribute schema:

wfs_fields(url, "Public_OpenDataWFS:LIST_CADASTRAL_PARCELS",
           version = "2.0.0", srs = "EPSG:28355")
#> # A tibble: 16 × 7
#>    name       type       subtype   width precision is_nullable is_unique
#>    <chr>      <chr>      <chr>     <int>     <int> <lgl>       <lgl>    
#>  1 gml_id     OFTString  OFSTNone      0         0 FALSE       FALSE    
#>  2 OBJECTID   OFTInteger OFSTNone      0         0 TRUE        FALSE    
#>  3 CID        OFTInteger OFSTNone      0         0 TRUE        FALSE    
#>  4 VOLUME     OFTString  OFSTNone      8         0 TRUE        FALSE    
#>  5 FOLIO      OFTInteger OFSTInt16     0         0 TRUE        FALSE    
#>  6 PID        OFTInteger OFSTNone      0         0 TRUE        FALSE    
#>  7 POT_PID    OFTInteger OFSTNone      0         0 TRUE        FALSE    
#>  8 LPI        OFTString  OFSTNone      7         0 TRUE        FALSE    
#>  9 CAD_TYPE1  OFTString  OFSTNone     60         0 TRUE        FALSE    
#> 10 CAD_TYPE2  OFTString  OFSTNone     60         0 TRUE        FALSE    
#> 11 TENURE_TY  OFTString  OFSTNone     60         0 TRUE        FALSE    
#> 12 FEAT_NAME  OFTString  OFSTNone     60         0 TRUE        FALSE    
#> 13 STRATA_LEV OFTString  OFSTNone     60         0 TRUE        FALSE    
#> 14 COMP_AREA  OFTReal    OFSTNone      0         0 TRUE        FALSE    
#> 15 MEAS_AREA  OFTReal    OFSTNone      0         0 TRUE        FALSE    
#> 16 UFI        OFTString  OFSTNone     12         0 TRUE        FALSE

Reading features

wfs_read() fetches features into a tibble. Use bbox to spatially filter and max_features to cap the request:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

parcels <- wfs_read(
  url,
  layer = "Public_OpenDataWFS:LIST_CADASTRAL_PARCELS",
  bbox = wfs_example_bbox("sandy_bay"),
  srs = "EPSG:28355",
  max_features = 20
)
#> Reading 'Public_OpenDataWFS:LIST_CADASTRAL_PARCELS': 464 features available, geometry column 'SHAPE'
#>   20 features returned (capped at max_features = 20, 444 more available; use max_features = Inf to get all)
parcels
#> # A tibble: 20 × 18
#>        FID gml_id    OBJECTID    CID VOLUME FOLIO    PID POT_PID LPI   CAD_TYPE1
#>    <int64> <chr>        <int>  <int> <chr>  <int>  <int>   <int> <chr> <chr>    
#>  1    2688 LIST_Cad…     2688 1.10e6 59944      4 5.68e6       0 GRF64 Private …
#>  2    2877 LIST_Cad…     2877 1.02e6 128641     1 1.82e6       0 FMR21 Private …
#>  3    3425 LIST_Cad…     3425 1.15e6 126957     1 3.27e6       0 FJA84 Authorit…
#>  4    4682 LIST_Cad…     4682 1.02e6 128641     2 1.82e6       0 FMR21 Private …
#>  5    6102 LIST_Cad…     6102 1.25e6 252507     1 7.57e6       0 FAC83 Private …
#>  6    6287 LIST_Cad…     6287 1.02e6 127656     3 1.79e6       0 FBF18 Private …
#>  7    6395 LIST_Cad…     6395 1.50e6 173638     1 9.17e6       0 GAV70 Private …
#>  8    8685 LIST_Cad…     8685 1.41e6 231575     1 0            0 NA    Casement 
#>  9   10819 LIST_Cad…    10819 1.33e6 152956     1 2.82e6       0 GDH05 Private …
#> 10   11600 LIST_Cad…    11600 1.09e6 58505      2 5.68e6       0 HSF79 Private …
#> 11   13561 LIST_Cad…    13561 1.02e6 127656     2 1.79e6       0 FBF18 Private …
#> 12   14096 LIST_Cad…    14096 1.02e6 127656     1 1.79e6       0 FBF18 Private …
#> 13   16589 LIST_Cad…    16589 1.15e6 NA        NA 0            0 NA    Casement 
#> 14   17354 LIST_Cad…    17354 1.10e6 59944      2 5.68e6       0 GRF64 Private …
#> 15   18846 LIST_Cad…    18846 1.01e6 60155     27 5.57e6       0 GAV68 Private …
#> 16   19312 LIST_Cad…    19312 1.01e6 226698     6 5.56e6       0 GAV57 Private …
#> 17   19408 LIST_Cad…    19408 1.01e6 7826       1 1.48e6       0 GAW20 Private …
#> 18   20464 LIST_Cad…    20464 1.01e6 60754      5 5.60e6       0 GAV95 Private …
#> 19   20516 LIST_Cad…    20516 1.01e6 59722     72 5.58e6       0 GAS73 Private …
#> 20   20534 LIST_Cad…    20534 1.01e6 218685     1 5.57e6       0 GAM93 Private …
#> # ℹ 8 more variables: CAD_TYPE2 <chr>, TENURE_TY <chr>, FEAT_NAME <chr>,
#> #   STRATA_LEV <chr>, COMP_AREA <dbl>, MEAS_AREA <dbl>, UFI <chr>,
#> #   geometry <wk_wkb>

The geometry column is a wk::wkb vector – wk’s interchange format for spatial data in R:

class(parcels$geometry)
#> [1] "wk_wkb"  "wk_vctr"
wk::wk_plot(parcels$geometry)

Multiple layers, same area

Fetch several layers from the same service for the same bbox:

bbox <- wfs_example_bbox("sandy_bay")

tasveg <- wfs_read(
  url,
  layer = "Public_OpenDataWFS:TASVEG_4.0",
  bbox = bbox, srs = "EPSG:28355"
)
#> Reading 'Public_OpenDataWFS:TASVEG_4.0': 10 features available, geometry column 'SHAPE'
#>   10 features returned

lga <- wfs_read(
  url,
  layer = "Public_OpenDataWFS:LIST_LOCAL_GOVERNMENT_AREAS",
  bbox = bbox, srs = "EPSG:28355"
)
#> Reading 'Public_OpenDataWFS:LIST_LOCAL_GOVERNMENT_AREAS': 2 features available, geometry column 'SHAPE'
#>   2 features returned

Spatial operations with geos

wk geometry works directly with geos:

library(geos)

# Area of each parcel
parcels |>
  mutate(area_m2 = geos_area(geometry))
#> # A tibble: 20 × 19
#>        FID gml_id    OBJECTID    CID VOLUME FOLIO    PID POT_PID LPI   CAD_TYPE1
#>    <int64> <chr>        <int>  <int> <chr>  <int>  <int>   <int> <chr> <chr>    
#>  1    2688 LIST_Cad…     2688 1.10e6 59944      4 5.68e6       0 GRF64 Private …
#>  2    2877 LIST_Cad…     2877 1.02e6 128641     1 1.82e6       0 FMR21 Private …
#>  3    3425 LIST_Cad…     3425 1.15e6 126957     1 3.27e6       0 FJA84 Authorit…
#>  4    4682 LIST_Cad…     4682 1.02e6 128641     2 1.82e6       0 FMR21 Private …
#>  5    6102 LIST_Cad…     6102 1.25e6 252507     1 7.57e6       0 FAC83 Private …
#>  6    6287 LIST_Cad…     6287 1.02e6 127656     3 1.79e6       0 FBF18 Private …
#>  7    6395 LIST_Cad…     6395 1.50e6 173638     1 9.17e6       0 GAV70 Private …
#>  8    8685 LIST_Cad…     8685 1.41e6 231575     1 0            0 NA    Casement 
#>  9   10819 LIST_Cad…    10819 1.33e6 152956     1 2.82e6       0 GDH05 Private …
#> 10   11600 LIST_Cad…    11600 1.09e6 58505      2 5.68e6       0 HSF79 Private …
#> 11   13561 LIST_Cad…    13561 1.02e6 127656     2 1.79e6       0 FBF18 Private …
#> 12   14096 LIST_Cad…    14096 1.02e6 127656     1 1.79e6       0 FBF18 Private …
#> 13   16589 LIST_Cad…    16589 1.15e6 NA        NA 0            0 NA    Casement 
#> 14   17354 LIST_Cad…    17354 1.10e6 59944      2 5.68e6       0 GRF64 Private …
#> 15   18846 LIST_Cad…    18846 1.01e6 60155     27 5.57e6       0 GAV68 Private …
#> 16   19312 LIST_Cad…    19312 1.01e6 226698     6 5.56e6       0 GAV57 Private …
#> 17   19408 LIST_Cad…    19408 1.01e6 7826       1 1.48e6       0 GAW20 Private …
#> 18   20464 LIST_Cad…    20464 1.01e6 60754      5 5.60e6       0 GAV95 Private …
#> 19   20516 LIST_Cad…    20516 1.01e6 59722     72 5.58e6       0 GAS73 Private …
#> 20   20534 LIST_Cad…    20534 1.01e6 218685     1 5.57e6       0 GAM93 Private …
#> # ℹ 9 more variables: CAD_TYPE2 <chr>, TENURE_TY <chr>, FEAT_NAME <chr>,
#> #   STRATA_LEV <chr>, COMP_AREA <dbl>, MEAS_AREA <dbl>, UFI <chr>,
#> #   geometry <wk_wkb>, area_m2 <dbl>

Working with different service types

readwfs auto-detects the service type from the URL. The same wfs_read() interface works for all of them:

# WFS (auto-detected)
wfs_layers("https://example.com/arcgis/services/Data/MapServer/WFSServer")
#> NULL

# OAPIF (auto-detected from /collections in URL)
wfs_layers("https://example.com/ogc/features/collections")
#> NULL

# ArcGIS REST (auto-detected from arcgis + MapServer)
wfs_layers("https://example.com/arcgis/rest/services/Data/MapServer",
           driver = "ESRIJSON")
#> NULL

You can also force the driver explicitly with driver = "WFS", driver = "ESRIJSON", or driver = "OAPIF".