Part 2 - Working with Geometries

Creating geometries

The arcgis.geometry module defines geometry types for working with geographic features in a GIS. It provides functions which use geometric types as input and output as well as functions for easily converting geometries between different representations.

Several functions in this module accept geometries represented as Python dictionary objects. To get started, we import the following classes and functions from the geometry module.

Input
from arcgis.gis import GIS
from arcgis.geocoding import geocode
from arcgis.geometry import lengths, areas_and_lengths, project
from arcgis.geometry import Point, Polyline, Polygon, Geometry
import pandas as pd
Input
gis = GIS('home')

Object Model Diagram (OMD) of the geometry module

The picture below illustrates how the geometry module is organized (which showcase only the important objects, sub-modules, and properties via the OMD):

The arcgis.geometry module contains classes and utility functions to construct geometry objects such as Points, Polylines, Polygons and perform spatial operations on them. The arcgis.geometry.Geometry class is a base class from which specific geometry classes inherit. Users generally do not instantiate the base class, but work directly with one or more child classes such as Point, Polyline, Polygon etc. The base class provides geometry operations such as clip(), difference(), buffer()... which the child classes inherit. The arcgis.geometry.functions sub-module contains an alternate set of spatial operations. There difference between these two sets of operations is explained in the section 'Two patterns of applying spatial operations'.

The SpatialReference class is used to represent spatial reference and coordinate system information. The arcgis.geometry.filters module provides functions that can be used to define spatial relationships to be used in queries against feature layers and imagery layers.

The layout diagram further shows how Feature and spatially enabled DataFrame objects from arcgis.features module return components of the geometry module.

Creating Point objects

A point contains x and y fields along with a spatialReference field. A point can also contain m and z fields as well, representing the vertical and linear referencing system coordinates. A point is empty when its x field is present and has the value null or the string NaN. An empty point has no location in space.

Input
pt = Point({"x" : -118.15, "y" : 33.80, 
            "spatialReference" : {"wkid" : 4326}})
pt
Output

As shown above, you can create a Point geometry using a dictionary. The x and y key value pairs in this example contain longitude and latitude respectively. The spatialReference dictionary with the wkid kvp specifies the coordinate system in which x and y are in. Thus, you could have passed it X and Y values from a projected coordinate system and constructed the same point by specifying the appropriate wkid.

When using the Jupyter Notebook (or ArcGIS Notebook) interface, you can query a geometry and get a visual representation as shown in the cell earlier. Alternately you can check the validity of a geometry by querying the is_valid() method.

Input
pt.is_valid()
Output
True
Input
print(pt.is_empty)
False

To get the geometry type, use the type property:

Input
print(pt.type)
Point

To get the object type in Python, use the built-in type function:

Input
type(pt)
Output
arcgis.geometry._types.Point
Input
map0 = gis.map("Port of Long Beach")
map0.basemap = "satellite"
map0.zoom = 6
map0
Input
pt_sym = {
    "type": "esriSMS",
    "style": "esriSMSDiamond",
    "color": [255,140,0,255],        
    "size": 14,
    "angle": 0,
    "xoffset": 0,
    "yoffset": 0,
    "outline": {
        "color": [255,140,0,255],
        "width": 1}
}
map0.draw(pt, symbol=pt_sym)

Creating Polyline objects

A polyline contains an array of paths or curvePaths and a spatialReference. For polylines with curvePaths, see the sections on JSON curve object and Polyline with curve. Each path is represented as an array of points, and each point in the path is represented as an array of numbers. A polyline can also have boolean-valued hasM and hasZ fields. An empty polyline is represented with an empty array for the paths field. Nulls and/or NaNs embedded in an otherwise defined coordinate stream for polylines/polygons is a syntax error.

First, let us try initializing an invalid polyline object, and see how it is handled.

Input
line = {
  "paths" : [[[-97.06138],[-97.06133,32.836],[-97.06124,32.834],[-97.06127,32.832]],
             [[-97.06326,32.759],[-97.06298,32.755]]],
  "spatialReference" : {"wkid" : 4326}
}
polyline = Polyline(line)
Input
print(polyline.spatialReference)
{'wkid': 4326}
Input
print(polyline.is_valid())
False

The geometry we just created is invalid, yet the API allowed you to still create it. This because, the API does not check for validity in an eager fashion. Rather it checks when is_valid() is called, in an on-demand fashion. This is for performance reasons.

Input
map1 = gis.map()
map1.basemap = "dark-gray"
map1.center = {  'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
                 'x': -10801694.467855213,
                 'y': 3868771.7699954524}
