Skip to content

Instantly share code, notes, and snippets.

@jimjam-slam
Last active November 2, 2022 21:53
Show Gist options
  • Save jimjam-slam/82730dd84589bda8674f0e7145195eb6 to your computer and use it in GitHub Desktop.
Save jimjam-slam/82730dd84589bda8674f0e7145195eb6 to your computer and use it in GitHub Desktop.
Get ASGS geometry from ABS ArcGIS REST API #experiments
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
Copy link

Deguerre commented Nov 1, 2022

You want either STATE_CODE_2021=2 or STATE_NAME_2021='Victoria' and don't forget the quotes.

@jimjam-slam
Copy link
Author

jimjam-slam commented Nov 1, 2022

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

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!

@jimjam-slam
Copy link
Author

jimjam-slam commented Nov 2, 2022

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'

@jimjam-slam
Copy link
Author

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)

image

@Deguerre
Copy link

Deguerre commented Nov 2, 2022

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

@Deguerre
Copy link

Deguerre commented Nov 2, 2022

Well done nonetheless!

@jimjam-slam
Copy link
Author

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).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment