Get ASGS geometry from ABS ArcGIS REST API #experiments

query.r

library(sf)
library(httr2)
library(ggplot2)

# based on https://community.esri.com/t5/gis-blog/
# accessing-arcgis-rest-services-using-r/ba-p/898451

api_query <- url_parse("https://geo.abs.gov.au")

api_root <- "/arcgis/rest/services"
service <- "ASGS2021"               # ASGS[2011:2022], some others
product <- "STE"                    # geography code
layer <- "0"                        # full layer
suffix <- "MapServer/0/query"       # extra stuff

api_query$path <- paste(api_root, service, product, "MapServer", layer, "query", sep = "/")
api_query$query <- list(
  where = "STATE_NAME_2021='Victoria'",
  outFields = "*",
  returnGeometry = "true",
  f = "geojson")

# other request options:
# 1) `where` can be any legal SQL WHERE clause. for example:
#   where = "STATE_NAME_2021 IN ('Victoria', 'Tasmania')"
#   where = "STATE_CODE_2021>'0'"
# 2) `geometry` can get stuff within (or touching) a bounding box
# full docs: https://geo.abs.gov.au/arcgis/sdk/rest/index.html#/Query_Map_Service_Layer/02ss0000000r000000

# get the boundary
api_query |>
  url_build() |>
  request() |>
  req_perform() ->
api_request

# push the boundary into sf
api_request |>
  resp_body_string() |>
  st_read() ->
boundaries

ggplot() + geom_sf(data = boundaries)

@Deguerre You're totally right! Can't believe I missed that. These are working in the browser:

• https://geo.abs.gov.au/arcgis/rest/services/ASGS2021/STE/MapServer/0/query?where=STATE_NAME_2021%3D%27Victoria%27&outFields=%2A&returnGeometry=true&f=geojson
• https://geo.abs.gov.au/arcgis/rest/services/ASGS2021/STE/MapServer/0/query?where=STATE_CODE_2021%3D%272%27&outFields=%2A&returnGeometry=true&f=geojson

Not in the above snippet, though; I'm still getting a 403 Forbidden. I'll have to use some of httr's debugging tools to see whether it's a user agent issue or something else. Thanks very much!

I'd love to modify it to return all features too, but STATE_NAME_2021='* (escaped) or just removing the where clause entirely seems to return an empty FeatureCollection.

EDIT: Aha, going to just use an inequality operator to get them all (although clearly I have to be careful about using inequality on strings!): STATE_CODE_2021>'0'

And there we go:

# get the boundary
api_query |>
  url_build() |>
  request() |>
  req_perform() ->
api_request

# push the boundary into sf
api_request |>
  resp_body_string() |>
  st_read() ->
vic_boundary

ggplot() + geom_sf(data = vic_boundary)

Of course, it'd be even easier just to download the shape files and convert them in something like QGIS.

Well done nonetheless!

Haha, indeed. I'd prefer not to commit massive shapefiles to version control since it's for a pipeline (although I forgot that {strayr} and {absmapsdata} can download shapefiles too; that probably would've been a better solution in hindsight).