Extracting MORE Features from Map Services

Back in August 2015 I posted information about how you can extract features from a map service.  Since then, I have had many contact me about modifying the code so it can extract features beyond the record limit set in the map service.  So today I decided to work on one that does!

To test the script, I headed over to the map services provided by the State of California GIS.  Specifically, the one for wildfires:

http://services.gis.ca.gov/arcgis/rest/services/Environment/Wildfires/MapServer

When you scroll down a quick look will reveal that there is a maximum record count of 1000:

extractpt2-1

At the top, click on the “ArcGIS Online map viewer” link to take a look at the data.  I changed my basemap to Dark Gray Canvas so the symbology really pops.  Here is the area around Santa Clarita:

extractpt2-2

What a mess!  There sure has been a lot of fires in this area.  Clicking on the Show Table icon for the Wildfires layer will reveal how many records (wildfire polygons) I will need to extract:

extractpt2-3

I will extract all 17,550 polygons with the following Python script:

import arcpy
import urllib2
import json

# Setup
arcpy.env.overwriteOutput = True
baseURL = "http://services.gis.ca.gov/arcgis/rest/services/Environment/Wildfires/MapServer/0"
fields = "*"
outdata = "H:/cal_data/data.gdb/testdata"

# Get record extract limit
urlstring = baseURL + "?f=json"
j = urllib2.urlopen(urlstring)
js = json.load(j)
maxrc = int(js["maxRecordCount"])
print "Record extract limit: %s" % maxrc

# Get object ids of features
where = "1=1"
urlstring = baseURL + "/query?where={}&returnIdsOnly=true&f=json".format(where)
j = urllib2.urlopen(urlstring)
js = json.load(j)
idfield = js["objectIdFieldName"]
idlist = js["objectIds"]
idlist.sort()
numrec = len(idlist)
print "Number of target records: %s" % numrec

# Gather features
print "Gathering records..."
fs = dict()
for i in range(0, numrec, maxrc):
  torec = i + (maxrc - 1)
  if torec > numrec:
    torec = numrec - 1
  fromid = idlist[i]
  toid = idlist[torec]
  where = "{} >= {} and {} <= {}".format(idfield, fromid, idfield, toid)
  print "  {}".format(where)
  urlstring = baseURL + "/query?where={}&returnGeometry=true&outFields={}&f=json".format(where,fields)
  fs[i] = arcpy.FeatureSet()
  fs[i].load(urlstring)

# Save features
print "Saving features..."
fslist = []
for key,value in fs.items():
  fslist.append(value)
arcpy.Merge_management(fslist, outdata)
print "Done!"

I put in a few comments so you could see how things are done where.  This script also gives you some feedback as it runs.  I will explain the script.

First, it loads the arcpy module as well as the urllib2 and json modules (we need those to open the url and read in the json output).  In the Setup section, I set my baseURL variable to the wildfire map service layer 0 (you would change this line for the map service layer you are interested in), the fields variable to “*” to grab all the attributes, and the outdata variable to where I want to store my features, in this case in a file geodatabase and a feature class name of testdata (change this line to where you want to save your data, you can also specify a shapefile name instead).

Next, the script extracts the record limit for the map service and saves it to the maxrc variable.  Next, the script grabs the id field name (which in this case is OBJECTID) and all the record id numbers are stored in the idlist variable.  The id list is sorted and then the total number of records value is stored in the numrec variable.

Next, all the feature records are gathered up.  This is done by taking the total number of records and stepping through them in chunks defined by the record limit.  Since this map service has a limit of 1000, the first 1000 records are selected by the id field (OBJECTID) and saved to a temporary featureset fs[i].  It then cycles to the next 1000 records and stores those features, etc., until it reaches that last set.  Note the if statement “if torec > numrec” is a fail safe when we reach the end and our to record number we calculate is higher than the maximum number of records, and if so, the torec value is set to the last record number in the list.  After that, the loop is finished.

Finally, we loop through all the temporary featuresets to put them in a list and then pass that list to merge them all together to create our outdata.  Whew!

Here is what my output looked like when I ran the Python script:

extractpt2-4

Note that as it gathers the features it shows you the queries it uses to do it, in 1000 record chunks.  Also note that sometimes there are missing record numbers, probably because those records with certain OBJECTIDs were either edited and changed numbers or were deleted.  So that is why you see in the first grab OBJECTIDs from 1 to 1003.  It is not grabbing 1003 records, it is just that the first 1000 record numbers in the list start with OBJECTID 1 and end in OBJECTID 1003.  Get it?

The script took a little over 3 minutes to run.  Taking a look using ArcCatalog shows the gathered features:

extractpt2-5

And looking at the table shows I grabbed all 17,550 records and the attributes:

extractpt2-6

I hope this updated script will help you out with those map services that have record limits.  Enjoy!  -mike

5 thoughts on “Extracting MORE Features from Map Services

  1. This is good stuff Mike. How does the python script work with an app? In what way is it incorporated? Do you have any samples?

    Thank you
    Brandon Price

    • Hi Brandon. This python script does use the arcpy module, so it was written for use with ArcGIS. I have a Windows 7 machine, so I just run it at the DOS prompt. You could also just double click the script and it will run too. You could modify the script to not use arcpy at all, you would just need to load other modules that would allow you to save to a shapefile. I’m not sure if this script is app worthy, but I guess you could create one that would allow the user to input a map service URL and output format. I know you could turn this into a geoprocessing tool in ArcGIS with a few modifications. -mike

      • Hey Mike:

        Thanks for the clarification. This is still a big benefit for me for ArcGIS Desktop. I was just wondering because I have an app that I might want to add a download widget to later on down the line.

        Brandon

  2. Thank you so much for sharing this! We are moving all of our services to AGOL, but have needs for local geoprocessing. This is just what I needed. I had to do some tweaks to get it to run using ArcGIS Pro’s Python 3 installation. I’d be happy to share those modifications if you’d like.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s