map1.zoom = 11.0
map1
Input
map1.draw(polyline) # would not draw anything

That did not draw anything. Now, try a valid polyline geometry. In the example below, we specify coordinates in latitude and longitude and specify the coordinate system appropriately as GCS WGS 84 using the wkid 4326.

Input
line1 = {
  "paths" : [[[-97.06138,32.837],[-97.06133,32.836],[-97.06124,32.834],[-97.06127,32.832]],
             [[-97.06326,32.759],[-97.06298,32.755]]],
  "spatialReference" : {"wkid" : 4326}
}
polyline1 = Polyline(line1)
print(polyline1.is_valid())
True
Input
polyline1
Output
Input
map2 = gis.map()
map2.basemap = "dark-gray"
map2.center = {  'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
                 'x': -10801694.467855213,
                 'y': 3868771.7699954524}
map2.zoom = 11.0
map2
Input
sym_poly_aoi = {
  "type": "esriSFS",
  "style": "esriSFSSolid",
  "color": [0,0,0,0],
    "outline": {
     "type": "esriSLS",
     "style": "esriSLSSolid",
     "color": [0,255,0,255],
     "width": 3}
}
map2.draw(polyline1, symbol = sym_poly_aoi)

Creating Polygon objects

A polygon contains an array of rings or curveRings and a spatialReference. Each ring is represented as an array of points. The first point of each ring is always the same as the last point. Each point in the ring is represented as an array of numbers. A polygon can also have boolean-valued hasM and hasZ fields. An empty polygon is represented with an empty array for the rings field. Nulls and/or NaNs embedded in an otherwise defined coordinate stream for polylines/polygons is a syntax error.

Exterior rings are oriented clockwise, while holes are oriented counter-clockwise. Rings can touch at a vertex or self-touch at a vertex, but there should be no other intersections.

Input
polygon1 = Polygon({'spatialReference': {'latestWkid': 4326}, 
                'rings': [[[-97.06587202923951, 32.75656343500563], [-97.07033522518535, 32.75454232619796],
                           [-97.07179434702324, 32.75443405154119], [-97.073596791488, 32.75475887587208],
                           [-97.07501299810983, 32.75475887587208], [-97.07492716677937, 32.75616643554153],
                           [-97.07595713555828, 32.75602207118053], [-97.07115061698558, 32.75887321736912],
                           [-97.06930525730476, 32.75890930713694], [-97.06479914614289, 32.75739351976198],
                           [-97.06587202923951, 32.75656343500563]]]
                })
Input
polygon1
Output
Input
print(polygon1.type)
Polygon
Input
type(polygon1)
Output
arcgis.geometry._types.Polygon
Input
map3 = gis.map()
map3.basemap = "topo"
map3.zoom = 16
map3.center = {'x': -97.05815464365813, 'y': 32.75494892021667, 
               "spatialReference" : {"wkid" : 4326}}
map3
Input
map3.draw(polygon1)

Creating geometries with a different spatial reference

In the examples above, all geometries were created using latitude and longitude. Hence we used a spatial reference of 4326 which corresponds to the WGS84 GCS. In the cell below, we will create the same polygon object, but with a PCS.

Input
polygon1_proj = Polygon({'spatialReference': {'wkid': 3857}, 
                  'rings': [[[-10806331.461044524, 3862983.679353406], [-10806446.116586955, 3862964.5700963344], 
                             [-10805911.057388945, 3863341.9779235027], [-10805705.632875424, 3863346.7552377703], 
                             [-10805204.01487729, 3863146.1080385167], [-10805323.447733987, 3863036.229810354],  
                             [-10805820.288417853, 3862768.7002113485], [-10805982.717102963, 3862754.368268545], 
                             [-10806183.364302218, 3862797.3640969563], [-10806341.01567306, 3862797.3640969563], 
                             [-10806331.461044524, 3862983.679353406]]]
                })

The PCS used here is Web Mercator Auxiliary Sphere with a wkid of 3857. The coordinates is in meters.

Input
map3_proj = gis.map()
map3_proj.basemap = "topo"
map3_proj.zoom = 16
map3_proj.center = {'x': -97.05815464365813, 'y': 32.75494892021667, 
               "spatialReference" : {"wkid" : 4326}}
map3_proj
Input
map3_proj.draw(polygon1_proj)

As you can see from both the maps, we were able to create the same geometries using different coordinate systems and they would line up when plotted one over the other.

Text representations of geometries

At any point in time, you can view the text representations of your geometry objects in three forms as shown below:

