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. thefrom_geo_coordinate_string
orto_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.
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.
arcpy
, then to visualize the geometries as SVG graphics, replace the cells that query the object with geom_object.as_arcpy
from arcgis.gis import GIS
from arcgis.geometry import Polygon, Geometry, Point, Polyline
from arcgis.geocoding import geocode
import arcpy
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):
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]]]}
geom1 = Polygon(geom1_json)
geom2 = Polygon(geom2_json)
geom1
geom2
geom_union = geom1.union(geom2)
geom_union
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.
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]]]}
geom3 = Polygon(geom3_json)
geom3
geom_diff = geom_union.difference(geom3)
geom_diff
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()
.
geom_sdiff = geom_union.symmetric_difference(geom3)
geom_sdiff
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
: requiredarcgis.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).
# scroll to the top of this notebook where these two geometries are visualized
geom_intersection = geom1.intersect(second_geometry=geom2, dimension=4)
geom_intersection
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
.
# True, because the intersection of geom1 and geom2 is partial shape of geom1, and is not
# equivalent to either one
geom1.overlaps(geom2)
# False, because the intersection of geom1 and geom_intersection is geom_intersection,
# and that is equivalent to the second geometry
geom1.overlaps(geom_intersection)
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.
geom1.equals(geom1)
geom_intersection.equals(geom_sdiff)
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 thegeneralize()
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 ofbuffer()
operations against different geomtry types can be found in Fig 3, which also shows that results ofbuffer()
are always polygon objects.
To illustrate generalization, let us construct a complex polyline object
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
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.
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.
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()