I am attempting to access ArcGIS Map Services and query a dataset to only return records that intersect the boundary of a custom geometry I created. I have queried and downloaded data before using this script, but never using geometry. Do i simply add my geometry_json to the Where clause?
#My geometry (a simply polygon)
geometry_json = {
"geometryType" : "esriGeometryPolygon",
"spatialReference" : {
"wkid" : 4326,
"latestWkid" : 4326
},
"features" : [{
"geometry":{
"curveRings":[
[[-123.77335021899995,49.353870065000081],{"a":[[-123.77335021899995,49.353870065000081],[-123.88831213030559,49.338864996255367],0,1]}]]}}]
}
The script to access and download the data
import arcpy # Import the arcpy module, which provides access to ArcGIS functionality.
import urllib.request # Import the urllib.request module to work with URLs.
import json # Import the json module for JSON data handling.
# Setup
arcpy.env.overwriteOutput = True
baseURL = r"https://gisp.dfo-mpo.gc.ca/arcgis/rest/services/FGP/MSDI_Dynamic_Current_Layer/MapServer/0" # Base URL
fields = "*" # Return all fields
outdata = r"out/path" #output location
# Get record extract limit
urlstring = baseURL + "?f=json" # Construct the URL to request JSON metadata from the service.
j = urllib.request.urlopen(urlstring) # Open the URL and get the response.
js = json.load(j) # Load the JSON response into a Python dictionary.
maxrc = int(js["maxRecordCount"]) # Extract the maximum record count from the JSON response.
print ("Record extract limit: %s" % maxrc) # Print the maximum record count.
# Get object ids of features
where = "1=1" # Define a where clause to retrieve all records.
urlstring = baseURL + "/query?where={}&returnIdsOnly=true&f=json".format(where) # Construct the URL to request object IDs.
j = urllib.request.urlopen(urlstring) # Open the URL and get the response.
js = json.load(j) # Load the JSON response into a Python dictionary.
idfield = js["objectIdFieldName"] # Extract the name of the object ID field.
idlist = js["objectIds"] # Extract the list of object IDs.
idlist.sort() # Sort the list of object IDs.
numrec = len(idlist) # Get the number of records.
print ("Number of target records: %s" % numrec) # Print the number of records.
# Gather features
print ("Gathering records...") # Print a message indicating feature gathering.
fs = dict() # Create an empty dictionary to store features.
for i in range(0, numrec, maxrc): # Loop over the object IDs in chunks based on the maximum record count.
torec = i + (maxrc - 1) # Calculate the end index of the chunk.
if torec > numrec: # Adjust the end index if it exceeds the number of records.
torec = numrec - 1
fromid = idlist[i] # Get the starting object ID of the chunk.
toid = idlist[torec] # Get the ending object ID of the chunk.
where = "{} >= {} and {} <= {}".format(idfield, fromid, idfield, toid) # Define a where clause for the chunk.
print (" {}".format(where)) # Print the where clause.
urlstring = baseURL + "/query?where={}&returnGeometry=true&outFields={}&f=json".format(where,fields) # Construct the URL to request features.
fs[i] = arcpy.FeatureSet() # Create a new empty FeatureSet object.
fs[i].load(urlstring) # Load features from the URL into the FeatureSet.
# Save features
print ("Saving features...") # Print a message indicating feature saving.
fslist = [] # Create an empty list to store FeatureSets.
for key,value in fs.items(): # Iterate over the dictionary of FeatureSets.
fslist.append(value) # Append each FeatureSet to the list.
arcpy.Merge_management(fslist, outdata) # Merge all FeatureSets into a single feature class at the specified output location.
print ("Done!") # Print a message indicating completion.
