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.functionsdirectly without initiating any geometry objects, e.g. thefrom_geo_coordinate_stringorto_geo_coordinate_stringmethods, 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()
a7. Find the nearest point¶
Imagine that you are jogging on a trail(drawn as the light blue line segment on the map) inside Central Park and looking for a water fountain. You know the water fountain is located inside the pavilion (as marked as the blue pin). You want to know: What is the shortest path from the trail to pavilion?
The method query_point_and_distance(second_geometry, use_percentage=False) can find the point on a polyline nearest to the second_geometry and the distance between these points. It also returns information about which side of the line the second_geometry is on as well as the distance along the line where the nearest point occurs.
can be used to find the point on the polyline nearest to the in_point and the distance between those points, and return information about the side of the line the in_point is on as well as the distance along the line where the nearest point occurs.
access_point = Point({'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
'x': -8234818.501757936, 'y': 4978337.398475052})
access_polyline = Polyline({'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
'paths': [[[-8234920.019686119, 4978251.645698531], [-8234885.38415768, 4978242.091069995],
[-8234849.554300674, 4978245.674055696], [-8234808.947129401, 4978255.228684231],
[-8234773.117272396, 4978251.645698531], [-8234751.619358192, 4978239.702412861],
[-8234736.093086823, 4978215.815841525]]]})
map1 = gis.map()
map1.basemap = "satellite"
map1.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
'xmin': -8235024.105413707,
'ymin': 4978145.828180327,
'xmax': -8234435.301430117,
'ymax': 4978384.69389375}
map1
map1.draw(access_polyline,
symbol = {
"type": "simple-line",
"color": "red",
"width": "4px",
"style": "short-dot"
})
map1.draw(access_point,
symbol = {"angle":0,"xoffset":2,"yoffset":8,"type":"esriPMS",
"url":"http://static.arcgis.com/images/Symbols/Basic/RedShinyPin.png",
"contentType":"image/png","width":24,"height":24})
Now let us use
query_point_and_distance(second_geometry, use_percentage=False)to find the point on the polyline (access_polyline) nearest to the second_geometry (access_point). The acces_res variable assigned to the results contains:
- the PointGeometry object for the nearest point,
- the x-, and y- coordinates of this point
- the distance between those points
access_res = access_polyline.query_point_and_distance(access_point)
access_res
access_res[0].type
access_res[0].JSON
With the information we can now draw the nearest point on the map view as well, with the selected symbology (which is the green pin marker).
a_p = Point({"x":-8234808.9471294014,"y":4978255.2286842307,"spatialReference":{"wkid":102100,"latestWkid":3857}})
map1.draw(a_p,
symbol = {"angle":0,"xoffset":2,"yoffset":8,"type":"esriPMS",
"url":"http://static.arcgis.com/images/Symbols/Basic/GreenShinyPin.png",
"contentType":"image/png","width":24,"height":24})
a8. contains()¶
The method contains(second_geometry, relation=None) indicates if the base geometry contains the comparison geometry.
access_polyline.contains(a_p)
The contains method accepts an optional parameter relation, which can be used to specify the spatial relationship type.
BOUNDARY- Relationship has no restrictions for interiors or boundaries.CLEMENTINI- Interiors of geometries must intersect. Specifying CLEMENTINI is equivalent to specifying None. This is the default.PROPER- Boundaries of geometries must not intersect.
access_polyline.contains(a_p, relation = "PROPER")
access_polygon = Polygon({'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
'rings': [[[-8234920.019686119, 4978251.645698531], [-8234885.38415768, 4978242.091069995],
[-8234849.554300674, 4978245.674055696], [-8234808.947129401, 4978255.228684231],
[-8234773.117272396, 4978251.645698531], [-8234751.619358192, 4978239.702412861],
[-8234736.093086823, 4978215.815841525], [-8234663.657066563, 4978349.222320599],
[-8234737.705437716, 4978225.012149632], [-8234935.963979837, 4978263.2306637755]]]})
access_polygon
access_polygon.contains(access_polyline)
access_polygon.contains(access_polyline, relation = "BOUNDARY")
a9. Clip a geometry object¶
The method clip(envelope) constructs the intersection of the geometry and the specified extent. Note that, if ArcPy is not installed, none is returned. The one and required input parameter is,
envelope: required tuple. The tuple must be in the form of a tuple (XMin, YMin, XMax, YMax) where each value represents the lower left bound and upper right bound of the extent.
# (XMin, YMin, XMax, YMax)
x = access_polygon.centroid[0]
y = access_polygon.centroid[1]
envelope_tuple = (x-10, y-10, x+10, y+10)
clipped = access_polygon.clip(envelope_tuple)
clipped.JSON
clipped.as_arcpy
You can find the source geometry (dark shaded polygon) and the clipped geometry (the contained polygon with white diagnoal crosses) in the map below:
map2 = gis.map()
map2.basemap = "satellite"
map2
map2.height= "650px"
location = geocode("Central Park, New York")[0]
map2.extent = location['extent']
map2.draw(access_polygon)
map2.draw(clipped,
symbol = {
"type" : "esriSFS",
"style" : "esriSFSDiagonalCross",
"outline" : "blue",
"color": "lightblue"
})
Up till this point, you have tried creating geometry objects and performed spatial operations such as union and intersect, from the built-in methods of these geometry objects using the OOP pattern. Doing so executed the operations using local geometry engines (ArcPy or Shapely). Next, let us explore the other way in performing spatial operations.
b. Pattern using arcgis.geometry.functions: uses server-side geometry engine¶
Besides calling methods off the geometry object as discussed in section (a), users can also call spatial operations directly from arcgis.geometry.functions. In the following section, you will see spatial operations such as union and intersect being performed with global functions.
from arcgis.geometry.functions import union, intersect, difference, generalize, buffer, relation
b1. Union¶
Compared to what is used in section (a1) in performing union operation (a.k.a. geom_union = geom1.union(geom2)), calling the union() method directly on a geometry resource service would require these input parameters:
spatial_ref- The well-known ID of the spatial reference or a spatial reference json object.geometries- The array of geometries to be unioned.gis- The active GIS instance.
geom1
geom2
geom_union_b = union(spatial_ref={"wkid":102100,"latestWkid":3857},
geometries = [geom1, geom2],
gis = gis)
geom_union_b
b2. intersect() V.S. overlaps()¶
The intersect function arcgis.geometry.functions.intersect(spatial_ref, geometries, geometry, gis=None) is performed on a geometry service resource. This function constructs the set-theoretic intersection between an array of geometries and another geometry. The dimension of each resultant geometry is the minimum dimension of the input geometry in the geometries array and the other geometry specified by the geometry parameter. Inputs include:
spatial_ref- The well-known ID or a spatial reference JSON object for the input geometries.geometries- An array of points, multipoints, polylines, or polygons. The structure of each geometry in the array is the same as the structure of the JSON geometry objects returned by the ArcGIS REST API.geometry- A single geometry of any type with a dimension equal to or greater than the elements of geometries.
Let's see how it compares to what we used in section (a4) - geom_intersection = geom1.intersect(second_geometry=geom2, dimension=4):
geom_intersection_b = intersect(spatial_ref={"wkid":102100,"latestWkid":3857},
geometries=[geom1],
geometry=geom2,
gis=gis)
geom_intersection_b
geom_intersection_b[0]
Unlike what's stated in section (a4) that a simple first_geometry.overlaps(second_geometry) is just sufficient to perform an overalapping operation which then returns 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, there is no direct match of overlaps function in the arcgis.geometry.functions.
However, we can write our own function as in my_overlap in situations when OOP pattern is not applicable, and it includes two operations - first, it creates the interection 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.
def my_overlap(spatial_ref, geometry1, geometry2, gis):
intersection = intersect(spatial_ref,
geometries=[geometry1],
geometry=geometry2,
gis=gis)
if intersection is None:
return False
elif intersection[0].equals(geometry1) or intersection[0].equals(geometry2):
return False
else:
return True
# False, because the intersection of geom1 and geom_union is geom1, and that is equivalent to the first geometry
my_overlap(3857, geom1, geom_union, gis)
# True, because the intersection of geom1 and geom_intersection is partial shape of geom_intersection, and is not
# equivalent to either one
my_overlap(3857, geom1, geom_intersection, gis)
b3. Difference¶
In order to construct the set-theoretic difference between each element of an array of geometries and another geometry (a.k.a. the so-called difference geometry), users need to perform a difference function on a geomtry service resource, which requires these inputs:
geometries- An array of points, multipoints, polylines or polygons. The structure of each geometry in the array is the same as the structure of the JSON geometry objects returned by the ArcGIS REST API.spatial_ref- The well-known ID of the spatial reference or a spatial reference JSON object for the input geometries.geometry- A single geometry of any type and of a dimension equal to or greater than the elements of geometries. The structure of geometry is the same as the structure of the JSON geometry objects returned by the ArcGIS REST API. The use of simple syntax is not supported.
In other words, when performing arcgis.geometry.functions.difference(geometries=[A1, A2, ..., AN], spatial_ref, geometry=B, gis=None), let B be the difference geometry. For each geometry in [A1, A2, ..., AN], in the input geometry array, it constructs a new geometry An-B (that n is numbered from 1 to N).
display_svg(geom_union)
geom_diff_b = difference(geometries=[geom1, geom2, geom_union],
spatial_ref=3857,
geometry=geom3,
gis=gis)
geom_diff_b
len(geom_diff_b)
display_svg(geom_diff_b[1]), display_svg(geom_diff_b[2])
""" Compared to geom_diff derived from `geom_diff = geom_union.difference(geom3)`
"""
geom_diff_b[2].equals(geom_diff)
b4. Symmetric difference¶
Again there is no direct match of symmetric_difference function in the arcgis.geometry.functions. In situations when OOP pattern of this call is not applicable, you can construct your own customized function that creates the geometry which is the union of two geometries minus the instersection of those geometries.
def my_symmetric_difference(spatial_ref, geometry1, geometry2, gis):
union_b = union(spatial_ref,
geometries = [geometry1, geometry2],
gis = gis)
intersect_b = intersect(spatial_ref,
geometries=[geometry1],
geometry=geometry2,
gis=gis)
symm_diff_b = difference(geometries=[union_b],
spatial_ref=spatial_ref,
geometry=intersect_b[0],
gis=gis)
return symm_diff_b[0]
display_svg(geom3)
geom_sdiff_b = my_symmetric_difference(3857, geom_union, geom3, gis)
geom_sdiff_b
geom_sdiff_b.equals(geom_sdiff)
b5. Equals¶
There is no direct match of equals() function in the arcgis.geometry.functions. In situations when OOP pattern of this call is not applicable, you can construct your own customized function that checks if two geometries are equal.
def my_equals(spatial_ref, geometry1, geometry2, gis):
def my_1direction_diff(spatial_ref, geom1, geom2, gis):
diff_b = difference(geometries=[geom1],
spatial_ref=spatial_ref,
geometry=geom2,
gis=gis)
# print(diff_b[0])
if diff_b[0] is None:
return True
elif "rings" in diff_b[0] and len(diff_b[0]["rings"])==0:
return True
else:
return False
if geom1.type == geom2.type:
res1 = my_1direction_diff(spatial_ref, geometry1, geometry2, gis)
res2 = my_1direction_diff(spatial_ref, geometry2, geometry1, gis)
return res1 and res2
else:
return False
# Compared to "geom_intersection.equals(geom_diff)"
my_equals(3857, geom_intersection, geom_diff, gis)
# Compared to "geom_intersection.equals(geom_sdiff)"
my_equals(3857, geom_sdiff, geom_intersection, gis)
b6. generalize() V.S. buffer()¶
In section (a6), we have already discussed the difference between these two spatial operations:
generalize(max_offset)- Creates a new simplified geometry using a specified maximum offset tolerance.buffer(distance)- Constructs a polygon at a specified distance from the geometry.
Now let us explore how these two spatial operations can be done by executing the functions on the server-side.
The function arcgis.geometry.functions.generalize(spatial_ref, geometries, max_deviation, deviation_unit, gis=None) is performed on a geometry service resource. The generalize function simplifies the input geometries using the Douglas-Peucker algorithm with a specified maximum deviation distance. The output geometries will contain a subset of the original input vertices.
One of the input parameters, the deviation_unit represents the unit for maximum deviation. If a unit is not specified, the units are derived from spatial_ref. For a list of valid units, see esriSRUnitType Constants. 9001 here means Meters.
geom_generalized_b = generalize(spatial_ref={"wkid":102100,"latestWkid":3857},
geometries=[polyline1],
max_deviation=50,
deviation_unit="9001",
gis=gis)
geom_generalized_b[0]
geom_generalized2_b = generalize(spatial_ref={"wkid":102100,"latestWkid":3857},
geometries=[polyline1],
max_deviation=100,
deviation_unit="9001",
gis=gis)
geom_generalized2_b[0]
The buffer function arcgis.geometry.functions.buffer(geometries, in_sr, distances, unit, out_sr=None, buffer_sr=None, union_results=None, geodesic=None, gis=None) is performed on a geometry service resource. The result of this function is buffered polygons at the specified distances for the input geometry array. Options are available to union buffers and to use geodesic distance:
union_results- If True, all geometries buffered at a given distance are united into a single (gis,possibly multipart) polygon, and the united geometry is placed in the output array. The default is False.geodesic- Set geodesic to True to buffer the input geometries using geodesic distance. Geodesic distance is the shortest path between two points along the ellipsoid of the earth. If geodesic is set to False, the 2D Euclidean distance is used to buffer the input geometries. The default value depends on the geometry type, unit and buffer_SR.
geom_buffered_b = buffer(geometries = [polyline1],
in_sr=3857,
distances=50,
unit="9001",
out_sr=None,
buffer_sr=None,
union_results=None,
geodesic=None,
gis=gis)
geom_buffered_b[0]
The unit parameter indicates the units for calculating each buffer distance. If not specified, the units are derived from bufferSR. If bufferSR is not specified, the units are derived from in_sr. In the example above, we are calling the buffer method with unit set to be a constant string "9001". Equivalently we can use a constant Enum, e.g. LengthUnits.METER, or the integer 9001 to specify the same requirement.
from arcgis.geometry import LengthUnits
geom_buffered_b = buffer(geometries = [polyline1],
in_sr=3857,
distances=50,
unit=LengthUnits.METER,
out_sr=None,
buffer_sr=None,
union_results=None,
geodesic=None,
gis=gis)b7. Distance¶
From section (a7) we have figured out that a_p is the closest point on the trail from the pavilion, is there a way to confirm that the distance from a_p to the access_point is in fact the shortest?
The distance function arcgis.geometry.functions.distance(spatial_ref, geometry1, geometry2, distance_unit='', geodesic=False, gis=None) is performed on a geometry service resource. It reports the 2D Euclidean or geodesic distance between the two geometries. We can take advantage of the distance() method in validation.
from arcgis.geometry import distance
radius = distance( spatial_ref=3857,
geometry1=access_point,
geometry2=a_p,
distance_unit='9001',
geodesic=False,
gis=gis)
radius['distance']
access_point_buffer = buffer( geometries = [access_point],
in_sr=3857,
distances=radius['distance'],
unit="9001",
out_sr=3857,
buffer_sr=3857,
union_results=None,
geodesic=None,
gis=gis)
access_point_buffer[0].spatialReference = {"wkid":102100,"latestWkid":3857}
map3 = gis.map()
map3.basemap = "satellite"
map3
map3.height= "650px"
location = geocode("Central Park, New York")[0]
map3.extent = location['extent']
map3.draw(access_polyline,
symbol = {
"type": "simple-line",
"color": "lightblue",
"width": "4px",
"style": "short-dot"
})
map3.draw(access_point_buffer[0],
symbol = {
"type" : "esriSFS",
"style" : "esriSFSDiagonalCross",
"outline" : "blue",
"color": "lightblue"
})
map3.draw(access_point,
symbol = {"angle":0,"xoffset":2,"yoffset":8,"type":"esriPMS",
"url":"http://static.arcgis.com/images/Symbols/Basic/BlueShinyPin.png",
"contentType":"image/png","width":24,"height":24})
map3.draw(a_p,
symbol = {"angle":0,"xoffset":2,"yoffset":8,"type":"esriPMS",
"url":"http://static.arcgis.com/images/Symbols/Basic/GreenShinyPin.png",
"contentType":"image/png","width":24,"height":24})
It is obvious that the buffer created only intersects with the access_polyline at one point - a_p, hence we can validate that the shortest path from trail to pavilion is to leave at a_p, which well corresponds to the results seen in section (a7).
b8. Determine the spatial relation between two geometries¶
The relation function arcgis.geometry.functions.relation(geometries1, geometries2, spatial_ref, spatial_relation='esriGeometryRelationIntersection', relation_param='', gis=None) is performed on a geometry service resource. This function determines the pairs of geometries from the input geometry arrays that participate in the specified spatial relation. Both arrays are assumed to be in the spatial reference specified by spatial_ref, which is a required parameter. Geometry types cannot be mixed within an array. The relations are evaluated in 2D. In other words, z coordinates are not used.
The spatial relationship to be tested between the two input geometry arrays is defined by spatial_relation, and can be chosen from these values:
esriGeometryRelationCrossesriGeometryRelationDisjointesriGeometryRelationInesriGeometryRelationInteriorIntersectionesriGeometryRelationIntersectionesriGeometryRelationLineCoincidenceesriGeometryRelationLineTouchesriGeometryRelationOverlapesriGeometryRelationPointTouchesriGeometryRelationTouchesriGeometryRelationWithinesriGeometryRelationRelation
Take the Point Geometry a_p and Polyline Geometry access_polyline created previously. Let us explore if the specified spatial relation exists between the input geometry arrays. Here, a_p is on the access_polyline.
"""esriSpatialRelationIn is the same as IRelationalOperator::Within. esriSpatialRelationWithin also allows polylines
that are strictly on the boundaries of polygons to be considered “in” the polygon. This case is disallowed by the former
relation (and disallowed by IRelationalOperator::Within).
"""
relation([a_p], [access_polyline],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationWithin",
gis = gis)
relation([a_p], [access_polyline],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationIn",
gis = gis)
"""`esriGeometryRelationInteriorIntersection` same as esriSpatialRelationIntersection, but excludes intersections that
occur only at boundaries (two touching polygons for example). The interior of a point is considered to be the point
itself and the interior of a polyline excludes the endpoints of all its parts.
"""
relation([a_p], [access_polyline],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationIntersection",
gis = gis)
relation([a_p], [access_polyline],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationInteriorIntersection",
gis = gis)
"""esriSpatialRelationPointTouch and esriSpatialRelationLineTouch are boundary intersectionsclassified by dimension of
intersection. esriSpatialRelationTouch is the union of those two and is equivalent to the ArcObjects IRelationalOperator
Touch method.
"""
relation([a_p], [access_polyline],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationPointTouch",
gis = gis)
The difference of esriGeometryRelationWithin and esriGeometryRelationIn, and that of esriGeometryRelationIntersection and esriGeometryRelationInteriorIntersection are not as significant in the previous example. Next, let us look at a different use case of access_polyline and access_polygon that the former is overlappig with partial edges of the latter.
"""esriSpatialRelationIn is the same as IRelationalOperator::Within. esriSpatialRelationWithin also allows polylines
that are strictly on the boundaries of polygons to be considered “in” the polygon. This case is disallowed by the former
relation (and disallowed by IRelationalOperator::Within).
"""
relation([access_polyline], [access_polygon],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationWithin",
gis = gis)
relation([access_polyline], [access_polygon],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationIn",
gis = gis)
"""`esriGeometryRelationInteriorIntersection` same as esriSpatialRelationIntersection, but excludes intersections that
occur only at boundaries (two touching polygons for example). The interior of a point is considered to be the point
itself and the interior of a polyline excludes the endpoints of all its parts.
"""
relation([access_polyline], [access_polygon],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationIntersection",
gis = gis)
relation([access_polyline], [access_polygon],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationInteriorIntersection",
gis = gis)
"""esriSpatialRelationPointTouch and esriSpatialRelationLineTouch are boundary intersectionsclassified by dimension of
intersection. esriSpatialRelationTouch is the union of those two and is equivalent to the ArcObjects IRelationalOperator
Touch method.
"""
relation([access_polyline], [access_polygon],
spatial_ref={'spatialReference': {'latestWkid': 3857, 'wkid': 102100}},
spatial_relation="esriGeometryRelationLineTouch",
gis = gis)
As shown above, though access_polyline is considered to be esriGeometryRelationWithin inside access_polygon, but it is not esriGeometryRelationIn inside access_polygon. Similarly, access_polyline is esriGeometryRelationIntersection intersecting with access_polygon, but not esriGeometryRelationInteriorIntersection.
b9. Clip a geometry object¶
There is no direct match in the arcgis.geometry.functions that is equivalent to clips(), but we can take advantage of intersect() and write our own customized function that creates the envelope object to clip with, and performs intersect() with source geometry to generate the equivalent clipped result.
def my_clip(spatial_ref, geometry1, envelope_tuple, gis):
x_min, y_min, x_max, y_max = envelope_tuple
envelope_json = {"rings":[[[x_min,y_min],[x_min,y_max],
[x_max,y_max],[x_max,y_min],
[x_min,y_min]]],
"spatialReference":spatial_ref}
geom_intersection_b = intersect(spatial_ref=spatial_ref,
geometries=[geometry1],
geometry=Polygon(envelope_json),
gis=gis)
return geom_intersection_b[0]
clipped_b = my_clip(spatial_ref={"wkid":102100,"latestWkid":3857},
geometry1=access_polygon,
envelope_tuple=envelope_tuple,
gis=gis
)
clipped_b.as_arcpy
You can find the source geometry (dark shaded polygon) and the clipped geometry (the contained polygon with white diagnoal crosses) in the map below:
map4 = gis.map()
map4.basemap = "satellite"
map4
map4.height= "650px"
location = geocode("Central Park, New York")[0]
map4.extent = location['extent']
map4.draw(access_polygon)
map4.draw(clipped_b,
symbol = {
"type" : "esriSFS",
"style" : "esriSFSDiagonalCross",
"outline" : "blue",
"color": "lightblue"
})
b10. Areas and Lengths¶
The areas_and_lengths function calculates areas and perimeter lengths for each Polygon object specified in the input array. The input parameters include:
polygons- The array of polygons whose areas and lengths are to be computed.length_unit- The length unit in which the perimeters of polygons will be calculated.- If
calculation_typeisplanar, then length_unit can be anyesriUnitsconstant. - If
calculation_typeis notplanar, thenlength_unitmust be a linearesriUnitsconstant, such asesriSRUnit_Meter(i.e.9001|LengthUnits.METER) oresriSRUnit_SurveyMile(i.e.9035|LengthUnits.SURVEYMILE). - If
length_unitis not specified, the units are derived fromspatial_refparameter. Ifspatial_refis not specified, the units are in meters. - For a list of valid units, see `esriSRUnitType Constants` and `esriSRUnit2Type Constants`.
- If
area_unit- The area unit in which areas of polygons will be calculated.- If
calculation_typeisplanar, thenarea_unitcan be anyesriAreaUnitsconstant. - If
calculation_typeis notplanar, thenarea_unitmust be a linear esriUnits constant such asAreaUnits.SQUAREMETERS(i.e.{"areaUnit": "esriSquareMeters"}) orAreaUnits.SQUAREMILES(i.e.{"areaUnit": "esriSquareMiles"}). - If
area_unitis not specified, the units are derived fromspatial_ref. Ifspatial_refis not specified, then the units are meters. - For a list of valid units, see `esriAreaUnits Constants`. The list of valid
esriAreaUnits constantsinclude, esriSquareInches | esriSquareFeet | esriSquareYards | esriAcres | esriSquareMiles | esriSquareMillimeters | esriSquareCentimeters | esriSquareDecimeters | esriSquareMeters | esriAres | esriHectares | esriSquareKilometers.
- If
calculation_type- The type defined for the area and length calculation of the input geometries. The type can be one of the following values:planar- Planar measurements use 2D Euclidean distance to calculate area and length. This should only be used if the area or length needs to be calculated in the givenSpatialReference. Otherwise, usepreserveShape.geodesic- Use this type if you want to calculate an area or length using only the vertices of thePolygonand define the lines between the points as geodesic segments independent of the actual shape of thePolygon. A geodesic segment is the shortest path between two points on an ellipsoid.preserveShape- This type calculates the area or length of the geometry on the surface of the Earth ellipsoid. The shape of the geometry in its coordinate system is preserved.
There are other optional parameters to be specify in the input arguments. For more information, check out the API ref. Now let's look at two use cases below:
The traditional use case¶
For API version <= 1.9.0, length_unit shall be input as constant integer value (as found in esriSRUnitType Constants or `esriSRUnit2Type Constants`) while area_unit only accepts a dict object where the a list of valid units can be found in `esriAreaUnits Constants`.
from arcgis.geometry.functions import areas_and_lengths
areas_and_lengths(polygons =[access_polygon],
length_unit = 9001,
spatial_ref=3857,
area_unit = {"areaUnit": "esriSquareMeters"},
calculation_type = "planar")
The enhanced use case¶
After 1.9.0, there is an easier way to represent length_unit and area_unit. Users can use enum declared in the arcgis.geometry.functions.AreaUnits and arcgis.geometry.functions.LengthUnits as input, e.g.
- the
length_unitparameter can be represented byLengthUnits.BRITISH1936FOOT | LengthUnits.GOLDCOASTFOOT | LengthUnits.INTERNATIONALCHAIN | LengthUnits.INTERNATIONALLINK | LengthUnits.INTERNATIONALYARD | LengthUnits.METER | LengthUnits.FOOT | LengthUnits.SURVEYFOOT | LengthUnits.CLARKEFOOT | LengthUnits.FATHOM | etc. - while the
area_unitparameter can be represented byAreaUnits.SQUAREINCHES | AreaUnits.SQUAREFEET | AreaUnits.SQUAREYARDS | AreaUnits.ACRES | AreaUnits.SQUAREMILES | AreaUnits.SQUAREMILLIMETERS | AreaUnits.SQUARECENTIMETERS | AreaUnits.SQUAREDECIMETERS | AreaUnits.SQUAREMETERS | AreaUnits.ARES | AreaUnits.HECTARES | AreaUnits.SQUAREKILOMETERS | etc.
from arcgis.geometry.functions import LengthUnits, AreaUnits
res = areas_and_lengths(polygons =[access_polygon],
length_unit = LengthUnits.METER,
area_unit = AreaUnits.SQUAREMETERS,
calculation_type = "planar",
future = True)
res
Also note that, when future is set to True, the areas_and_lengths function will return a GeometryJob object which indicates the job is processed asynchronously. When the job is done, querying its result will bring out a dict object of areas, and lengths for the list of polygons.
res.result()
Conclusion¶
In this notebook, we explored the two approaches to performing spatial operations - (a) Object Oriented Programming (OOP) pattern - (calling operations directly off Geometry object). This pattern utilizes local geometry engines (ArcPy or Shapely) and performs the computation on the machine running the Python kernel (b) calling functions from arcgis.geometry.functions which sends the spatial data over to a Geometry Service running on the active GIS connection, performs the operations remotely and returns the results as Geometry objects.
In most cases, approach (a) is simplier and more staight forward, and is recommended for users getting to use the ArcGIS API for Python. Also, another advantage of (a) is that it does not consume credits and is performant for larger datasets. However, (a) uses local geom engines such as shapely or arcpy, in situations where neither shapely nor arcpy is installed on the execution environment, approach (b) is the alternative.