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

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

  3. Hi Micheal,
    Thank you so much for all your efforts on this code, super cool! Unfortunately I am running into an issue right at the end of the script, I receive the following error: “Runtime error Traceback (most recent call last): File “”, line 49, in File “c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py”, line 4217, in Merge raise e ExecuteError: ERROR 999998: Unexpected Error.”

    Line 49 in the code is “arcpy.Merge_management(fslist, outdata)”
    And line 4217 in my management.py file is “raise e” (I can provide the previous lines if that’s helpful to you).
    I know it’s likely impossible to troubleshoot all errors one could expect but I’m not sure if something jumps out at you. I’m not python expert so I appreciate any help you can offer.
    I should also note, I copied and pasted your code exactly as you have it, within the python window in ArcGIS 10.3; only changing the path to the map service and where I desired to save my output data.
    Thank you again!
    ~D

    • Hi Danielle. Either there is an anomaly in the map service data, or maybe the name you use for the outdata variable is not valid by ArcGIS standards. For example, does you data name have spaces in it? Sometimes that is a problem. Is it going to a .SHP file? If so, is it hitting some limits with the number of features or fields? Send me the map service path you are using and I can give it a try as long as I have access to it. -mike

      • Another thing that comes to mind is your computer. Maybe the environment on your computer is causing a problem. Does it have enough memory for all those temporary feature datasets before the merge? Is your hard drive fast enough? Have you tried using another computer to see if the code works? Have you tried upgrading ArcGIS Desktop to see if the problem goes away? I am using ArcGIS Desktop 10.5.1. Just some thoughts. -mike

      • Mike,
        Thank you much! This is the service I’m using: https://gis.nassaucountyny.gov/arcgis/rest/services/ncgis/MapServer/13

        I’m fairly certain it’s not locked down as I’m able to run a single query and download batches of a 1000 but there any many records (407,888) so doing this action 400+ times is prohibitive and makes your script all the more appealing 🙂

        I’ve been sure to exclude any spaces in my naming conventions. I’ve attempted another, smaller map service as well but with the same result. Additionally I also attempted using your map service and simply changed the output path (again no spaces) to a geodatabase of the same name; still no luck there.

        I also used another PC, as you suggest below, same result, but it’s also on the same computer network with the same ArcGIS build specs.
        As far as my PC; I have a 64-bit OS with 8GB of RAM (7.88GB usable) and a 3.40GHz processor.

        ~D

      • Ok, so I think I know the problem … too much data! With 407,888 records, that means you will need enough memory to store it all, then basically merge all 408 temporary feature datasets into one. I don’t think you have enough memory OR merge cannot handle it all. I am going to try on my end with a hefty computer and see what happens to me!

      • Well it worked for me, go figure! I had 12gb of memory, but noticed it only consumed about half a gig when the script was almost done. I am using ArcGIS 10.5.1, so maybe that helped me? I saved the parcels to a geodatabase and it is only 88mb in size. Anyway, if you want the data I downloaded, feel free to email me at mcarson@burbankca.gov and I will send you a link.

      • Well imagine that! I appreciate you giving it a whirl! I definitely think it’s the Desktop versioning since I’ve tried map services with far fewer records and had the same result.
        I’ll download the latest version and report back just in case others have the same issue.
        Thank you for all of your efforts Mike!

  4. I successfully updated my ArcGIS Desktop Version to 10.6 and still no luck with running the script. Wanted to let folks know if they run into the same error I posted above.

    • Hi Danielle. I was thinking about your issue with the script. I still think it is a memory issue for you. Do you run the script inside ArcMap using the Python window? If so, try running the python script at the Windows DOS prompt instead (Windows cmd window).

      Another option is the script below. It will save each feature extract instead of using memory, then it merges all the features together. Just modify the outdata variable and see if that works. Do note that it will leave a bunch of temporary feature classes in the geodatabase, each representing the 1000 feature extract. Give it a try.

      import arcpy
      import urllib2
      import json

      # Setup
      arcpy.env.overwriteOutput = True
      baseURL = “https://gis.nassaucountyny.gov/arcgis/rest/services/ncgis/MapServer/13”
      fields = “*”
      outdata = “C:/temp/data.gdb/ncgis_parcels”

      # 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…”
      merge_list = []
      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 = arcpy.FeatureSet()
        fs.load(urlstring)
        tempdata = outdata + str(i)
        print " Copying features…"
        arcpy.CopyFeatures_management(fs, tempdata)
        merge_list.append(tempdata)

      # Save features
      print "Saving features…"
      arcpy.Merge_management(merge_list, outdata)
      print "Done!"

  5. We were having the same issue. Have you tried turning off the background processing in Geoprocessing Tab>Geoprocessing Options and then uncheck enable background processing.

    Oh and thanks Mike that script helped a ton!

    • Thanks Bryan. I have background geoprocessing on and it works for me, I tested it. However, I do not run the script inside ArcMap using the Python window. I run it at the Windows DOS prompt. If Python is setup right in Windows, you can even double click the .py script and it will run too.

      • Mike, thanks for the additional info. I will definitely keep that in mind. My comment above was more in reference to possibly help Danielle out. I was not detailed enough in my previous response.

        I was getting the same error Danielle had running the script through the python window in ArcCatalog which was the “Runtime error Traceback (most recent call last): File “”, line 49…..: ERROR 999998: Unexpected Error.”” But as soon as I disabled background processing (Geoprocessing Tab>Geoprocessing Options and then uncheck enable background processing) the script you provided worked exactly how it is supposed to.

        Again, thank you Mike for sharing the script with us.

    • Hi Kegan. It is not too difficult, you just need to authenticate yourself with ArcGIS Server and get a token, then you use that token when you make HTTPS requests to secure map and feature services.

      For example, to manually do this if an ArcGIS Server is setup for tokens (see http://maps.scag.ca.gov/scaggis/rest/services) you will see in the upper right corner a “Get Token” link. If you click on that, you will be taken to the tokens page (https://maps.scag.ca.gov/scaggis/tokens/) to enter a username and password to get a token from the server. Note you can specify the Client (I use RequestIP to limit my token to my current IP address), Expiration (I usually keep it to an hour), and Format (I use JSON). Once you get the token, you have access to the secured services at that site.

      So how to do this in Python? First at the top of your code where you import the arcpy module, you need to import the requests module.

      import requests

      Then the code block to get the token looks like this:

      tokenurl = "https://YOUR_ARCGIS_SERVER/tokens/"
      username = "YOUR_USERNAME"
      password = "YOUR_PASSWORD"
      params = {'username':username, 'password':password, 'client':'requestip', 'expiration':60, 'f':'json'}
      tokenresponse = requests.post(tokenurl, data=params, verify=False)
      data = tokenresponse.json()
      token = str(data["token"])

      Note for the params I set the client to “requestip”, the expiration to 60 (in minutes, so 1 hour), and my output format to “json”.

      With the token now set in the token variable, you use that variable when you make requests to secure map and feature services. You do this by adding “&token=” to the calling URL with the token value after the equal sign. For example, using the code sample in this post, I would change the urlstring variable and add the token parameter:

      urlstring = baseURL + "?f=json&token={}".format(token)

      and

      urlstring = baseURL + "/query?where={}&returnIdsOnly=true&f=json&token={}".format(where,token)

      and also

      urlstring = baseURL + "/query?where={}&returnGeometry=true&outFields={}&f=json&token={}".format(where,fields,token)

      That’s it!

      • Hi Michael

        Thank you so much for the advice! I had to use Webapp means of authentification to generate a token (the requestIP generated token did not seem to work for some reason). I eventually got the process to work since my dataset had about 14 million verticies it kept giving me the following error:

        “RuntimeError: RecordSetObject: Cannot open table for Load”,

        Turns out i had to reduce my MaxRecordCount to increments of 500 and it worked!

        Keep up the good work
        Kegan

  6. I keep getting an error when I run this code in ArcGIS Pro:
    import arcpy
    import urllib2
    import json

    # Setup
    arcpy.env.overwriteOutput = True
    baseURL = “https://gis2.idaho.gov/arcgis/rest/services/DPR/Idaho_Trails_MapC/MapServer/4/query”
    fields = “NAME,OFFICE,HUNTING_UNIT”
    outdata = “Q:\Temp\Mark_Temp\Idaho_Project\Idaho Trails Map\Trailheads.gdb/Automobiles2”

    # 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 get a error message:
    File "”, line 16
    print “Record extract limit: %s” % maxrc
    ^
    SyntaxError: Missing parentheses in call to ‘print’. Did you mean print(“Record extract limit: %s” % maxrc)?

    I’ll change these out but will still give me a syntax error.

    When I run the same code in ArcMap 10.3.1 I recieve this error:
    Runtime error
    Traceback (most recent call last):
    File “”, line 15, in
    KeyError: ‘maxRecordCount’

    I’m not sure what I’m doing wrong, any advice?

    • Hi Mark. The code was made to run in Python that comes with ArcGIS Desktop (ArcMap), not ArcGIS Pro. ArcGIS Pro uses Python 3.6.x while ArcGIS Desktop uses 2.7.x, so the code would need to be modified to work in Python 3 and ArcGIS Pro. As for the error you are getting using ArcMap, I would have to give it a try and see what is going on. -mike

    • Mark, remove the “/query” from your baseURL variable. It is messing up the code when it tries to query the map service to get the max record count. If you want to run in ArcGIS Pro, you will need to modify all the print statements to the new print function print( ). -mike

  7. Hi Michael,

    Thank you so much for the two posts to extract features from a map services.

    I’ve been trying to get the information from a mapserver with Python 2.7.3 from ArcGIS 10.2 but I’m getting an error when I run it:

    Record extract limit: 1000
    Number of target records: 11175
    Gathering records…
    ESRI_OID >= 1 and ESRI_OID <= 1000

    Traceback (most recent call last):
    File "C:/Users/ecarrillo/Downloads/Extraer_IGME_BDMIN_Indicios_II_(ID-0)_bis.py", line 48, in
    fs.load(urlstring)
    File “C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\arcobjects\arcobjects.py”, line 173, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))
    RuntimeError: RecordSetObject: Cannot open table for Load

    How can I solve it.

    Thanks,

    Eloy

      • Hi Michael,

        Thanks for the quick answer. M baseURL is:
        http://mapas.igme.es/gis/rest/services/Cartografia_Tematica/IGME_BDMIN_Indicios/MapServer/1
        I don’t think the map service is secure because normaly, I’ve applied the algorithm of the previous post to extract the first 1000 values, and it’s work perfectly. Then for extract the other values I used to change the line where=”1=1″ by where=”objectid>=1000+and+objectid<1999&f=json" and I subsequently varying from 1000 to 1000 until have all the elements. The ol roblem it's when you have a lot of elements.

        Thanks,

        Eloy

      • Eloy, it looks like that server is very slow at times. I get a “504 Gateway time-out” error. When that happens no features are returned and the python script errors out. There is nothing much you can do unless they fix their server’s response to internet traffic. -mike

      • Eloy, since I do not like to let the computer win, after some experimenting, I was able to get the script to work with that map service. You will need to make the following modifications for it to work.

        First, change the maxrc variable to 50 and comment out the rest:

        maxrc = 50 #int(js["maxRecordCount"])

        I found gathering 50 records at a time from the map service instead of 1000 helped the map service keep up with little delay.

        Next, you need to add an oidlist variable and oids variable to the code, as well as modify the urlstring variable, all in the “Gather features” code block. Here is what I did:

        # 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)
          oidlist = idlist[i:torec + 1]
          print "  {}".format(where)
          oids = ','.join(str(e) for e in oidlist)
          print oids
          urlstring = baseURL + "/query?objectIds={}&returnGeometry=true&outFields={}&f=json".format(oids,fields)
          fs[i] = arcpy.FeatureSet()
          fs[i].load(urlstring)

        Instead of selecting records using a range of ids in a field, you specify the exact ids to pull out separated by a comma (like 1,2,3,4,5,6). Note the urlstring variable was changed to query the objectIds parameter instead of using the where parameter in the original code. For some reason this map service hates pulling records using the where parameter, but works using the objectIds. I also included a print statement so you can see the value of the oids variable when it is passed to the map service query.

        Give this a try. It worked for me. Since this pulls out by 50 record chunks, it does take a while to run, but you get all 11175 point features. -mike

      • Hi Michael,

        Thanks you very much, the script works perfectly.

        Keep up the good work and thanks you again

        Eloy

      • Hi, Michael, im getting other error after modify the script.

        Traceback (most recent call last):
        File “C:\Python27\ArcGIS10.5\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py”, line 326, in RunScript
        exec codeObject in __main__.__dict__
        File “D:\JuanJulio\Cursos\PennState\GetFeatures_v2.py”, line 23, in
        idfield = js[“objectIdFieldName”]
        KeyError: ‘objectIdFieldName’

        This is the feature service im trying to download:

        http://geocatmin.ingemmet.gob.pe/arcgis/rest/services/SERV_GEOLOGIA_REGIONAL/MapServer/5

        I hope you could help me. Maybe im missing something in the script.

        Thank you.

      • Hi Juan. That map service is pretty slow, and I get the error too. You can take a look at the layer in the map service using the ArcGIS Online map viewer:

        http://www.arcgis.com/home/webmap/viewer.html?url=http%3A%2F%2Fgeocatmin.ingemmet.gob.pe%2Farcgis%2Frest%2Fservices%2FSERV_GEOLOGIA_REGIONAL%2FMapServer&source=sd

        You can also open the attribute table to take a look at the data.

        There are only 15 polygons for that layer, however those polygons are very large so you have a massive amount of shape points per polygon. You would need to modify the python script to download each OBJECTID individually instead of a range of object IDs to see if that works for that particular map service. I just modified this section in the code:

        # Gather features
        print "Gathering records..."
        fs = dict()
        for i in idlist:
          print "  {}".format(i)
          urlstring = baseURL + "/query?objectIds={}&returnGeometry=true&outFields={}&f=json".format(i,fields)
          fs[i] = arcpy.FeatureSet()
          fs[i].load(urlstring)

        I’m running mine now. It is only on record 5 and has been running for about an hour! Good luck! -mike

  8. Hi Michael,

    I am trying to extract features from a Map Server using the python script. I am thinking it is not working due to I am using QGIS. I am not at all familiar with python and scripting. I was not able to find out how many feature are in the layers but it will not let me download them a json file. If you could help that would be great.

    The map server is in Georgian but that should not matter.
    Here is the link http://gisappsn.reestri.gov.ge/ArcGIS/rest/services/WebMap/MapServer

    I need to extract the Layer number 1 and 2. 1 is equal to Sectors and 2 is equal to quarters.

    Either shapefile or geojson would be fine.

    OH. the spatial reference of the data is 32638 and I would like the output to be 4326.

    Thanks so much if you can help.

    Dana

  9. Mike,

    This worked so well I had to leave a reply! Thank you so much, this is an amazing piece of code and between this post and your other “Extracting Features from Map Services” walk-through, its incredibly easy to implement even for novice coders such as myself. Thank you so much for contributing to the community and sharing this!
    -Kyle

      • Hi Mike, Happy New Year to you as well!

        I’ve successfully used the script multiple times now for different data sets and it has worked flawlessly with the following exception. There was one data set which produced the following error and I was wondering if there is something in the script that could be modified to fix the error:

        Python 2.7.14 (v2.7.14:84471935ed, Sep 16 2017, 20:19:30) [MSC v.1500 32 bit (Intel)] on win32
        Type “help”, “copyright”, “credits” or “license” for more information.
        >>> execfile(r’C:\MapService_Extract\script.py’)
        Record extract limit: 1000
        Number of target records: 11380
        Gathering records…
        OBJECTID >= 7372781 and OBJECTID <= 7373780
        Traceback (most recent call last):
        File "”, line 1, in
        File “C:\MapService_Extract\script.py”, line 42, in
        fs[i].load(urlstring)
        File “C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\arcobjects\arcobje
        cts.py”, line 175, in load
        return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args
        )))
        RuntimeError: RecordSetObject: Cannot load a table into a FeatureSet
        >>>

        The map server URL with the data set I was attempting to extract is here:
        http://www.cobbgis.org:81/arcgisext/rest/services/water/waterlivemapwm/MapServer/43

        It seems odd because many of the layers hosted by this same map server extract just fine but the storm water layers seem to be the ones returning the error.

        Thanks for your help!
        -Kyle

      • Hi Kyle. I have run into a few of these before. Looking at the high values of the record numbers … and the name of the map service … the data seems to be edited live which means old records that have been deleted or changed tend to stick around until they reconcile the data. That leaves you with a bunch of blank records not associated with any features, so there is nothing to save, thus the python script errors.

        The script would need to be modified to check if any features were selected first within the object id range (a feature count greater than zero), then if so, save the features and move on to the next group. If not, don’t save anything but move on to the next group. If you make any changes to the script to test for something selected, please share in this comments section! -mike

  10. Hi Miike. I tried just adding an exception to see if the script would go to the next record batch and gather any records that could be successfully loaded into the feature set, but it never seemed to find any records that could be loaded which makes me wonder if the issue really is blank records or if there is something else going on. I have very little understanding of arcpy (and my general python knowledge is limited as well) so please bare with me if my change to the code is not doing exactly what I think it’s doing. I am now using a different layer for testing which has far fewer total records:
    http://www.cobbgis.org:81/arcgisext/rest/services/water/waterlivemapwm/MapServer/39

    This was what I changed in the script:

    # 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 {} >> execfile(r’C:\MapService_Extract\script_long.py’)
    Record extract limit: 1000
    Number of target records: 11380
    Gathering records…
    OBJECTID >= 7384381 and OBJECTID = 7385381 and OBJECTID = 7386381 and OBJECTID = 7387381 and OBJECTID = 7388381 and OBJECTID = 7389381 and OBJECTID = 7390381 and OBJECTID = 7391381 and OBJECTID = 7392381 and OBJECTID = 7393381 and OBJECTID = 7394381 and OBJECTID = 7395381 and OBJECTID <= 7395760
    No records gathered, moving to next set!
    Saving features…
    Traceback (most recent call last):
    File "”, line 1, in
    File “C:\MapService_Extract\script_long.py”, line 52, in
    arcpy.Merge_management(fslist, outdata)
    File “C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\management.py”, li
    ne 4487, in Merge
    raise e
    arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
    ERROR 000400: Duplicate inputs are not allowed
    Failed to execute (Merge).

    Any thoughts would be appreciated! Thanks,
    -Kyle

    • Oops, I apologize, my comment seems to have formatted oddly, my change to the code was this:

      # 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()
      try:
      fs[i].load(urlstring)
      except RuntimeError:
      print "No records gathered, moving to next set!"

  11. Hi Michael,
    It is a great tool. I tried in 10.1 and works perfectly. However in 9.3 version it does not, Do you know why?.Shoud I change some commands?. Thanks a lot!

  12. Hi Michael,

    I tried running this script and it looked to be almost successful, but it failed towards the end with a SyntaxError: Invalid syntax, but I’m not sure where the invalid syntax is:

    (lots and lost of previous text)
    OBJECTID >= 50001 and OBJECTID >> # Save features
    … print “Saving features…”
    Saving features…
    >>> fslist = []
    >>> for key,value in fs.items():
    … fslist.append(value)
    … arcpy.Merge_management(fslist, outdata)
    File “”, line 3
    arcpy.Merge_management(fslist, outdata)
    ^
    SyntaxError: invalid syntax
    >>> print “Done!”

    Any suggestions?

    • Hi Chris. I see prompts “>>>”, are you entering each command by hand? Or are you running it inside the ArcGIS python window? Try running it in a windows command window instead and see what happens. It seems to error on the merge_management command.

      • Looking a little closer at your code, looks like you made the arcpy.Merge_managment command part of the for key,value loop. Take that out of the loop (remove the spaces before it) and the script should run without error. -mike

  13. Love this script, use it loads. Thank you!
    Is there anyway to run it against the entire Map Server? Ending the url in mapserver instead or mapserver/39? I’d like to be able to create a local gdb mimicking the map server without have to run the script 50 times. Any help is greatly appreciated! Thanks again!!

    • Hi Brita. You would need to program your python script to read in all the layer numbers for a particular map service, then iterate through the layers to capture the features in each. I do not have time to figure that out right now, but you can check out the API reference to see if you can. Just go to the map service URL that lists all the layers in the map service, then click on the API Reference link in the upper right, then click on “All Layers and Tables” to read about how to get a list. You would need to read in this list and grab the layer IDs which you would then use to iterate through each with the script. Good luck and do post your solution! -mike

  14. i was try this script, but i got error
    Traceback (most recent call last):
    File “D:\CBI\Master\download_data.py”, line 43, in
    fs[i].load(urlstring)
    File “C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py”, line 172, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))
    RuntimeError: RecordSetObject: Cannot open table for Load

    please help me

  15. i was try, but i got error

    Record extract limit: 1000
    Number of target records: 34
    Gathering records…
    FID >= 0 and FID <= 33
    Traceback (most recent call last):
    File "D:\CBI\Master\download_data.py", line 43, in
    fs[i].load(urlstring)
    File “C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py”, line 172, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))
    RuntimeError: RecordSetObject: Cannot open table for Load

    can help me to ressolve this?
    thanks

  16. i was try and succed, but didnt still end because i got error again

    Record extract limit: 1000
    Number of target records: 166240
    Gathering records…
    FID >= 0 and FID = 1000 and FID = 2000 and FID = 3000 and FID = 4000 and FID = 5000 and FID = 6000 and FID = 7000 and FID = 8000 and FID = 9000 and FID <= 9999
    Traceback (most recent call last):
    File "D:\CBI\Master\coba.py", line 42, in
    fs[i].load(urlstring)
    File “C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py”, line 172, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))
    RuntimeError: RecordSetObject: Cannot open table for Load

    can you help me, how to ressolve this? and If i just want to download 10.000 record data not all, what code i must to change?

    Thanks for your help

    • Nasmi, unfortunately you are dealing with a slow map service because the polygon data is massive and some polygons have thousands of shape points. Take a look at previous questions about how to tackle map services with large amounts of features. It shows how to modify your code. -mike

  17. Dear Mike

    When i try to access a small data
    http://geoportal.menlhk.go.id/arcgis/rest/services/test/Admin_Provinsi/MapServer/0

    i got error too

    Record extract limit: 1000
    Number of target records: 34
    Gathering records…
    FID >= 0 and FID <= 33
    Traceback (most recent call last):
    File "D:\CBI\Master\coba.py", line 42, in
    fs[i].load(urlstring)
    File “C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py”, line 172, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))
    RuntimeError: RecordSetObject: Cannot open table for Load

    can you help this error

    • Nasmi, this is another example of a slow map service and huge polygons that have thousands of shape points. Since there are only 34 records, in this case I would modify the maxrc variable to equal 1:

      maxrc = 1

      It should step through each of the 34 polygons to create your data. Other than that, I cannot help any further. -mike

  18. Hello, that is a great code and a great inspiration, but I want your help. When I run the code shows me this error:
    Record extract limit: 1000
    Number of target records: 3702326
    Gathering records…
    POG_CODE >= 1 and POG_CODE <= 1000
    Runtime error
    Traceback (most recent call last):
    File "”, line 42, in
    File “c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\arcobjects.py”, line 175, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))
    RuntimeError: RecordSetObject: Cannot load a table into a FeatureSet

    I tried lots of things but nothing could you please help me?
    Thank you again, have a nice day!

    • Hi Tzirineto. You really want to extract over 3.7 million records? Good luck! What map service are you trying to use? Others had the same issues as you. The answer to your issue will be in comments above yours. Just read them to find out what to do for your situation. -mike

      • You cannot exceed the max record count of the map service. If it is set to 1000, you can only extract 1000 or less records at a time. I suggest you change the maxrc variable in the Python script to something like 200, just change the code to maxrc = 200, or whatever number you want. Maybe extracting smaller groups of records will help, but will take a very long time in your situation since you are dealing with over 3.7 million records. I suggest you try contacting the organization that maintains the map service and see if you can get a copy of the data from them instead. -mike

  19. Hi Michael,

    Can you help me extract the table from the REST map service and save it as a tabular data in Geodatabases or table in ArcGIS using Arcpy?

    • Hi Eva. It can, but will take a long time since that map service only allows 1000 records at a time. Also everything is saved into memory before merged together and saved, so potentially you could use up your memory on your computer if you don’t have a whole bunch! The other option is instead of saving to memory, you modify the code to save to a geodatabase for each 1000 record query. In one of the comments I talk about how to do that. Good luck! -mike

      • Sorry for the late reply, but thank you for the response Michael. I was able to use this and it helped immensely! This was great. Thank you for providing this!

  20. Hi Michael,

    Thank you for your wonderful code, this was super helpful. I was wondering if you could help me, I am running into some issues with running the code for one of the layers I wanted to download. Here’s the map service link:

    https://arcgis.sd.gov/arcgis/rest/services/DENR/NR69_Observation_Wells_Public/MapServer/1

    The error that I’m getting is:
    Runtime error
    Traceback (most recent call last):
    File “”, line 23, in
    KeyError: ‘objectIdFieldName’

    The first part of the code I’m able to run and it returns the extract limit (2000). There’s over 5 million records on this so I changed the code to save to a gdb to save on memory usage.

    I saw someone had a similar error to mine in the comments, I changed the code to what you had suggested for downloading each object ID individually but still came with this error. I consulted with some friends and they suspect that this error might be due to not having unique Object IDs. I’m fairly new to python so forgive me if any of my questions seem elementary.

    Thanks,
    Annie

    • Hi Annie. So that layer in the map service has 5,506,335 points! When the code tries to gather all the object ids for each of those points, their server times out. There are just too many records for their server to return. You have two options. First, you could contact the owner of that server/data and ask them if they would cut you a copy of the data. If you don’t want to do that (or you know they would not share with you), then you have to be a little creative with the code. Let me know if you want to go the creative code route. -mike

      • Hi Mike,

        Thanks for your reply, I think I want to get a little creative with the code, I would have to go through several people and emails to ask for permission otherwise. I also would like to know more about how to solve this using code. Let me know if there’s any way I can help. No problem to drop this if this turns out to be a huge time sink for you, I’m grateful you’re taking the time to help me!

        Thanks again,
        Annie

      • OK, creative code it is! Since you will be dealing with over 5.5 million records, make sure to modify your code to save to a geodatabase instead of memory … see my post in the comments above on May 16, 2018.

        Now the fun part. We found that trying to get their server to give us a list of all record numbers was making the query timeout, so we need to find the first record number and last record number for ourselves. Visiting https://arcgis.sd.gov/arcgis/rest/services/DENR/NR69_Observation_Wells_Public/MapServer/1 we scroll down to the bottom of that page and there is a field named OBJECTID, which as an ArcGIS user you know that holds the unique record numbers for any GIS data. At the bottom of the page note there is a Supported Operations section. Click the Query link. This is actually the page that the script uses to query the data. To find the first record number, enter in the Where input field the expression OBJECTID = 5506402, Return IDs Only to True, and Return Count Only to False. Click the Query (GET) button again. Boom! Scroll down to the bottom and you will see the last record number. As of this writing it was 5528141 for me, but it might be different for you since they keep editing the data that number will change. Make a note of it.

        Now you need to modify your code. In the “Get object ids of features” section, comment out the code from “where =” to “idlist =” (that is put a # before each line) and change the line idlist = js[“objectIds”] to idlist = range(1, 5528141 + 1). Your last record number might be different, so put that in.

        That should do it. If you want to test, try changing the idlist variable to something like idlist = range(1, 4000 + 1). That would pull out the first 4000 records, should run really fast and you can see the results. If it looks good, change the idlist variable back to the full record count and sit back for a few days I bet! Should it blow up during the run, you can look at the geodatabase to see what set of records numbers were saved and start the range from there in your code to continue! Good luck! And tell us how long it took! -mike

      • Hi Mike,

        So sorry for my late reply, thank you so much for your suggested edits. I’ll let you know how it goes, thanks again for your help!

    • Hi Mike,

      Here’s my update on how the code went. So I did what you had advised to do, comment out parts of #Get Object Ids in the code and replace with the id list range. I got the error “NameError: name ‘idfield’ is not defined”, so I put back in idfield = js[“objectIdFieldName”]

      I got another error that said, KeyError: ‘objectIdFieldName’. I tried replacing with OBJECTID and that didn’t work, and I tried troubleshooting this for a while before finally giving up.

      I have another work-around that finally did work, and that was using Feature Class to Feature Class tool in ArcGIS Pro. After logging in ArcGIS Pro, I put in the url for my feature layer as the Input and ran the tool to output shapefiles (https://arcgis.sd.gov/arcgis/rest/services/DENR/NR69_Observation_Wells_Public/MapServer/1/). This worked and it took me approximately 1 hour to run the whole thing. Note that I had to try this out a couple of times, the tool kept failing and gave me the error “000210: Cannot create output ” and for some reason it worked one day, I’m not sure why to be honest. Anyways thanks for your incredible help again, the code may not have worked for this layer, but it helped me a lot with my other layers.

    • Hi Farin. The script works as designed and will handle map services with more than 3000 features for sure. Unfortunately, that map service you are trying to access is very slow/has gigantic polygon features and sometimes does not respond to queries, so the code will error out. You can look at the other comments posted here about modifying the code to deal with slow map services (for example, change the maxrc variable = 2), or better yet, you might want to contact whoever owns the map service and data and ask if they will give you a copy.

      To handle any errors that might occur during feature extraction, you can modify the code fs[i].load(urlstring) to the following:

      try:
        fs[i].load(urlstring)
      except:
        print "    ERROR: Could not extract records"

      At least that would catch the error and alert you, then continue on to try and fetch the rest. Good luck. -mike

  21. How do you modify the code so that the output it outputs spatial coodinates, such as WGS84?

    Also whenever I list GLOBALID as my fields, it will not output. I tried setting arcpy.env.preserveGlobalIds = True but it has no effect.

    • Hi Gerry.

      To output the geographic data in a different spatial reference, just add the &outSR argument to the query string for the urlstring variable. For WGS84, add &outSR=4326 to the query string.

      As for GLOBALID fields, I do not believe they can be output unless it is a Feature Service.

  22. Hello Milke,

    Thanks for the article. It has been really helpful. Is there a way to select parcels by location from a map server. I have a set of points where I would like to select parcels(from map server) within 500 feet of the point. The map server has about 5.5 Million parcels and maxCount of 5K, so I assume its a bit tedious to download them as a feature set locally and run Select By Location. Wondering if you’ve come across similar situation. Appreciate any ideas.

Leave a Reply to Eloy J. Carrillo Cancel 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 )

Connecting to %s