Input
print("Geometry represented in JSON:\n", polygon1.JSON, 
      "\n\nGeometry represented in WKB:\n", polygon1.WKB, 
      "\n\nGeometry represented in WKT:\n", polygon1.WKT)
Geometry represented in JSON:
 {'type': 'MultiPolygon', 'coordinates': [(((-97.06587202923951, 32.75656343500563), (-97.07033522518535, 32.75454232619796), (-97.07179434702324, 32.75443405154119), (-97.073596791488, 32.75475887587208), (-97.07501299810983, 32.75475887587208), (-97.07492716677937, 32.75616643554153), (-97.07595713555828, 32.75602207118053), (-97.07115061698558, 32.75887321736912), (-97.06930525730476, 32.75890930713694), (-97.06479914614289, 32.75739351976198), (-97.06587202923951, 32.75656343500563)),)]} 

Geometry represented in WKB:
 b'\x01\x06\x00\x00\x00\x01\x00\x00\x00\x01\x03\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x83\xd3P?7DX\xc0lY\x15\x12\xd7`@@b\xfbP_\x80DX\xc0\xe7;\xcb\xd7\x94`@@) QG\x98DX\xc0\xd9-\x85K\x91`@@\x04"Q\xcf\xb5DX\xc0?\x1eX\xf0\x9b`@@\x06jQ\x03\xcdDX\xc0?\x1eX\xf0\x9b`@@\x8f\xb9P\x9b\xcbDX\xc0\xec}\xcf\x0f\xca`@@\xb9GQ{\xdcDX\xc0%c\xcbT\xc5`@@\xc5BQ\xbb\x8dDX\xc0\x914\xf1\xc1"a@@&dQ\x7foDX\xc0?d\xaf\xf0#a@@\x84_Q\xab%DX\xc0?\xcaVE\xf2`@@\x83\xd3P?7DX\xc0lY\x15\x12\xd7`@@' 

Geometry represented in WKT:
 MULTIPOLYGON (((-97.06587202923951 32.75656343500563, -97.07033522518535 32.75454232619796, -97.07179434702324 32.75443405154119, -97.073596791488 32.75475887587208, -97.07501299810983 32.75475887587208, -97.07492716677937 32.75616643554153, -97.07595713555828 32.75602207118053, -97.07115061698558 32.75887321736912, -97.06930525730476 32.75890930713694, -97.06479914614289 32.75739351976198, -97.06587202923951 32.75656343500563)))

You can convert from one spatial reference to another through a process called 'projection'. The arcgis.geometry.project() function can be used for this. This function sends to geometries (in JSON format seen above) to the Geometry Service configured with your GIS. The service performs the projection and sends it back to your Python client.

If you have large geometries, projecting them locally (either using ArcPy or Shapely backend) is recommended over sending them over to the geometry service. For more details on this, checkout part 3 - spatial operations on this guide.

Input
geom1_reprojected = project(geometries = [polygon1_proj], in_sr = 3857, out_sr = 4326)[0]
geom1_reprojected.type
Output
'Polygon'

You can inspect the coordinates by printing them as shown below. Compare those with the geometries printed earlier.

Input
print("Geometry represented in JSON:\n", geom1_reprojected.JSON, 
      "\n\nGeometry represented in WKB:\n", geom1_reprojected.WKB, 
      "\n\nGeometry represented in WKT:\n", geom1_reprojected.WKT)
Geometry represented in JSON:
 {'type': 'MultiPolygon', 'coordinates': [(((-97.06587202923951, 32.75656343500563), (-97.07033522518535, 32.75454232619796), (-97.07179434702324, 32.75443405154119), (-97.073596791488, 32.75475887587208), (-97.07501299810983, 32.75475887587208), (-97.07492716677937, 32.75616643554153), (-97.07595713555828, 32.75602207118053), (-97.07115061698558, 32.75887321736912), (-97.06930525730476, 32.75890930713694), (-97.06479914614289, 32.75739351976198), (-97.06587202923951, 32.75656343500563)),)]} 

Geometry represented in WKB:
 b'\x01\x06\x00\x00\x00\x01\x00\x00\x00\x01\x03\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x83\xd3P?7DX\xc0lY\x15\x12\xd7`@@b\xfbP_\x80DX\xc0\xe7;\xcb\xd7\x94`@@) QG\x98DX\xc0\xd9-\x85K\x91`@@\x04"Q\xcf\xb5DX\xc0?\x1eX\xf0\x9b`@@\x06jQ\x03\xcdDX\xc0?\x1eX\xf0\x9b`@@\x8f\xb9P\x9b\xcbDX\xc0\xec}\xcf\x0f\xca`@@\xb9GQ{\xdcDX\xc0%c\xcbT\xc5`@@\xc5BQ\xbb\x8dDX\xc0\x914\xf1\xc1"a@@&dQ\x7foDX\xc0?d\xaf\xf0#a@@\x84_Q\xab%DX\xc0?\xcaVE\xf2`@@\x83\xd3P?7DX\xc0lY\x15\x12\xd7`@@' 

