Part 3 - Spatial operations on geometries¶

Previously, in Part 1 of this guide series to arcgis.geometry module, you have seen the introduction to the module and some foundational concepts. In Part 2, you learned how to create geometry objects, their properties, and how to work with one, including its interactions with map widgets. In this part, let's continue to explore spatial operations of geometry objects using two different usage patterns.

Spatial Operations¶

What are spatial operations? Spatial operations are functions that "create new spatial data from specified input data", e.g. union, difference, symmetric difference, and intersection, and can be used for a variety of spatial analyses. For example, government authorities may use the intersect operation to determine whether a proposed road cuts through a restricted piece of land such as a nature reserve or a private property.

Chained operations are even more powerful. For example, in order to identify food deserts within an urban area, the analysis might begin by performing a union operation to the service areas of grocery stores, farmer's markets, and food co-ops. Then, taking the difference between this single geometry of all services areas and that of a polygon delineating a neighborhood would reveal the areas within that neighborhood where access to healthy, whole foods may not exist.

Two patterns of applying spatial operations¶

There are two ways of performing spatial operations

(a) Object Oriented Programming (OOP) pattern - here you create a geometry object first, then call methods off the geometry object, and
(b) calling functions from arcgis.geometry.functions directly without initiating any geometry objects, e.g. the from_geo_coordinate_string or to_geo_coordinate_string methods, which takes spatial data as input, analyzes the data, then produces output data that is the derivative of the analysis performed on the input data.

The major difference between these two is that the OOP pattern uses local geometry engines such as shapely or arcpy, while calling functions from arcgis.geometry.functions will use server-side geometry engine by sending the geometries over to the Geometry Service configured with your web GIS server. If your web GIS is ArcGIS Online, the latter pattern requires credits. Further, the latter pattern is not performant for larger datasets. Users are recommended to only use (b) if they do not have access to either of the local geometry engines. Further more, if you want to process geometries on the server, rather than applying (b), you can also publish the data as feature service and then run analyses using arcgis.feature module.

The rest of this page demonstrates a set of spatial operations performed using both these patterns.

a. OOP Pattern - Uses local geometry engine¶

Before spatial operations are being demonstrated, let us first import necesary libraries, and create a GIS instance that connects to ArcGIS Online.

Note: If you are running this guide in an environment without arcpy, you can comment out the line that imports it and proceed with the guide. The Python API will look for shapely library and use it for geometry operations.

Note: When using Jupyter/ArcGIS Notebooks interface, you can quickly visualize a geometry object by just querying it in a cell. This lets the notebook know to display the geometry as a SVG image. If you are running this guide in an environment with arcpy, then to visualize the geometries as SVG graphics, replace the cells that query the object with geom_object.as_arcpy

In [2]:

from arcgis.gis import GIS
from arcgis.geometry import Polygon, Geometry, Point, Polyline
from arcgis.geocoding import geocode
import arcpy

In [3]:

gis = GIS('home')

a1. Union¶

The first method, <first geom object>.union(<second_geometry>), can be called off from a geometry object to construct the geometry that is the set-theoretic union of the input geometries. You are going to see an example of how to create a union of two polygons (e.g. geometry1 and geometry2):

In [4]:

geom1_json = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100}, 'rings': [[[-8232647.749922129, 4978983.410862541], [-8232389.7749516675, 4978840.091434507], [-8232762.405464557, 4978161.712808477], [-8233001.2711779475, 4978295.477607976], [-8232647.749922129, 4978983.410862541]]]}
geom2_json = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100}, 'rings': [[[-8232619.086036522, 4978994.876241834], [-8232275.11940924, 4979644.590982256], [-8231988.480553171, 4979482.162297151], [-8232380.220323131, 4978822.892928192], [-8232619.086036522, 4978994.876241834]]]}

In [13]:

geom1 = Polygon(geom1_json)

In [14]:

geom2 = Polygon(geom2_json)

In [15]:

geom1

Out[15]:

In [16]:

geom2

Out[16]:

In [17]:

geom_union = geom1.union(geom2)
geom_union

Out[17]:

a2. Difference¶

The, difference(second_geometry), constructs the geometry that is composed only of the region unique to the base geometry but not part of the second geometry. The following illustration shows the results that are inside the source geometry, but not the second geometry.

In [18]:

geom3_json = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100}, 'rings': [[[-8233020.380435019, 4978303.121194171], [-8232303.783294846, 4979640.769189159], [-8232026.699067313, 4979497.449761124], [-8232791.069350163, 4978112.028623458], [-8233020.380435019, 4978303.121194171]]]}

