ArcGIS is a set of powerful propriety geospatial tools developed by
ESRI. It is widely used in the industry and academia. In this blog, I
will explore the arcgis R package by getting data from the
Australian tree crop map database, hosted on ArcGIS.
The Australian tree crop map is a database of tree crops in Australia. The map is the work of the Applied Agricultural Remote Sensing Centre at the University of New England.
The dashboard is hosted online at the following link: https://experience.arcgis.com/experience/6cde8c0467e542398fb0afd1dde48a73/
The data informing the dashboard is hosted on ArcGIS and can be accessed at the following link: https://www.arcgis.com/home/item.html?id=12bdb57eec7548dd90409603645cfd5e
Full documentation can be found at the following link on installation and setup. This is hosted on R-universe.
arcgis is a meta-package for arcgisutils,
arcgislayers, arcgisgeocode and
arcgisplaces. For me successful installation required the
indication of the location on the r-universe
repository.
install.packages("arcgis", repos = c("https://r-arcgis.r-universe.dev", "https://cloud.r-project.org"))Next lets load the library into the R environment, and the simple features
(sf) library.
library(arcgis)## Attaching core arcgis packages:
## → arcgisutils v0.3.1.9000
## → arcgislayers v0.3.1.9000
## → arcgisgeocode v0.2.2
## → arcgisplaces v0.1.1library(sf) # also needed for working with spatial data## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUElibrary(leaflet) # for plotting onto interactive maps
There are two feature layers as part of the Australian tree crop map database. A Points database and a polygon database, we will try querying each of these for the northern rivers area of New South Wales, Australia.
We set the point URL using the arc_open function, which stores the database location and some metadata about the feature layer.
atcm_pts <- arc_open("https://services5.arcgis.com/3foZbDxfCo9kcPwP/ArcGIS/rest/services/AARSC_treecrop_point/FeatureServer/0")
atcm_pts## <FeatureLayer>
## Name: atcm_20241024_points
## Geometry Type: esriGeometryPoint
## CRS: 7844
## Capabilities: QueryThis feature server contains point geometry, a coordinate reference system (CRS) 7844 and is capable of being queried.
We can obtain more metadata information from the feature server using
the list_fields() function.
list_fields(atcm_pts)## name type alias sqlType nullable editable domain
## 1 OBJECTID esriFieldTypeOID OBJECTID sqlTypeOther FALSE FALSE NA
## 2 commodity esriFieldTypeString commodity sqlTypeOther TRUE TRUE NA
## 3 source esriFieldTypeString source sqlTypeOther TRUE TRUE NA
## 4 year esriFieldTypeString year sqlTypeOther TRUE TRUE NA
## 5 hectares esriFieldTypeDouble hectares sqlTypeOther TRUE TRUE NA
## defaultValue length
## 1 NA NA
## 2 NA 255
## 3 NA 255
## 4 NA 255
## 5 NA NAThe data can be queried with the arc_select function,
which is effectively a wrapper for the api and uses syntax similar to
querying a SQL database.
arc_select(atcm_pts, where = c("commodity = 'avocado'"))## Registered S3 method overwritten by 'jsonify':
## method from
## print.json jsonlite## Simple feature collection with 2166 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 113.6694 ymin: -41.18715 xmax: 153.5533 ymax: -15.7996
## Geodetic CRS: GDA2020
## First 10 features:
## OBJECTID commodity source year hectares geometry
## 1 1 avocado Imagery 2021 22.962312 POINT (145.1045 -17.16699)
## 2 2 avocado Imagery 2021 1.241645 POINT (145.3262 -16.24474)
## 3 3 avocado Imagery 2021 0.495399 POINT (145.964 -17.60302)
## 4 4 avocado Imagery 2021 1.546013 POINT (145.3421 -17.02106)
## 5 5 avocado Imagery 2021 0.739951 POINT (145.342 -17.01513)
## 6 6 avocado Imagery 2021 4.507486 POINT (145.3132 -17.16463)
## 7 7 avocado Imagery 2021 0.962824 POINT (145.3401 -17.02179)
## 8 8 avocado Imagery 2021 10.928766 POINT (145.2054 -17.19938)
## 9 9 avocado Imagery 2021 3.613079 POINT (145.3156 -17.16684)
## 10 10 avocado Field 2022 9.129086 POINT (152.1197 -24.99829)The package functions returnsf simple features
objects.
From a 10km radius around a point. I will choose the centre of Alstonville in Northern New South Wales. This area is known to contain lots of Avocado and Macadamia farms.
# Define a point (e.g., coordinates for a location)
point_coords <- c(153.4404,-28.8421) # Example coordinates (longitude,latitude)
# Create a sf spatial point feature
# We will use the crs 7844 specified in the ATCM metadata
alston <- st_sfc(st_point(point_coords), crs = 7844)
# Buffer the point to create a polygon (e.g., 10 kilometres)
buffer_alston <- st_buffer(alston, 10000)
# transform the data to WGS84 so we can plot it on a leaflet map
buffer_alston <- st_transform(buffer_alston, 4326)
# plot onto a interactive leaflet map
leaflet(buffer_alston) |>
addProviderTiles(provider = "Esri.WorldImagery")|>
addMarkers(lng = point_coords[1], lat = point_coords[2],
popup = "Alstonville")|> # add the centre point as a marker
addPolygons(color = "green", fillOpacity = 0.2) # add the polygon buffer layerNow lets query and download only the data from inside the buffer
polygon. This is done as before with the arc_select
function. The filter_geom argument is used to specify
sf class objects such as bbox,
sfc or sfg to filter the database on the
arcgis server side. This limits the data download to only what is
required by the user. We also need to specify how this layer should be
queried with the predicate argument. We will specify
intersects as the predicate.
alston_farms <-
arc_select(atcm_pts,
filter_geom = buffer_alston,
predicate = "intersects")
head(alston_farms)## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 153.4086 ymin: -28.8758 xmax: 153.4478 ymax: -28.83944
## Geodetic CRS: GDA2020
## OBJECTID commodity source year hectares geometry
## 1 75 avocado Field 2022 1.712882 POINT (153.4445 -28.8758)
## 2 76 avocado Field 2022 17.237826 POINT (153.44 -28.85337)
## 3 80 avocado Field 2022 2.003204 POINT (153.4476 -28.8748)
## 4 278 avocado Imagery 2022 0.796033 POINT (153.4478 -28.87576)
## 5 279 avocado Imagery 2022 0.551528 POINT (153.4312 -28.86598)
## 6 280 avocado Imagery 2022 0.412121 POINT (153.4086 -28.83944)Next we can plot this on an interactive leaflet map with the WGS84 coordinate reference system.
# transform the data to WGS84 so we can plot it on a leaflet map
alston_farms <- st_transform(alston_farms, 4326)
leaflet(alston_farms) |>
addProviderTiles(provider = "Esri.WorldImagery")|>
addCircleMarkers(color = "orange",
weight = 1, # circle marker edge width
fillOpacity = 0.8,
radius = 4,
popup = alston_farms$commodity)|>
addMarkers(lng = point_coords[1], lat = point_coords[2],
popup = "Alstonville")
Lets do this again for the polygon feature layer.
atcm_poly <- arc_open("https://services5.arcgis.com/3foZbDxfCo9kcPwP/ArcGIS/rest/services/AARSC_treecrop_polygon/FeatureServer/0")
list_fields(atcm_poly)## name type alias sqlType nullable
## 1 OBJECTID esriFieldTypeOID OBJECTID sqlTypeOther FALSE
## 2 commodity esriFieldTypeString commodity sqlTypeOther TRUE
## 3 source esriFieldTypeString source sqlTypeOther TRUE
## 4 year esriFieldTypeString year sqlTypeOther TRUE
## 5 hectares esriFieldTypeDouble hectares sqlTypeOther TRUE
## 6 Shape__Area esriFieldTypeDouble Shape__Area sqlTypeDouble TRUE
## 7 Shape__Length esriFieldTypeDouble Shape__Length sqlTypeDouble TRUE
## editable domain defaultValue length
## 1 FALSE NA NA NA
## 2 TRUE NA NA 255
## 3 TRUE NA NA 255
## 4 TRUE NA NA 255
## 5 TRUE NA NA NA
## 6 FALSE NA NA NA
## 7 FALSE NA NA NAFor a different method we will use a bounding box to query and the polygon layer.
# lets get the data from inside a bounding box of the buffer polygon
buffer_box <-
st_as_sfc(st_bbox(buffer_alston))
# request only the data in the box
alston_farms <-
arc_select(atcm_poly,
filter_geom = buffer_box,
predicate = "intersects")
# transform the data to WGS84 so we can plot it on a leaflet map
alston_farms <- st_transform(alston_farms, 4326)
leaflet(alston_farms) |>
addProviderTiles(provider = "Esri.WorldImagery")|>
addPolygons(color = "orange",
weight = 1,
fillOpacity = 0.8,
popup = alston_farms$commodity)|>
addPolygons(data = buffer_box, color = "green", weight = 2, fillOpacity = 0.2)I hope this helps your spatial data analysis in R. I know I will find this a useful reference for my work into the future. Feel free to reach out and let me know how you are using this package in different interesting ways.