Geometry represented in WKT:
 MULTIPOLYGON (((-97.06587202923951 32.75656343500563, -97.07033522518535 32.75454232619796, -97.07179434702324 32.75443405154119, -97.073596791488 32.75475887587208, -97.07501299810983 32.75475887587208, -97.07492716677937 32.75616643554153, -97.07595713555828 32.75602207118053, -97.07115061698558 32.75887321736912, -97.06930525730476 32.75890930713694, -97.06479914614289 32.75739351976198, -97.06587202923951 32.75656343500563)))

Working with geometry object

Properties of geometry objects include,

  • centroid: returns the center of the geometry
  • area: returns the area of a polygon feature, or None for all other feature types. The area is in the units of the spatial reference.
  • length: returns the length of the linear feature, or zero for point and multipoint feature types. The length units is the same as the spatial reference.
  • extent: returns the extent of the geometry as a tuple containing (xmin, ymin, xmax, ymax).
  • as_arcpy: returns the class as an arcpy SpatialReference object
  • as_shapely: returns a shapely geometry object

You will see that for geometry objects in different spatial references, the properties return varied results (transformed through the spatial reference).

Input
polygon1_proj.centroid, geom1_reprojected.centroid
Output
((-10805839.610595193, 3863047.3497144855),
 (-97.0705088053709, 32.756647420138925))
Input
polygon1_proj.area, geom1_reprojected.area
Output
(431509.386639783, 2.92840477522626e-05)

The variable geom1_reprojected has the same polygon in GCS spatial reference. Accessing area or length while in GCS may not yield meaningful results as the units are decimal degrees. Hence, the recommendation is to project the geometry to a PCS (such as polygon1_proj in this example which uses Web Mercator Aux Sphere projection with units in meters) and access those properties in projected units.

Input
polygon1_proj.length, geom1_reprojected.length
Output
(2955.798323907705, 0.025605496007056323)
Input
print(polygon1_proj.extent)
print(geom1_reprojected.extent)
(-10806446.116586955, 3862754.368268545, -10805204.01487729, 3863346.7552377703)
(-97.07595713555828, 32.75443405154119, -97.06479914614289, 32.75890930713694)

Geometry Engines

The Python API uses either shapely or arcpy as back-ends (engines) for processing geometries. The API is identical no matter which you use. However, at any point in time, only one engine will be used. The section below shows you how to access your geometries as their corresponding back-end objects.

By default, the geometry module looks for arcpy. If not present in your current conda environment, it looks for shapely. The cell below shows how to access the current object as a native shapely geometry object.

Input
try:
    polygon1_proj_shapely_obj = polygon1_proj.as_shapely
    print(type(polygon1_proj_shapely_obj))
except OSError as e:
    print(e)
<class 'shapely.geometry.multipolygon.MultiPolygon'>

If your environment did not have shapely, you would be presented with an error similar to shown below:

Input
try:
    polygon1_proj.as_shapely
except OSError as e:
    print(e)
[WinError 126] The specified module could not be found

If your condaenvironment has arcpy, you can access the geometry as an arcpy geometry as shown below:

Input
geom1_reprojected.as_arcpy
Output

Besides using properties as_shapely and as_arcpy to detect the geometry engine installed, the built-in function _check_geometry_engine returns a tuple of boolean values representing if arcpy and/or shapely geometry engine is installed.

Input
pt._check_geometry_engine()
(False, True)
Input
pt._HASARCPY
False
Input
pt._HASSHAPELY
True

Creating geometries interactively using the map widget

The examples shown above create geometry objects by specifying the full coordinates. Another way to create simple and small geometries is interactively using the map widget. The map widget supports drawing geometries. In addition, it supports registering callbacks with which you can extract the geometry you just drew.

The following example uses the arcgis.geometry module to compute the length of a path that the user draws on the map. The particular scenario is of a jogger who runs in the Central Park in New York (without gizmos like GPS watches to distract him/her), and wants a rough estimate of their daily runs based on the path taken. The notebook starts out with a satellite map of Central Park in New York:

Input
map5 = gis.map()
map5.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'xmin': -8235706.664189668,
 'ymin': 4977993.551288029,
 'xmax': -8233351.448255569,
 'ymax': 4978949.014141619}
map5.basemap = 'satellite'
map5

Input
# Define the callback function that computes the length.
drawn_line = None
def calc_dist(map1, g):
    print("Computing length of drawn polyline...")
    length = lengths(g['spatialReference'], [g], "", "geodesic")['lengths']
    print("Length: " + str(length[0]) + " m.")
    # update the geometry
    global drawn_line
    drawn_line = g

# Set calc_dist as the callback function to be invoked when a polyline is drawn on the map
drawn_line = map5.on_draw_end(calc_dist)

We want the reader to draw a freehand polyline to indicate the paths that the runner takes for his/her runs. When the drawing operation ends, we use the GIS's Geometry service to compute the length of the drawn path. We can do this by adding an event listener (i.e. a callback function) to the map widget that gets called when drawing is completed (i.e. on_draw_end). The event listener then computes the geodesic length of the drawn geometry using the geometry service and prints it out:

Input
map5.draw("polyline")

Now draw a freehand sketch on the map above, representing the jogger's path. Finish the sketch by double-clicking the mouse. You should notice the length printed below the map.

You can inspect the geometry that you just drew as it was stored in a variable.

Input
drawn_line
Output
{'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'paths': [[[-8235387.89786598, 4978256.38321168],
   [-8235366.399951775, 4978437.921153862],
   [-8235239.801123674, 4978550.188039159],
   [-8235175.307381056, 4978660.066267322],
   [-8235034.3766101515, 4978636.1796959825],
   [-8234848.061353702, 4978490.47161081]]]}

Computing area of a polygon drawn on the map

Similar to the line shown above, you could draw a polygon and compute properties such as area.

Input
map6 = gis.map()
map6.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'xmin': -8235706.664189668,
 'ymin': 4977993.551288029,
 'xmax': -8233351.448255569,
 'ymax': 4978949.014141619}
map6.basemap = 'satellite'
map6
Output

Now draw a polygon on the map representing the route taken by the jogger

Input
# Define the callback function that computes the area.
drawn_polygon = None
def calc_area(map1, g):
    print("Computing area of drawn polygon...")
    area0 = Polygon(g).area
    print("Area (sr=3857): " + str(area0) + " sq m.")
    area = project(geometries=[Polygon(g)], in_sr=g['spatialReference'], out_sr=4326)[0].area
    print("Area (sr=4326): " + str(area) + " degrees")
    global drawn_polygon
    drawn_polygon = g

# Set calc_dist as the callback function to be invoked when a polygon is drawn on the map
map6.on_draw_end(calc_area)
Input
map6.draw("polygon")

Now draw a polygon on map6, finish it by double clicking the mouse pointer. The area in two different coordinate systems (PCS and GCS) are printed. As discussed earlier, to get meaningful results for lengths and areas, it is best to project the geometry to a PCS. The cell below shows the geometry you just drew on the map.

Input
drawn_polygon
Output
{'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'rings': [[[-8234513.649354945, 4978225.887995541],
   [-8234432.43501239, 4978194.8354528],
   [-8234434.823669524, 4978173.337538593],
   [-8234406.159783917, 4978156.616938656],
   [-8234456.32158373, 4978075.402596101],
   [-8234489.762783606, 4978096.9005103065],
   [-8234508.872040678, 4978087.34588177],
   [-8234585.309068965, 4978116.009767379],
   [-8234513.649354945, 4978225.887995541]]]}

Extracting geometries from existing feature layers

Besides creating geometries from scratch, we can also access geometries from existing feature classes, and feature layers - either remotely on the web, or stored locally. Let's first explore how to get geometries from a web feature layer.

Get geometries from a spatially enabled DataFrame object

Next, we will be looking at a path where geometry objects are extracted from Spatially Enabled DataFrame (SeDF) which is obtained from FeatureLayer (in short, the path can be represented as FeatureLayer -> SeDF -> geometry objects).

Input
dc_fl = FeatureLayer('https://maps2.dcgis.dc.gov/dcgis/rest/services/DCGIS_DATA/Administrative_Other_Boundaries_WebMercator/MapServer/10')
dc_df = pd.DataFrame.spatial.from_layer(dc_fl)
display(dc_df.head())
OBJECTID CITY_NAME STATE_CITY CAPITAL WEB_URL AREAKM AREAMILES Shape_Length Shape_Area SHAPE
0 1 Washington 1150000 Y http://www.dc.gov 177.47 68.52 67608.276922 1.774562e+08 {"rings": [[[-8584936.334474642, 4712272.26069...
Input
type(dc_df.spatial)
Output
arcgis.features.geo._accessor.GeoAccessor
Input
display(dc_df.spatial)
Output

