In arcpy: Interface to 'ArcGIS' 'Python' Modules

The arcpy module provides classes and methods a variety of tasks, including manipulating geometries. These classes and methods can be accessed in R as well using functionality provided by the reticulate package.

This document demonstrates an example task where the geometry of polygons in a shapefile is shifted 5000 linear units in the y-direction and 2500 linear units in the x-direction. This contrived task demonstrates an approach to manipulating the arcobjects class instances.

library(arcpy)

# set workspace
arcpy$env$workspace = tempdir()
arcpy$env$scratchWorkspace = tempdir()

# copy feature for example
fc = arcpy$management$CopyFeatures(system.file("CA_Counties",
  "CA_Counties_TIGER2016.shp", package = "arcpy"), "CA_Counties")

# get the width and height of the layer
arcpy$Describe(fc)$extent$XMin

## [1] -13857275

arcpy$Describe(fc)$extent$YMin

## [1] 3832931

The reticulate package provides functions for manipulating python objects. These tools include iterators and the %as% operator for aliasing. The code below uses %as% to create an instance of the fc layer, which we can use to iterate over the polygons.

# create iterable instance of the SHAPE@ objects in the layer
with(arcpy$da$UpdateCursor(fc, "SHAPE@") %as% rows, {
  # move to first row
  row = iter_next(rows)
  # iterate over the features
  while (!is.null(row)) {
    # create an arcpy array for storing features
    arr_feat = arcpy$Array()
    # get the parts of the current feature
    feat = row[[1]]$getPart()
    # get the first part of the current feature
    part = iter_next(feat)
    # iterate over the parts of this feature
    while (!is.null(part)) {
      # create another arcpy array to store feature parts
      arr_part = arcpy$Array()
      # get the first point of this part
      pnt = iter_next(part)
      # iterate over the points in this part
      while (!is.null(pnt)) {
        # get the point X coordinate and subtract 5000
        x = pnt$X + 2500
        # get the point Y coordinate and add 5000
        y = pnt$Y + 5000
        # create a new point
        arr_part$add(arcpy$Point(x, y))
        # iterate to next point
        pnt = iter_next(part)
      }
      # add the modified part to the feature array
      arr_feat$add(arr_part)
      # iterate to next part
      part = iter_next(feat)
    }
    # create a new polygon from the feature array
    polygon = arcpy$Polygon(arr_feat)
    # overwrite the original polygon with the new polygon
    row[[1]] = polygon
    # update the cursor
    rows$updateRow(row)
    # iterate to next feature
    row = iter_next(rows)
  }
})

# recalculate the width and height of the layer
arcpy$management$RecalculateFeatureClassExtent(paste(fc))
arcpy$Describe(fc)$extent$XMin
arcpy$Describe(fc)$extent$YMin

## <Result ''>
## [1] -13854775
## [1] 3837931

It's also possible to mix and match iterators and R functions, often with significantly faster results. Here's an alternative approach that makes judicious use of reticulate's iterate() function in combination with lapply() and the da_read() and da_update() functions:

translate_geom = function(g, dx, dy) {
  arr_feat = arcpy$Array()
  feat = g$getPart()
  parts = iterate(feat, function(part) {
    arr_part = arcpy$Array()
    pnts = iterate(part, function(pnt) {
      x = pnt$X + dx
      y = pnt$Y + dy
      arcpy$Point(x, y)
    })
    lapply(pnts, function(part) arr_part$add(part))
    arr_part
  })
  lapply(parts, function(feat) arr_feat$add(feat))
  arcpy$Polygon(arr_feat)
}

fc.d = da_read(fc, "SHAPE@")
fc.d[["SHAPE@"]] = lapply(fc.d[["SHAPE@"]][1], translate_geom,
  dx = 2500, dy = 5000)

da_update(fc, fc.d)

arcpy$management$RecalculateFeatureClassExtent(fc)
arcpy$Describe(fc)$extent$XMin
arcpy$Describe(fc)$extent$YMin

## <Result ''>
## [1] -13471139
## [1] 4787916

Handling python objects can be challenging in some cases, but reticulate provides virtually all the tools needed to manipulate arcpy classes.

arcpy documentation built on May 29, 2024, 9:58 a.m.