library(VicmapR)
#> 
#> Attaching package: 'VicmapR'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(sf)
#> Linking to GEOS 3.10.0, GDAL 3.3.2, PROJ 8.1.1
library(leaflet)


#check sf installation
sf::sf_extSoftVersion()
#>           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
#>       "3.10.0"        "3.3.2"        "8.1.1"         "true"         "true" 
#>           PROJ 
#>        "8.1.1"

How to query data

Searching for data

In order to begin a query of the WFS server a spatial layer must be selected. To know which layers are available use the listLayers() function, which will return ~ 690 layers to choose from.

available_layers <- listLayers()

head(available_layers, 10)
#>                                         Name
#> 1     datavic:WATER_ISC2010_BANKFULL_WIDTH_R
#> 2          datavic:WATER_ISC2010_BARE_GROUND
#> 3    datavic:WATER_ISC2010_CHANNEL_TRANSECTS
#> 4        datavic:WATER_ISC2010_FRAGMENTATION
#> 5          datavic:WATER_ISC2010_LARGE_TREES
#> 6          datavic:WATER_ISC2010_TOE_OF_BANK
#> 7          datavic:WATER_ISC2010_TOP_OF_BANK
#> 8  datavic:WATER_ISC2010_VEGETATION_OVERHANG
#> 9     datavic:WATER_ISC2010_VEGETATION_WIDTH
#> 10        datavic:WATER_ISC2010_WATER_BODIES
#>                                                                      Title
#> 1  2010 Index of Stream Condition - Bank Full Width Reach polygon features
#> 2            2010 Index of Stream Condition - Bare Ground polygon features
#> 3          2010 Index of Stream Condition - Channel Transect line features
#> 4          2010 Index of Stream Condition - Fragmentation polygon features
#> 5            2010 Index of Stream Condition - Large Trees polygon features
#> 6               2010 Index of Stream Condition - Toe of Bank line features
#> 7               2010 Index of Stream Condition - Top of Bank line features
#> 8    2010 Index of Stream Condition - Vegetation Overhang polygon features
#> 9          2010 Index of Stream Condition - Vegetation Width line features
#> 10            2010 Index of Stream Condition - Water Body polygon features

Vicmap promise

VicmapR introduces a new class called vicmap_promise, which is an extension to the httr::url class. Essentially this object is how the vicmap query is stored before data is collected. That is to say vicmap_promise is essentially a promise of what data will be retrieved.

In order to generate a new promise the vicmap_query function can be used to select the layer. The promise prints a sample of the data (max = 6 rows) as well as the dimensions (nrow and ncol).

# query the watercourse layer
vicmap_query(layer = "datavic:VMHYDRO_WATERCOURSE_DRAIN")
#> • Using collect() on this object will return 191506 features and 16
#> • fields
#> • At most six rows of the record are printed here
#> ────────────────────────────────────────────────────────────────────────────────
#> Simple feature collection with 6 features and 15 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 142.7675 ymin: -35.06905 xmax: 143.324 ymax: -35.04559
#> Geodetic CRS:  GDA94
#> # A tibble: 6 × 16
#>   id                  PFI    UFI FEATURE_TYPE_CODE NAME  NAMED_FEATURE_ID ORIGIN
#>   <chr>             <int>  <int> <chr>             <chr> <chr>            <chr> 
#> 1 VMHYDRO_WATERC… 8553127 2.55e6 watercourse_chan… <NA>  <NA>             2     
#> 2 VMHYDRO_WATERC… 8553130 2.55e6 watercourse_chan… <NA>  <NA>             2     
#> 3 VMHYDRO_WATERC… 8553143 2.55e6 watercourse_chan… <NA>  <NA>             2     
#> 4 VMHYDRO_WATERC… 8553149 2.55e6 watercourse_chan… <NA>  <NA>             2     
#> 5 VMHYDRO_WATERC… 8553158 2.55e6 watercourse_chan… <NA>  <NA>             2     
#> 6 VMHYDRO_WATERC… 8553168 2.55e6 watercourse_chan… <NA>  <NA>             2     
#> # … with 9 more variables: CONSTRUCTION <chr>, USAGE <chr>, HIERARCHY <chr>,
#> #   FEATURE_QUALITY_ID <int>, CREATE_DATE_PFI <dttm>, SUPERCEDED_PFI <chr>,
#> #   CREATE_DATE_UFI <dttm>, OBJECTID <int>, geometry <LINESTRING [°]>

Adding arguments to the query

The vicmap_promise object can be easily added to through piping in of additional functions (e.g. head(), filter() and select()).

The resulting query can be displayed using the show_query() function, which will list the WFS parameters.

vicmap_query(layer = "datavic:VMHYDRO_WATERCOURSE_DRAIN") %>%
  head(50) %>% #return only 50 rows
  filter(HIERARCHY == "L") %>% # filter the column 'HIERACHY' to values of 'L'
  select(HIERARCHY, PFI) %>% # select columns 'HIERARCHY' and 'PFI'
  show_query()