If you are currently on arcpy or shapely geometry engines, you should see a preview representing the shape of the geometry or the SeDF object above. If your machine is without geometry engines, then display(dc_df.spatial) would print the JSON representation of the geometry object, e.g.

{'rings': [[[-8584936.334474642, 4712272.260693985],
[-8584936.299852336, 4712272.294257972],
[-8576161.160866661, 4721094.044860408],
[-8561487.41107954, 4706346.502976524],
[-8575944.799488472, 4691870.205657273],
[-8575944.771677606, 4691904.203184697],
[-8575947.014354223, 4692094.371141138],
[-8575956.007432679, 4692617.3203350045],
[-8575958.298881542, 4692731.805968994],
[-8575962.881262433, 4693148.057541993],
[-8575941.7592397, 4693295.905004359],
...]],
'spatialReference': {'wkid': 102100, 'latestWkid': 3857}}

Besides looking at the SeDF object, we can also explore each entry it contains, and its boundingbox property bbox.

Input
poly1 = dc_df.iloc[0].SHAPE
type(poly1)
Output
arcgis.geometry._types.Polygon
Input
display(poly1)
Output
Input
print(type(dc_df.spatial.bbox))
display(dc_df.spatial.bbox)
<class 'arcgis.geometry._types.Polygon'>
Output

To show the geographic relationship between the geometric shapes of the SeDF object and its bbox property, we can draw both geometries on the map:

Input
dc_map = gis.map('Washington DC')
dc_map
Output
Input
dc_map.draw(dc_df.iloc[0].SHAPE, 
            symbol={
                      "type": "esriSFS",
                      "style": "esriSFSSolid",
                      "color": [115,76,0,255],
                        "outline": {
                         "type": "esriSLS",
                         "style": "esriSLSSolid",
                         "color": [110,110,110,255],
                         "width": 1}
                    })
dc_map.draw(dc_df.spatial.bbox)

Get geometries from Feature objects

Besides the first approach (FeatureLayer -> SeDF -> geometry), geometries can also be extracted from FeatureSet or Feature objects (in short, this path can be represented as FeatureLayer -> FeatureSet -> Feature -> geometry objects).

This approach is similar to what we have done in previous section, such as pp_df = pd.DataFrame.spatial.from_layer(pp_fl) and pp_df.head().

Input
# tax incentives
pp_fl = FeatureLayer('https://maps2.dcgis.dc.gov/dcgis/rest/services/DCGIS_DATA/Administrative_Other_Boundaries_WebMercator/MapServer/24')
pp_fset = pp_fl.query()
pp_fset
Output
<FeatureSet> 98 features
Input
pp_sdf = pp_fset.sdf
pp_sdf.head()
Output
OBJECTID NAME WEB_URL Shape_Length Shape_Area FEDTRACTNO TRACTNO AREASQMI POPDENSITY TOTAL TRACT_NAME TRACT GEOID ACRES ID CITY STATUS SHAPE
0 1 Supermarket Tax Incentives https://dmped.dc.gov/page/supermarket-tax-ince... 3629.654865 7.690694e+05 19.01 19.1 0.296939 13524.645551 4016 2010 Tract 19.01 001901 11001001901 190.041207 15 Brightwood Qualified {"rings": [[[-8574620.698597332, 4717365.76244...
1 2 Supermarket Tax Incentives https://dmped.dc.gov/page/supermarket-tax-ince... 3936.332156 7.871354e+05 88.04 88.4 0.303915 7936.437428 2412 2010 Tract 88.04 008804 11001008804 194.505405 80 Ivy City Qualified {"rings": [[[-8569465.936656462, 4707937.83407...
2 3 Supermarket Tax Incentives https://dmped.dc.gov/page/supermarket-tax-ince... 2991.033699 2.994830e+05 89.03 89.3 0.115631 22770.708307 2633 2010 Tract 89.03 008903 11001008903 74.003848 87 Kingman Park Qualified {"rings": [[[-8569236.061578225, 4708392.59219...
3 4 Supermarket Tax Incentives https://dmped.dc.gov/page/supermarket-tax-ince... 3952.278822 8.896368e+05 10.02 10.2 0.343491 10020.649467 3442 2010 Tract 10.02 001002 11001001002 219.834053 107 Mclean Gardens Qualified {"rings": [[[-8579661.69146017, 4712166.487627...
4 5 Supermarket Tax Incentives https://dmped.dc.gov/page/supermarket-tax-ince... 4800.866903 1.077533e+06 34 34 0.416038 10448.569315 4347 2010 Tract 34 003400 11001003400 266.264201 91 Ledroit Park Qualified {"rings": [[[-8573522.865645228, 4711365.69090...
Input
for pp_feature in pp_fset:
    if pp_feature.attributes["POPDENSITY"]>55000:
        print(pp_feature.attributes)
{'OBJECTID': 92, 'NAME': 'Supermarket Tax Incentives', 'WEB_URL': 'https://dmped.dc.gov/page/supermarket-tax-incentives', 'Shape_Length': 1725.3222152283226, 'Shape_Area': 171893.27291680302, 'FEDTRACTNO': '28.01', 'TRACTNO': '28.1', 'AREASQMI': 0.06636841, 'POPDENSITY': 56849.33539918, 'TOTAL': 3773, 'TRACT_NAME': '2010 Tract 28.01', 'TRACT': '002801', 'GEOID': '11001002801', 'ACRES': 42.47578244, 'ID': 93, 'CITY': 'Lewis Subdivision', 'STATUS': 'Qualified'}
Input
print(type(pp_sdf.spatial))
display(pp_sdf.spatial)
<class 'arcgis.features.geo._accessor.GeoAccessor'>
Output

