Skip to content

Latest commit

 

History

History
274 lines (193 loc) · 11.7 KB

README.md

File metadata and controls

274 lines (193 loc) · 11.7 KB

esri2sf

Scraping Geographic Features from ArcGIS Server

Still many geographic data is delivered through ESRI’s ArcGIS Server. It is not easy to utilize the geographic data in the GIS servers from data analysis platform like R or Pandas. This package enables users to scrape vector data in ArcGIS Server from R through the server’s REST API. It download geographic features from ArcGIS Server and saves it as Simple Features.

How it works

This program sends a request to an ArcGIS Server and gets json responses containing coordinates of geometries of which format is not the same as geojson. So it converts the json into simple feature geometries from the response. Then it combines attribute data to the geometries to create sf dataframe. Often ArcGIS servers limits the maximum number of rows in the result set. So this program creates 500 features per request and automatically re-send requests until it gets all features in the dataset.

Install

Use remotes to install this package. This package has dependencies on dplyr, sf, httr, jsonlite, rstudioapi, DBI, RSQLite, crayon and rmarkdown, knitr, and pbapply are suggested.

library(remotes)
install_github("yonghah/esri2sf")

How to use

What you need is the URL of REST service you want. You can get the URL by viewing the URL widget on the service’s webpage (see image below), by asking a GIS admin, or looking at the javascript code of a webpage where it creates a feature layer.

REST Service screenshot

Point data

library("esri2sf")
url <- "https://services.arcgis.com/V6ZHFr6zdgNZuVG0/arcgis/rest/services/Landscape_Trees/FeatureServer/0"
df <- esri2sf(url)

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPoint

## Service Coordinate Reference System: 3857

## Output Coordinate Reference System: 4326

plot(df)

## Warning: plotting the first 9 out of 56 attributes; use max.plot = 56 to plot
## all

Polyline data

You can filter output fields. This may take a minute since it gets 18000 polylines.

url <- "https://services.arcgis.com/V6ZHFr6zdgNZuVG0/arcgis/rest/services/Florida_Annual_Average_Daily_Traffic/FeatureServer/0"
df <- esri2sf(url, outFields=c("AADT", "DFLG"))

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolyline

## Service Coordinate Reference System: 3857

## Output Coordinate Reference System: 4326

plot(df)

Polygon data

You can filter rows as well by giving a where condition.

url <- "https://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/3"
df <- esri2sf(url, 
              where = "STATE_NAME = 'Michigan'", 
              outFields = c("POP2000", "pop2007", "POP00_SQMI", "POP07_SQMI"))

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolygon

## Service Coordinate Reference System: 4269

## Output Coordinate Reference System: 4326

plot(df)

Tabular data

You can download non-spatial tables of the ‘Table’ layer type using esri2df().

df <- esri2df('https://sampleserver1.arcgisonline.com/ArcGIS/rest/services/WaterTemplate/WaterDistributionInventoryReport/MapServer/5', objectIds = paste(1:50, collapse = ","))

## Layer Type: Table

df

## # A tibble: 50 × 6
##    OBJECTID FACILITYID FCLASS ASSETCOND CONDDATE REPLSCORE
##       <int> <chr>      <chr>  <lgl>     <lgl>        <int>
##  1        1 1          wMain  NA        NA               0
##  2        2 2          wMain  NA        NA              45
##  3        3 3          wMain  NA        NA              45
##  4        4 4          wMain  NA        NA              45
##  5        5 5          wMain  NA        NA              15
##  6        6 6          wMain  NA        NA              15
##  7        7 1          wMain  NA        NA               0
##  8        8 2          wMain  NA        NA              45
##  9        9 3          wMain  NA        NA              45
## 10       10 4          wMain  NA        NA              45
## # … with 40 more rows

crs parameter example

When specifying the CRS parameter, any transformation that happens will be done within ESRI’s REST API. Caution should be taken when specifying an output crs that requires a datum transformation as ESRI will automatically apply a default transformation (with no feedback as to which one) which could end up adding small unexpected errors into your data. By default, esri2sf() will transform any datasource to WGS 1984 (EPSG:4326), but it may be safer to set crs = NULL which will return the data in the same CRS as it is being hosted as in the Feature/Map Service. That way you can control any transformation manually in a known, reproducible manner.

url <- "https://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/3"
where <- "STATE_NAME = 'Michigan'"
outFields <- c("POP2000", "pop2007", "POP00_SQMI", "POP07_SQMI")

#default crs = 4326
esri2sf(url, where = where, outFields = outFields) 

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolygon

## Service Coordinate Reference System: 4269

## Output Coordinate Reference System: 4326

## Simple feature collection with 83 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.4082 ymin: 41.69749 xmax: -82.41984 ymax: 48.1738
## Geodetic CRS:  WGS 84
## First 10 features:
##    POP2000 POP2007 POP00_SQMI POP07_SQMI                          geoms
## 1     2301    2324        4.1        4.2 MULTIPOLYGON (((-88.49753 4...
## 2    36016   36791       34.6       35.3 MULTIPOLYGON (((-88.50068 4...
## 3     7818    7444        5.9        5.6 MULTIPOLYGON (((-88.98743 4...
## 4     8746    8760        9.5        9.5 MULTIPOLYGON (((-88.67172 4...
## 5    64634   64904       34.6       34.7 MULTIPOLYGON (((-87.6137 45...
## 6    17370   17057       15.2       14.9 MULTIPOLYGON (((-88.9853 46...
## 7     7024    7159        7.6        7.7 MULTIPOLYGON (((-85.85923 4...
## 8     9862   10188       10.5       10.9 MULTIPOLYGON (((-87.11047 4...
## 9     8903    8781        7.3        7.2 MULTIPOLYGON (((-86.45828 4...
## 10   13138   12750       10.8       10.5 MULTIPOLYGON (((-88.9252 46...

#No transformation (recommended)
esri2sf(url, where = where, outFields = outFields, crs = NULL)

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolygon

## Service Coordinate Reference System: 4269

## Simple feature collection with 83 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.4082 ymin: 41.69749 xmax: -82.41984 ymax: 48.1738
## Geodetic CRS:  NAD83
## First 10 features:
##    POP2000 POP2007 POP00_SQMI POP07_SQMI                          geoms
## 1     2301    2324        4.1        4.2 MULTIPOLYGON (((-88.49753 4...
## 2    36016   36791       34.6       35.3 MULTIPOLYGON (((-88.50068 4...
## 3     7818    7444        5.9        5.6 MULTIPOLYGON (((-88.98743 4...
## 4     8746    8760        9.5        9.5 MULTIPOLYGON (((-88.67172 4...
## 5    64634   64904       34.6       34.7 MULTIPOLYGON (((-87.6137 45...
## 6    17370   17057       15.2       14.9 MULTIPOLYGON (((-88.9853 46...
## 7     7024    7159        7.6        7.7 MULTIPOLYGON (((-85.85923 4...
## 8     9862   10188       10.5       10.9 MULTIPOLYGON (((-87.11047 4...
## 9     8903    8781        7.3        7.2 MULTIPOLYGON (((-86.45828 4...
## 10   13138   12750       10.8       10.5 MULTIPOLYGON (((-88.9252 46...

Also since the addition of the WKT1_ESRI output from sf::st_crs() in sf version 1.0-1, you can enter common CRS format (any that sf::st_crs() can handle) into the crs parameters and it will be able to convert to the ESRI formatted WKT needed for the outSR field in the REST query. Below are examples of the variety of input types that you can use with the crs parameters. All examples are just different formulations of the ESRI:102690 CRS.

#ESRI Authority Code
df1 <- esri2sf(url, where = where, outFields = outFields, crs = "ESRI:102690")

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolygon

## Service Coordinate Reference System: 4269

## Output Coordinate Reference System: ESRI:102690

#PROJ string
df2 <- esri2sf(url, where = where, outFields = outFields, crs = "+proj=lcc +lat_1=42.1 +lat_2=43.66666666666666 +lat_0=41.5 +lon_0=-84.36666666666666 +x_0=4000000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs")

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolygon

## Service Coordinate Reference System: 4269

## Output Coordinate Reference System: +proj=lcc +lat_1=42.1 +lat_2=43.66666666666666 +lat_0=41.5 +lon_0=-84.36666666666666 +x_0=4000000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs

#OGC WKT
df3 <- esri2sf(url, where = where, outFields = outFields, crs = 'PROJCS["NAD_1983_StatePlane_Michigan_South_FIPS_2113_Feet",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",13123333.33333333],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-84.36666666666666],PARAMETER["Standard_Parallel_1",42.1],PARAMETER["Standard_Parallel_2",43.66666666666666],PARAMETER["Latitude_Of_Origin",41.5],UNIT["Foot_US",0.30480060960121924],AUTHORITY["EPSG","102690"]]')

## Layer Type: Feature Layer

## Geometry Type: esriGeometryPolygon

## Service Coordinate Reference System: 4269

## Output Coordinate Reference System: PROJCS["NAD_1983_StatePlane_Michigan_South_FIPS_2113_Feet",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",13123333.33333333],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-84.36666666666666],PARAMETER["Standard_Parallel_1",42.1],PARAMETER["Standard_Parallel_2",43.66666666666666],PARAMETER["Latitude_Of_Origin",41.5],UNIT["Foot_US",0.30480060960121924],AUTHORITY["EPSG","102690"]]

Their similarity on the output CRS can be proven by the following function that calculates the mean difference in X-Y coordinates at each point. All are very close to 0.

coord_diff <- function(df1, df2) {
  suppressWarnings({
    c(
      "x" = mean(sf::st_coordinates(sf::st_cast(df1, "POINT"))[,1] - sf::st_coordinates(sf::st_cast(df2, "POINT"))[,1]),
      "y" = mean(sf::st_coordinates(sf::st_cast(df1, "POINT"))[,2] - sf::st_coordinates(sf::st_cast(df2, "POINT"))[,2])
    )
  })
}
coord_diff(df1, df2)

##             x             y 
##  1.827251e-08 -1.191372e-09

coord_diff(df1, df3)

##             x             y 
##  1.827251e-08 -1.191372e-09

coord_diff(df2, df3)

## x y 
## 0 0