在 Python 中使用几何图形进行 ArcGIS REST 地图服务数据查询

Question

我正在尝试访问 ArcGIS 地图服务并查询数据集，以仅返回与我创建的自定义几何图形边界相交的记录。在使用这个脚本之前我已经查询并下载了数据，但从未使用过几何。我是否只需将我的

geometry_json

添加到

Where

子句中？

#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]}]]}}]
}

访问和下载数据的脚本

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.

Answer 1

与 ArcGIS REST API 交互时，除了

arcgis

之外，您还可以使用

arcpy

包来节省大量时间/精力/错误处理。

arcpy

是 ArcMap/ArcGIS Pro 的界面，而

arcgis

是 ArcGIS REST API 的界面（无论是 ArcGIS Online 还是 ArcGIS Enterprise）。

这是一个简短的例子：

from arcgis.gis import GIS
from arcgis.geometry.filters import intersects

gis = GIS("pro")

geometry = {
    "rings" : [[
        [-117.751322686672, 34.1209151280806], 
        [-117.751320675015, 34.1209101319963], 
        [-117.751322686672, 34.1209068012733], 
        [-117.751326709986, 34.1209056910323], 
        [-117.751330733299, 34.1209062461528], 
        [-117.751333415508, 34.1209095768758], 
        [-117.751331403852, 34.1209151280806], 
        [-117.751327380538, 34.1209173485624], 
        [-117.751322686672, 34.1209151280806]
    ]]
}

item = gis.content.get("<item_id>")
layer = item.layers[<layer_id>]
features = layer.query(where="1=1", geometry_filter=intersects(geometry, 4326))

如果使用 ArcGIS Pro 的 Python 环境，您可以通过

gis = GIS("pro")

使用 ArcGIS Pro 中的活动门户对 GIS 对象进行身份验证。有关其他登录选项，请参阅 doc。