The ArcGIS Python API also provides the functionality to access individual features in a feature layer, to manipulate layers, and to run analyses from a large collection of GIS tools including typical GIS geoprocessing tools. Most of this functionality is available in the different submodules of arcgis.features. To give you a few examples, let’s continue with the city feature service we still have in variable neCitiesFS from the previous section. A feature service can actually have several layers but, in this case, it only contains one layer with the point features for the cities. The following code example accesses this layer (neCitiesFE.layers[0]) and prints out the names of the attribute fields of the layer by looping through the ...properties.fields list of the layer:
for f in neCitiesFS.layers[0].properties.fields:
print(f)
When you try this out, you should get a list of dictionaries, each describing a field of the layer in detail including name, type, and other information. For instance, this is the description of the STATEABB field:
{
"name": "STATEABB",
"type": "esriFieldTypeString",
"actualType": "nvarchar",
"alias": "STATEABB",
"sqlType": "sqlTypeNVarchar",
"length": 10,
"nullable": true,
"editable": true,
"domain": null,
"defaultValue": null
}
The query(...) method allows us to create a subset of the features in a layer using an SQL query string similar to what we use for selection by attribute in ArcGIS itself or with arcpy. The result is a feature set in ESRI terms. If you look at the output of the following command, you will see that this class stores a JSON representation of the features in the result set. Let’s use query(...) to create a feature set of all the cities that are located in Pennsylvania using the query string "STATEABB"='US-PA' for the where parameter of query(...):
paCities = neCitiesFS.layers[0].query(where='"STATEABB"=\'US-PA\'')
print(paCities.features)
paCities
Output:
[{"geometry": {"x": -8425181.625237303, "y": 5075313.651659228},
"attributes": {"FID": 50, "OBJECTID": "149", "UIDENT": 87507, "POPCLASS": 2, "NAME":
"Scranton", "CAPITAL": -1, "STATEABB": "US-PA", "COUNTRY": "USA"}}, {"geometry":
{"x": -8912583.489066456, "y": 5176670.443556941}, "attributes": {"FID": 53,
"OBJECTID": "156", "UIDENT": 88707, "POPCLASS": 3, "NAME": "Erie", "CAPITAL": -1,
"STATEABB": "US-PA", "COUNTRY": "USA"}}, ... ]
<FeatureSet> 10 features
Of course, the queries we can use with the query(...) function can be much more complex and logically connect different conditions. The attributes of the features in our result set are stored in a dictionary called attributes for each of the features. The following loop goes through the features in our result set and prints out their name (f.attributes['NAME']) and state (f.attributes['STATEABB']) to verify that we only have cities for Pennsylvania now:
for f in paCities:
print(f.attributes['NAME'] + " - "+ f.attributes['STATEABB'])
Output: Scranton - US-PA Erie - US-PA Wilkes-Barre - US-PA ... Pittsburgh - US-PA
Now, to briefly illustrate some of the geoprocessing functions in the ArcGIS Python API, let’s look at the example of how one would determine the parts of the Appalachian Trail that are within 30 miles of a larger Pennsylvanian city. For this we first need to upload and publish another data set, namely one that’s representing the Appalachian Trail. We use a data set that we acquired from PASDA for this, and you can download the DCNR_apptrail file here. As with the cities layer earlier, there is a very small chance that this will not display for you - but you should be able to continue to perform the buffer and intersection operations. The following code uploads the file to AGOL (don’t forget to adapt the path and add your initials or ID to the name!), publishes it, and adds the resulting trail feature service to your cityMap from above:
# change this variable to the name of the zip file and service
dcnr_file = "dcnr_apptrail_2003_???"
appTrail = rf'C:\PSU\Geog 489\Lesson 3\Notebooks\{dcnr_file}.zip'
# test if it exists in AGOL first
appTrailshp_test = gis.content.search(query=f'owner:{user_name} title:{dcnr_file}', item_type='Feature Layer Collection')
if appTrailshp_test:
appTrailFS = appTrailshp_test[0]
print(f"{appTrailFS.title} found in AGOL")
else:
root_folder = gis.content.folders.get()
# Ensure the folder is found
if root_folder:
# Add the shapefile to the folder.
appTrailshp = root_folder.add({"title": dcnr_file, "type": "Shapefile"}, appTrail).result()
# Print the result
print(f"Shapefile added to folder: {appTrailshp.title}")
appTrailFS = appTrailshp.publish()
print(f"{appTrailFS.title} published as a service")
else:
print("Folder not found")
display(appTrailFS)
cityMap.content.add(appTrailFS, {})
cityMap
Let's take a moment to talk about credits. Credits are what esri uses to charge users for certain processing tasks and storage in AGOL. Processes such as geocoding, clipping, routing, buffering, or storing data and content costs the organization credits. More details can be found in esri's documentation here, and it is important to protect your account credentials so someone does not drain your credits (among other things). The rest of the walkthrough performs tasks that costs credits, so we include preliminary code that checks if the result exist in AGOL before executing the geoprocessing step. When you first run the blocks below, it will create the data in AGOL since it most likely does not exist. Subsequent executions (who executes a script just once?) it will find and use the existing data from a previous run, saving your credits.
Next, we create a 30 miles buffer around the cities in variable paCities that we can then intersect with the Appalachian Trail layer. The create_buffers(...) function that we will use for this is defined in the arcgis.features.use_proximity module together with other proximity based analysis functions. We provide the feature set we want to buffer as the first parameter, but, since the function cannot be applied to a feature set directly, we have to invoke the to_dict() method of the feature set first. The second parameter is a list of buffer distances allowing for creating multiple buffers at different distances. We only use one distance value, namely 30, here and also specify that the unit is supposed to be ‘Miles’. Finally, we add the resulting buffer feature set to the cityMap above.
This operation costs credits so once it is created in AGOL, try not to delete it unless absolutely necessary.
The title is used for the output_name so it can be found in AGOL after initial creation.
#
from arcgis.features.use_proximity import create_buffers
import urllib3
title = 'bufferedCities_???'
# test if it exists in AGOL first
bufferedCities_test = gis.content.search(query=f'owner:{user_name} title:{title}', item_type='Feature Layer Collection')
if bufferedCities_test:
bufferedCities = bufferedCities_test[0]
print(f"{bufferedCities.title} found in AGOL")
else:
# Disable SSL warnings (or it will print a lot of them)
urllib3.disable_warnings(urllib3.exceptions.InsecureRequestWarning)
bufferedCities = create_buffers(paCities.to_dict(), [30], units='Miles', output_name=title)
# Print the result
print(f"Shapefile added to folder: {bufferedCities.title}")
cityMap.content.add(bufferedCities, {})
As the last step, we use the overlay_layers(...) function defined in the arcgis.features.manage_data module to create a new feature set by intersecting the buffers with the Appalachian Trail polylines. For this we have to provide the two input sets as parameters and specify that the overlay operation should be ‘Intersect’.
This operation costs credits so once it is created in AGOL, try not delete it unless absolutely necessary.
The title is used for the output_name so it can be found in AGOL after initial creation.
#
from arcgis.features.analysis import overlay_layers
# OR from arcgis.features.manage_data import overlay_layers which does the same
title = 'city_trailhead_???'
# test if it exists in AGOL first
trailPartsCloseToCity_test = gis.content.search(query=f'owner:{user_name} title:{title}', item_type='Feature Layer Collection')
if trailPartsCloseToCity_test:
trailPartsCloseToCity = trailPartsCloseToCity_test[0]
print(f"{trailPartsCloseToCity.title} found in AGOL")
else:
# Perform the overlay (Intersection)
trailPartsCloseToCity = overlay_layers(appTrailFS.layers[0], bufferedCities.layers[0], overlay_type='Intersect', output_name=title)
print(f"{trailPartsCloseToCity.title} created")
# Check the result
print(trailPartsCloseToCity)
We show the result by creating a new map ...
#
resultMap = gis.map('Pennsylvania')
... and just add the features from the resulting trailPartsCloseToCity feature set to this map:
#
resultMap.content.add(trailPartsCloseToCity, {})
resultMap
The result is shown in the figure below.