Because pp_sdf.iloc[0] is just one of the 98 features contained in the pp_sdf, its shape (as shown below) would only show as part of the shape previewed above.

Input
pp_sdf_row0 = pp_sdf.iloc[0].SHAPE
print(type(pp_sdf_row0))
display(pp_sdf_row0)
<class 'arcgis.geometry._types.Polygon'>
Output

Get geometries from a local layer

Previously, we have seen how geometries can be extracted from existing web feature layer, let's now take a look at how to read a local layer into a SeDF and get geometries from it.

To illustrate, we can save the SeDF into a local FeatureClass object (in the format of shapefile), then show how geometries can be extracted from the local layer:

Input
dc_shp_path = dc_df.spatial.to_featureclass("./dc_feature_class")
dc_shp_path
./dc_feature_class
Input
dc_df_from_fc = pd.DataFrame.spatial.from_featureclass(dc_shp_path)
dc_df_from_fc
Output
index objectid city_name state_city capital web_url areakm areamiles shape_leng shape_area SHAPE OBJECTID
0 0 1 Washington 1150000 Y http://www.dc.gov 177.47 68.52 67608.276922 1.774562e+08 {"rings": [[[-8584936.334474642, 4712272.26069... 0
Input
dc_df_from_fc.spatial
Output

Upon validation with the dc_df.spatial object, the shapes look identical. Or we can check its geometry as represented in JSON:

Input
dc_df_from_fc["SHAPE"]
Output
0    {"rings": [[[-8584936.334474642, 4712272.26069...
Name: SHAPE, dtype: geometry

Construct higher level objects from geometry objects

Now we have learnt how to extract geometries from FeatureLayer, FeatureSet, Feature objects, and local feature layer, let's turn the workflow upside down, and focus at composing Feature, FeatureSet, and SeDF objects from geometry objects:

Create Feature, FeatureSet and FeatureCollection objects from Geometry objects

With geometry objects (point or line or polygon) created, how do we create a Feature, and FeatureSet objects out of them? Let's first look at an example of two Point geometries created from scratch:

Input
pt_1 = Point({'x': -101.376, 'y': 31.119})
fire_1 = Feature(geometry=pt_1, attributes={"id": "EONET_4959",
                                            "title": "Ferguson Fire",
                                            "description": "",
                                            "link": "https://eonet.sci.gsfc.nasa.gov/api/v2.1/events/EONET_4959",
                                            "categories": {
                                                "id": 8,
                                                "title": "Wildfires"
                                            }})
fire_1
Output
{"geometry": {"x": -101.376, "y": 31.119}, "attributes": {"id": "EONET_4959", "title": "Ferguson Fire", "description": "", "link": "https://eonet.sci.gsfc.nasa.gov/api/v2.1/events/EONET_4959", "categories": {"id": 8, "title": "Wildfires"}}}
Input
pt_2 = Point({'x': -101.91, 'y': 30.859999999999999})
fire_2 = Feature(geometry=pt_2, attributes={"id": "EONET_4958",
                                            "title": "Smith Canyon Fire",
                                            "description": "",
                                            "link": "https://eonet.sci.gsfc.nasa.gov/api/v2.1/events/EONET_4958",
                                            "categories": {
                                                "id": 8,
                                                "title": "Wildfires"
                                            }})
fire_2
Output
{"geometry": {"x": -101.91, "y": 30.86}, "attributes": {"id": "EONET_4958", "title": "Smith Canyon Fire", "description": "", "link": "https://eonet.sci.gsfc.nasa.gov/api/v2.1/events/EONET_4958", "categories": {"id": 8, "title": "Wildfires"}}}
Input
fires_fset = FeatureSet(features = [fire_1, fire_2], 
                        geometry_type="Point", 
                        spatial_reference={'latestWkid': 4326, 'wkid': 102100})
fires_fset
Output
<FeatureSet> 2 features
Input
fires_fset.sdf
Output
id title description link categories OBJECTID SHAPE
0 EONET_4959 Ferguson Fire https://eonet.sci.gsfc.nasa.gov/api/v2.1/event... {'id': 8, 'title': 'Wildfires'} 1 {"type": "Point", "coordinates": [-101.376, 31...
1 EONET_4958 Smith Canyon Fire https://eonet.sci.gsfc.nasa.gov/api/v2.1/event... {'id': 8, 'title': 'Wildfires'} 2 {"type": "Point", "coordinates": [-101.91, 30....
Input
fires_fc = FeatureCollection.from_featureset(fset=fires_fset)
fires_fc
Output
<FeatureCollection>

Create an SeDF object with geometries

Besides Feature, FeatureSet, and FeatureCollection, SeDF objects can also be created with geometry objects.

Input
dc_recycle_geom1_dict = {"x":-77.1188753,"y":38.8614258,"spatialReference":{"wkid":4326}}
dc_recycle_geom2_dict = {"x":-77.08481,"y":38.8478148,"spatialReference":{"wkid":4326}}
Input
dc_recycle_geom1 = Geometry(dc_recycle_geom1_dict)
display(dc_recycle_geom1)
Input
dc_recycle_geom2 = Geometry(dc_recycle_geom2_dict)
display(dc_recycle_geom2)
Input
dc_recycle_dict = {"osm_id":{"0":"2607287308","1":"3752101908"},
                   "amenity":{"0":"recycling","1":"recycling"},
                   "source":{"0":"site visit","1":"GPSO"},
                   "recycling_type":{"0":"container","1":"container"}
                   }
dc_recycle_df = pd.DataFrame.from_dict(dc_recycle_dict)
dc_recycle_df.head()
Output
osm_id amenity source recycling_type
0 2607287308 recycling site visit container
1 3752101908 recycling GPSO container
Input
dc_recycle_df["X"] = [dc_recycle_geom1.x, dc_recycle_geom2.x]
dc_recycle_df["Y"] = [dc_recycle_geom1.y, dc_recycle_geom2.y]
dc_recycle_df.head()
Output
osm_id amenity source recycling_type X Y
0 2607287308 recycling site visit container -77.118875 38.861426
1 3752101908 recycling GPSO container -77.084810 38.847815
Input
dc_recycle_sedf = GeoAccessor.from_xy(dc_recycle_df, x_column="X", y_column="Y", sr=4326)
dc_recycle_sedf.head()
Output
osm_id amenity source recycling_type X Y SHAPE
0 2607287308 recycling site visit container -77.118875 38.861426 {"spatialReference": {"wkid": 4326}, "x": -77....
1 3752101908 recycling GPSO container -77.084810 38.847815 {"spatialReference": {"wkid": 4326}, "x": -77....
Input
dc_recycle_sedf.spatial
Output

The dc_recycle_sedf object would be showing a shape of rectangle in the previous cell if you are currently on arcpy or shapely geometry engines. If not, then you would see its JSON representation instead.

Input
dc_recycle_sedf.SHAPE
Output
0    {"spatialReference": {"wkid": 4326}, "x": -77....
1    {"spatialReference": {"wkid": 4326}, "x": -77....
Name: SHAPE, dtype: geometry

Conclusion

In this part of the guide series to arcgis.geometry module, you have seen the introduction to the module, how to create Geometry objects, what are the basic properties, and how to work with one, including its interactions with map widgets. Also, it explores how geometries can be extracted from existing web layers or local feature classes, and also provides introduction to how Feature, FeatureSet, FeatureCollection and SeDF objects can be constructed from geometries.

Next, part 3 will discuss the spatial operations of geometry objects, and how to process through geometry services, and part 4 will talk about geometry filters and how they being applied in query, search and mapping.

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.