In [21]:

geom3 = Polygon(geom3_json)
geom3

Out[21]:

In [22]:

geom_diff = geom_union.difference(geom3)
geom_diff

Out[22]:

a3. Symmetric difference¶

Similar to difference(), the method symmetric_difference(second_geometry) constructs the geometry that is the union of two geometries minus the intersection of those geometries. As told by the definition, the result of difference() should always be included/contained in the result of a symmetric_difference().

In [23]:

geom_sdiff = geom_union.symmetric_difference(geom3)
geom_sdiff

Out[23]:

a4. intersect() V.S. overlaps()¶

The intersect() constructs a geometry that is the geometric intersection of the two input geometries. Different dimension values can be used to create different shape types. The intersection of two geometries of the same shape type is a geometry containing only the regions of overlap between the original geometries, and its arguments include:

second_geometry: required arcgis.geometry.Geometry.
A second geometry dimension: required Integer. The topological dimension (shape type) of the resulting geometry, which can be:
- 1: A zero-dimensional geometry (point or multipoint).
- 2: A one-dimensional geometry (polyline).
- 4: A two-dimensional geometry (polygon).

In [28]:

# scroll to the top of this notebook where these two geometries are visualized
geom_intersection = geom1.intersect(second_geometry=geom2, dimension=4)
geom_intersection

Out[28]:

While intersect() returns a new geometry object, the overlaps() will return a boolean value which indicates if the intersection of the two geometries has the same shape type as one of the input geometries and is not equivalent to either of the input geometries.

According to the definition, we can tell that overlaps() contains two operations - first, it creates the intersection of first and second geometries; second, it checks if the intersection is None, or it being equivalent to either the first or the second geometry - if so, then overlaps() returns False, or else, returns True.

In [31]:

# True, because the intersection of geom1 and geom2 is partial shape of geom1, and is not
# equivalent to either one
geom1.overlaps(geom2)

Out[31]:

True

In [33]:

# False, because the intersection of geom1 and geom_intersection is geom_intersection, 
# and that is equivalent to the second geometry
geom1.overlaps(geom_intersection)

Out[33]:

False

a5. Equals¶

equals() compares the source geometry with the second to check if they are of the same shape type and define the same set of points in the plane. This is a 2D comparison only; M and Z values are ignored.

In [38]:

geom1.equals(geom1)

Out[38]:

True

In [35]:

geom_intersection.equals(geom_sdiff)

Out[35]:

False

a6. generalize() V.S. buffer()¶

Next, you will see the difference between these two spatial operations:

generalize(max_offset) - Creates a new simplified geometry using a specified maximum offset tolerance (as shown in Figs 1 and 2). The result of the generalize() operation against a polyline object is a new polyline, while that against a polygon object is a new polygon.
buffer(distance) - Constructs a polygon at a specified distance from the geometry. The buffering process merges buffer polygons that overlap. Negative distances greater than one-half the maximum interior width of a polygon result in an empty geometry. Illustration of buffer() operations against different geomtry types can be found in Fig 3, which also shows that results of buffer() are always polygon objects.

Fig 1. Generialize operation against a polyline object (Source: ArcMap Toolbox Help)

Fig 2. Generialize operation against a polygon object (Source: ArcMap Toolset Help)

Fig 3. Buffer operation (Source: ArcGIS Desktop Help)

To illustrate generalization, let us construct a complex polyline object

In [101]:

polyline1 = Polyline({'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'paths': [[[-8235190.176786223, 4978434.716352775], [-8235198.537086192, 4978449.048295577],[-8235190.176786223, 4978458.602924112],[-8235172.261857721, 4978455.019938411],[-8235166.290214886, 4978452.631281278], [-8235154.346929218, 4978452.631281278],[-8235145.9866292495, 4978456.214266978],[-8235145.9866292495, 4978465.768895513], [-8235145.9866292495, 4978471.7405383475],[-8235145.9866292495, 4978480.100838316],[-8235148.375286384, 4978484.878152583], [-8235150.763943518, 4978487.266809717],[-8235160.318572052, 4978488.461138284],[-8235165.09588632, 4978489.65546685], [-8235171.067529154, 4978492.044123984],[-8235178.233500555, 4978492.044123984],[-8235181.816486255, 4978490.849795417], [-8235185.399471956, 4978487.266809717],[-8235197.342757625, 4978478.906509749],[-8235206.89738616, 4978480.100838316],[-8235210.48037186, 4978495.627109685], [-8235210.48037186, 4978508.76472392], [-8235210.48037186, 4978509.959052487], [-8235210.48037186, 4978514.736366754], [-8235210.48037186, 4978520.7080095885], [-8235210.48037186, 4978521.902338156], [-8235210.48037186, 4978523.096666723], [-8235210.48037186, 4978524.29099529], [-8235210.48037186, 4978525.485323856], [-8235210.48037186, 4978529.068309557], [-8235206.89738616, 4978535.039952391], [-8235206.89738616, 4978536.234280958], [-8235205.703057593, 4978542.205923792], [-8235204.508729026, 4978543.400252359], [-8235203.314400459, 4978544.594580926], [-8235202.120071892, 4978545.788909492], [-8235196.148429058, 4978550.56622376], [-8235196.148429058, 4978552.954880893], [-8235196.148429058, 4978554.14920946], [-8235196.148429058, 4978555.343538027], [-8235196.148429058, 4978556.537866594], [-8235194.954100491, 4978557.732195161], [-8235194.954100491, 4978558.926523728], [-8235194.954100491, 4978560.120852294], [-8235193.759771924, 4978561.315180861], [-8235190.176786223, 4978569.67548083], [-8235165.09588632, 4978581.618766498], [-8235165.09588632, 4978592.3677236], [-8235161.512900619, 4978594.756380733], [-8235161.512900619, 4978595.9507093], [-8235156.735586352, 4978599.533695001], [-8235151.958272084, 4978600.728023568], [-8235150.763943518, 4978600.728023568], [-8235147.180957817, 4978603.116680701], [-8235135.237672148, 4978607.893994969], [-8235134.043343581, 4978607.893994969], [-8235132.849015014, 4978607.893994969], [-8235132.849015014, 4978609.088323535], [-8235131.654686447, 4978609.088323535], [-8235130.46035788, 4978609.088323535], [-8235130.46035788, 4978610.282652102], [-8235125.683043613, 4978611.476980669], [-8235105.379457977, 4978611.476980669], [-8235101.796472277, 4978611.476980669], [-8235100.60214371, 4978611.476980669], [-8235099.407815143, 4978611.476980669], [-8235098.2134865755, 4978611.476980669], [-8235097.0191580085, 4978611.476980669], [-8235095.8248294415, 4978611.476980669], [-8235087.464529474, 4978611.476980669], [-8235086.270200907, 4978611.476980669], [-8235085.07587234, 4978611.476980669], [-8235083.881543773, 4978611.476980669], [-8235082.687215206, 4978611.476980669], [-8235081.492886639, 4978611.476980669], [-8235075.521243805, 4978611.476980669], [-8235074.326915238, 4978611.476980669], [-8235063.577958137, 4978607.893994969], [-8235059.994972437, 4978605.505337835], [-8235051.6346724685, 4978591.173395033], [-8235051.6346724685, 4978589.979066466], [-8235050.4403439015, 4978587.590409333], [-8235034.914072532, 4978579.230109365], [-8235027.748101131, 4978558.926523728], [-8235026.553772564, 4978558.926523728], [-8235026.553772564, 4978555.343538027], [-8235028.942429698, 4978555.343538027], [-8235031.331086832, 4978555.343538027], [-8235045.663029634, 4978555.343538027], [-8235051.6346724685, 4978555.343538027], [-8235052.829001036, 4978544.594580926], [-8235052.829001036, 4978543.400252359], [-8235052.829001036, 4978531.45696669], [-8235052.829001036, 4978530.262638124], [-8235052.829001036, 4978526.679652423], [-8235052.829001036, 4978520.7080095885], [-8235049.2460153345, 4978520.7080095885], [-8235033.719743965, 4978520.7080095885], [-8235021.776458297, 4978520.7080095885], [-8235021.776458297, 4978521.902338156], [-8235021.776458297, 4978531.45696669], [-8235016.99914403, 4978537.428609525], [-8235012.221829762, 4978537.428609525], [-8235006.250186928, 4978537.428609525], [-8234994.306901259, 4978530.262638124], [-8234993.112572692, 4978529.068309557], [-8234991.918244125, 4978527.87398099], [-8234991.918244125, 4978526.679652423], [-8234989.529586992, 4978523.096666723], [-8234981.169287024, 4978514.736366754], [-8234981.169287024, 4978513.542038187], [-8234972.808987056, 4978499.210095385], [-8234979.974958457, 4978503.987409652], [-8234962.060029954, 4978489.65546685], [-8234960.865701388, 4978475.323524048], [-8234958.477044254, 4978470.5462097805], [-8234958.477044254, 4978465.768895513], [-8234958.477044254, 4978464.574566946], [-8234958.477044254, 4978463.380238379], [-8234958.477044254, 4978462.185909812], [-8234958.477044254, 4978452.631281278], [-8234958.477044254, 4978437.105009909], [-8234958.477044254, 4978435.910681342], [-8234964.448687088, 4978421.5787385395], [-8234965.643015655, 4978421.5787385395], [-8234977.586301323, 4978408.441124304], [-8234985.946601291, 4978402.46948147], [-8234988.335258425, 4978401.275152903], [-8234988.335258425, 4978400.080824336], [-8234989.529586992, 4978400.080824336], [-8234989.529586992, 4978398.886495769], [-8234996.695558393, 4978394.109181502], [-8234997.889886959, 4978394.109181502], [-8234999.0842155265, 4978394.109181502], [-8235001.4728726605, 4978390.526195802], [-8235008.638844062, 4978382.165895834], [-8235009.833172629, 4978382.165895834], [-8235024.165115431, 4978371.416938731], [-8235036.108401099, 4978367.833953031], [-8235038.497058233, 4978367.833953031], [-8235040.885715366, 4978367.833953031], [-8235054.023329602, 4978364.25096733], [-8235055.217658169, 4978355.890667362], [-8235055.217658169, 4978353.502010229], [-8235055.217658169, 4978352.307681662], [-8235055.217658169, 4978351.113353095], [-8235052.829001036, 4978343.947381694], [-8235040.885715366, 4978320.060810357], [-8235033.719743965, 4978318.86648179], [-8235027.748101131, 4978316.477824656], [-8235018.193472596, 4978312.894838956], [-8235016.99914403, 4978311.700510389], [-8234999.0842155265, 4978311.700510389], [-8234964.448687088, 4978323.643796057], [-8234960.865701388, 4978328.421110325], [-8234960.865701388, 4978329.615438892], [-8234956.08838712, 4978337.97573886], [-8234954.894058553, 4978339.170067427], [-8234954.894058553, 4978340.364395994], [-8234953.6997299865, 4978341.55872456], [-8234953.6997299865, 4978342.753053127], [-8234953.6997299865, 4978343.947381694], [-8234950.116744285, 4978348.724695961], [-8234944.145101451, 4978355.890667362], [-8234942.950772884, 4978355.890667362], [-8234942.950772884, 4978357.084995929], [-8234941.756444317, 4978358.279324496], [-8234941.756444317, 4978359.473653063], [-8234940.56211575, 4978360.66798163], [-8234922.647187248, 4978378.582910133]]]})
polyline1

Out[101]:

In [95]:

from IPython.display import display_svg, clear_output, display, HTML
import time

In the cell below, we loop through various offset distances and observe how that effects the amount of generalization. The higher the offset, the greater the generalization. The greater the generalization, the fewer the vertices, and as such, lesser is the complexity of the resulting geometry. Compact geometries take up less space (for storage) and memory as well.

In [100]:

time.sleep(2)
for offset in [0,10,15,20,25,30,35,40,45,50]:
    # apply generalization
    generalized_polyline = polyline1.generalize(max_offset=offset)
    
    # display logic to create an animation
    display(HTML(f'<p>Offset: {offset}</p>'))
    display_svg(generalized_polyline)
    time.sleep(1)
    clear_output()

Buffer In the cell below, we loop through various distance measures for the buffer tool. We start with a distance of 1m (the unit is meters as that is the units of the geometry's spatial reference), where the buffer result is quite similar to the source geometry. As the distance increases, the effect of buffer increases, producing thicker polygons. The tool generalizes the resulting buffer polygon to produce a smooth shape and dissolves overlaps. At a distance of 70m, the buffer looks like a simple circle.

In [109]:

time.sleep(2)
for distance in [1,10,15,20,25,30,35,40,45,50,55,60,65,70]:
    # apply buffer
    buffered_polyline = polyline1.buffer(distance=distance)
    
    # display logic to create an animation
    display(HTML(f'<p>Buffer distance: {distance}</p>'))
    display_svg(buffered_polyline)
    time.sleep(1)
    clear_output()

ArcGIS API for Python Guide