#> Note: method with signature 'DBIConnection#character' chosen for function 'dbQuoteIdentifier',
#>  target signature 'wfsConnection#ident'.
#>  "wfsConnection#ANY" would also be valid
#> Note: method with signature 'DBIConnection#character' chosen for function 'dbQuoteIdentifier',
#>  target signature 'wfsConnection#character'.
#>  "wfsConnection#ANY" would also be valid
#> Note: method with signature 'DBIConnection#character' chosen for function 'dbQuoteString',
#>  target signature 'wfsConnection#character'.
#>  "wfsConnection#ANY" would also be valid
#> <base url>
#> http://services.land.vic.gov.au/catalogue/publicproxy/guest/dv_geoserver/wfs
#> 
#> <body>
#> service: wfs 
#> version: 2.0.0 
#> request: GetFeature 
#> typeNames: datavic:VMHYDRO_WATERCOURSE_DRAIN 
#> outputFormat: application/json 
#> count: 50 
#> srsName: EPSG:4283 
#> CQL_FILTER: ("HIERARCHY" = 'L') 
#> propertyName: SHAPE,HIERARCHY,PFI 
#> 
#> <full query url>
#> http://services.land.vic.gov.au/catalogue/publicproxy/guest/dv_geoserver/wfs?service=wfs&version=2.0.0&request=GetFeature&typeNames=datavic%3AVMHYDRO_WATERCOURSE_DRAIN&outputFormat=application%2Fjson&count=50&srsName=EPSG%3A4283&CQL_FILTER=%28%22HIERARCHY%22%20%3D%20%27L%27%29&propertyName=SHAPE%2CHIERARCHY%2CPFI

In order to return a spatial data.frame object (sf) collect() must be used.

watercourse_data <- vicmap_query(layer = "datavic:VMHYDRO_WATERCOURSE_DRAIN") %>%
  head(50) %>% #return only 50 rows
  filter(HIERARCHY == "L") %>% # filter the column 'HIERACHY' to values of 'L'
  select(HIERARCHY, PFI) %>% # select columns 'HIERARCHY' and 'PFI'
  collect()

str(watercourse_data)
#> sf [50 × 6] (S3: sf/tbl_df/tbl/data.frame)
#>  $ id       : chr [1:50] "VMHYDRO_WATERCOURSE_DRAIN.fid--2e4b3d2f_17d0210b7b9_-2d47" "VMHYDRO_WATERCOURSE_DRAIN.fid--2e4b3d2f_17d0210b7b9_-2d46" "VMHYDRO_WATERCOURSE_DRAIN.fid--2e4b3d2f_17d0210b7b9_-2d45" "VMHYDRO_WATERCOURSE_DRAIN.fid--2e4b3d2f_17d0210b7b9_-2d44" ...
#>  $ PFI      : int [1:50] 8553127 8553130 8553143 8553149 8553158 8553168 8553183 8553184 8550948 8550950 ...
#>  $ UFI      : int [1:50] 2552171 2552174 2552187 2552193 2552202 2552212 2552227 2552228 2549992 2549994 ...
#>  $ HIERARCHY: chr [1:50] "L" "L" "L" "L" ...
#>  $ OBJECTID : int [1:50] 2973 2975 2978 2981 2987 2992 2996 2997 2021 2022 ...
#>  $ geometry :sfc_LINESTRING of length 50; first list element:  'XY' num [1:8, 1:2] 143 143 143 143 143 ...
#>  - attr(*, "sf_column")= chr "geometry"
#>  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
#>   ..- attr(*, "names")= chr [1:5] "id" "PFI" "UFI" "HIERARCHY" ...

Geometric filters

VicmapR translates numerous geometric filter functions available in the Victorian Government’s WFS Geoserver supports numerous geometric filters:

  • EQUALS
  • DISJOINT
  • INTERSECTS
  • TOUCHES
  • CROSSES
  • WITHIN
  • CONTAINS
  • OVERLAPS
  • DWITHIN
  • BEYOND
  • BBOX

These filters can be used within the filter() function by providing them an object of class sf/sfc/sfg/bbox. Below is a leaflet map with the melbourne rail network being read in with the use of three different types of filter functions: INTERSECTS(), BBOX() and DWITHIN().

#### Return objects that intersect melbourne ####
# Read in an example shape to restrict our query to using geometric filtering
melbourne <- sf::st_read(system.file("shapes/melbourne.geojson", package="VicmapR"), quiet = F) %>% 
  sf::st_transform(4283)
#> Reading layer `melbourne' from data source 
#>   `/Users/runner/work/_temp/Library/VicmapR/shapes/melbourne.geojson' 
#>   using driver `GeoJSON'
#> Simple feature collection with 1 feature and 8 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 144.897 ymin: -37.85066 xmax: 144.9913 ymax: -37.77545
#> Geodetic CRS:  WGS 84

# Return data that intersects melbourne
rail_intersects <- vicmap_query(layer = "datavic:VMTRANS_TR_RAIL") %>% # layer to query
  filter(INTERSECTS(melbourne)) %>% # more advanced geometric filter
  collect()

rail_bbox <- vicmap_query(layer = "datavic:VMTRANS_TR_RAIL") %>%
  filter(BBOX(sf::st_bbox(melbourne))) %>%
  collect()

rail_dwithin <- vicmap_query(layer = "datavic:VMTRANS_TR_RAIL") %>%
  filter(DWITHIN(melbourne %>% sf::st_centroid(), distance = 10000, units = "meters")) %>%
  collect()

leaflet(width = "100%") %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(data = melbourne, color = "grey", group = "Melbourne polygon") %>%
  addPolygons(data = sf::st_bbox(melbourne) %>% st_as_sfc(), color = "black", group = "Melbourne bbox") %>%
  addPolylines(data = rail_intersects, color = "Red", group = "INTERSECTS") %>% 
  addPolylines(data = rail_bbox, color = "Blue", group = "BBOX") %>%
  addPolylines(data = rail_dwithin, color = "Green", group = "DWITHIN") %>%
  addLayersControl(baseGroups = c("Melbourne polygon", "Melbourne bbox"), 
                   overlayGroups = c("INTERSECTS", "BBOX", "DWITHIN")) %>%
  hideGroup(c("BBOX", "DWITHIN"))