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.

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.

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.

To get the geometry type, use the type property:

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

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.

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.

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.

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.

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.

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

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:

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.

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

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).

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.

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.

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

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

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.

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:

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:

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.

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.

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

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.

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).

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.

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:

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().

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.

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:

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

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:

Create an SeDF object with geometries¶

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

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.

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.