Extracting Features from Map Services

Posted on by

UPDATE 3/28/2018: See Extracting MORE Features from Map Services to get around the maximum record count limit.

ArcGIS Server map services are a pretty cool thing.  Someone else is providing GIS data for you to consume in a map or application.  And if they keep it up to date over time, even better.

However, there are some times when you just need the data instead so you can edit/manipulate it for your GIS analysis needs.  If you cannot acquire a shapefile or geodatabase, there is a way to extract the features from a map service.

For you to extract the features, there must be three things in place:

  • You must have access to the REST interface for the map service
  • The layer in the map service that you want must be a feature layer
  • The layer must have a query option with JSON format supported

To demonstrate this, let’s head on over to the map services provided by Southern California Association of Governments (SCAG).  Specifically, the one for city boundaries:

http://maps.scag.ca.gov/scaggis/rest/services/City_Boundaries/MapServer

Since we have access to the REST interface, we can check out the different settings of the map service.  Note there are two layers, City Boundaries (layer 0) and County Boundaries (layer 1).

extract1

We are interested in the City Boundaries layer, so click on it.  This brings up detailed information about the data.

extract2

They were good to include when the data was updated in the Description.  Also note that it is a feature layer and we can query it.  If you scroll down to the bottom you will see the query option in the Supported Operations section.

extract3

Also listed at the bottom are the fields associated with the features.  We should take a look at these next to see what fields we want.  To do this, go back to the previous page and then click on the ArcGIS.com Map option in the View In section at the top:

extract4

This will open up a default map on ArcGIS Online and display the map service layers.

extract5

Click on City Boundaries to view the two layers, then click on the icon that looks like a spreadsheet just below the City Boundaries layer.

extract6

This will bring up the attribute table for the city boundaries.

extract7

After taking a look at the data, all I want are the CITY, ACRES, and COUNTY fields.

Now let’s go back to the REST interface for the City Boundaries layer.  Take a look at the MaxRecordCount setting.  It is set to 1000.  This means that only 1000 records (features) can be returned for each query.  This can be extra work for you if this map service has more than 1000 features.  We need to find out how many features there are.  To do this, scroll to the bottom and click on the Query link in the Supported Operations section.  This will bring up the query interface for the map service.  This allows you to do a lot with the map service and get different results.  I will not go into the details of each parameter.  You can find out more info about them from the ArcGIS REST API reference.  For now, just enter “1=1” (without the quotes) in the Where input field.

extract8

This is a little trick to return as many records as we can.  Also make sure to set Return Geometry to False and Return Count only to True.  With Return Count only set to True, we will actually get all records (an exeption to the MaxRecordCount setting).  Once those are set, click the Query (GET) button at the bottom of the query form.  It looks like nothing happened, but scroll down to the bottom to see the result.

extract9

Note the count is 191.  So we don’t have to be concerned with the 1000 record limit.  We will be able to grab all the features at once.

NOTE: For those of you that are curious about the query, take a look at the URL in your web browser.  This is the REST query syntax used to get the record count.  We will be using something similar in a Python script to extract the actual geometry of all features.  You can mess around with the query form to see what happens and the URL that is created to make that happen.

We will be using a Python script to extract the features from the map service and create a featureset in a geodatabase.  Here is the script:

import arcpy

arcpy.env.overwriteOutput = True

baseURL = "http://maps.scag.ca.gov/scaggis/rest/services/City_Boundaries/MapServer/0/query"
where = "1=1"
fields = "CITY,ACRES,COUNTY"

query = "?where={}&outFields={}&returnGeometry=true&f=json".format(where, fields)

fsURL = baseURL + query

fs = arcpy.FeatureSet()
fs.load(fsURL)

arcpy.CopyFeatures_management(fs, "H:/scag/data.gdb/city_boundaries")

So what does this do?  First it loads the arcpy module and sets the overwrite option to True in case we want to run it again and overwrite the featureset.  A baseURL variable is set to the query URL of the city boundaries map service, a where variable for our where clause for the query, and a fields variable with a comma delimited list of the fields that we want.

Next a query variable is set with the actual query needed to extract the features.  Note the query is in quotes with {} for values that are plugged in by the variables listed in format() at the end.  Also note the returnGeometry parameter is set to true and f (for response format) is set to json.  What this query will do is extract all the records with just the fields we want, return all the geometry (the polygons), and output the whole thing in JSON format.  With the baseURL and query variables set, we merge them together for one long query URL as the fsURL variable.

Finally, we create an empty featureset (the fs variable), load the returned JSON data from the map service query into it, then copy it as a featureset to our geodatabase.

It takes a second or two to run the script.  Check out the result!

extract10

And here are the attributes we wanted for all 191 records.

extract11

I do want to point out that the data you get will be in the same projection as the map service.  You can find out what projection the map service is in by looking at the Spatial Reference setting in the REST interface.  Note this map service is using 26911:

extract12

What the heck is that?  Go over to spatialreference.org and enter the number in the search field.  It will tell you this is the NAD83 UTM zone 11N projection.  If you wanted something else, like NAD83 California State Plane zone 5 US feet (2229), you would edit the Python script to include the output spatial reference parameter in your query, like this: “&outSR=2229”.  The data is projected on the fly for you and saved to your geodatabase!

I hope this helps you out when you are looking for data in map services.  -mike

32 thoughts on “Extracting Features from Map Services”

  1. Great post!. Some food for thought…an even easier way is to use the new ArcGIS Open Data Sites. If you have an ArcGIS Online Organizational Account, just setup an Open Data site specifically for downloading data from REST endpoints…takes about 10 minutes to get it up and running with no coding (don’t get me wrong, I love coding). Basically you just need to create an ArcGIS Online listing for each item you want to be able to download and then share the item with your Open Data Site Group/Groups. After that you can launch your “data grabber” Open Data site and view all the listing..now you get a big “Download” button that will automatically make batch requests to the REST endpoint to download all the data. You can even use the Open Data Site interface to filter the data before you download it. I have used this technique multiple times and it works like a charm!

    Reply ↓

  4. Mike,

    Have you been able to modify the script to export more than 1000 features from a rest service url? I understand the query limit for most rest service urls is 1000, but I wonder if you may have come across a loop/iteration for this script since it was last posted. Thanks!

    Reply ↓

    • Hi Daniel. I have not tried changing the script to handle more than 1000 features. You would need to query the map service to return how many features there are, then if over 1000, iterate through by 1000 chunks to grab all of them and then put them all together. Hmm, I feel another project might be in my future! -mike

      Reply ↓

      • import arcpy
        # (17,030 fstopo quad features)
        arcpy.env.overwriteOutput = True

        for i in range(0, 18):
        print i
        fro = i * 1000
        to = i * 1000 + 999
        baseURL = “https://apps.fs.usda.gov/arcx/rest/services/wo_nfs_gstc/GSTC_FSTopoQuad_01/MapServer/0/query”
        #where = “1=1”
        print ‘OBJECTID >= ‘ + str(fro) + ‘ and OBJECTID = ‘ + str(fro) + ‘ and OBJECTID <=' + str(to)
        fields = "*"
        query = "?where={}&outFields={}&returnGeometry=true&f=json".format(where, fields)
        fsURL = baseURL + query
        fs = arcpy.FeatureSet()
        fs.load(fsURL)
        Groupoffeatures = i * 1000

        arcpy.CopyFeatures_management(fs, "C:/Users/Desktop/data.gdb/FSTopoGrids" + str(Groupoffeatures))

      • Thanks for the code! You commented out the where variable setting but did not set it for your query for your OBJECTID statement. Also I have found that you cannot assume the OBJECTID’s start at 0 and increment by 1 up to the total number of records. Sometimes I have seen some starting at 10,000 and ending at 100,000. I think that has something to do with services that allow access to data that is edited live, but I am just guessing about that. -mike

      • I was able to use the code successfully trying to get some data from parks and rec. However I came across some datasets that have more than 1000 features, one has 32,000. How do you grab 1000 piece chunks at a time? And is there an easier way to do that than repeating it 32 times? If I was to use dirtsniffers code, what lines would I change, and what would I keep the same (I’m new to coding so I have no clue what I’m doing).

  5. Pingback: Web Map Services – M.I.

  6. I was following your code it works fine, but when I am using a service that has a token in the path it throws an error. Where would I put the word “query”?

    This is the end of the url
    /PROJECTS/FeatureServer/0?token=pBLNMLQZYJclJ23noy5tZ9vNmLMJP9cLBSldYyYd-y_N9P6BDiO9G0Hl855V_3smxSWej4kiBUrUrOeS7FPuD37_2pjae0FqguJzKRiVrKzXw9l56LoYO4FQ2GnZKKmt0iDl6iNxCWQQYyfbhW_gHZtKvXJM3x0LKmM8DguNfD3BE5PgoBwxYisZP54B0JXN2Zd9oMG8TiUzK0T-4ra7ukrfAe4I6YB75ssZOxeco3Jb4AWkAGAKn-gWGCVvXfvF

    Reply ↓

    • Hi Rudy. I have never had to deal with secured services using python, so I cannot help you that much. I would think everything is the same in the query URL, except you just add the token parameter anywhere in the URL. In your example above, try adding your query parameter at the end of your token parameter. Like in my example I would have “… token=xxxx&where=xxx …”. See if that works. Also “query” would go after the zero in your URL, like “… FeatureServer/0/query …”.

      Reply ↓

    • Hi Lee. I will assume you are running Windows and have ArcGIS Desktop installed on it. Just copy the script and using a text editor save it on your computer with the .PY extension, like myscript.py . Make any changes in the script that is needed. Next you can just double click on the script to run OR you can open a DOS window and navigate to the file’s location and then at the prompt just type the file name to run it.

      Reply ↓

      • Hi Michael, I so appreciate your reply. And I am working on a unviersityproject that needs me to get the bathymetric contour data from a map service (https://gis.ohiodnr.gov/arcgis/rest/services/OCM_Services/LakeErie_Bathymetry/MapServer). And there are over 1000 record on that.. I tried the script and but I didn’t make it. Could you help me on that?

        And here is my script:
        import arcpy
        # (17,030 fstopo quad features)
        arcpy.env.overwriteOutput = True

        for i in range(0, 18):
        print i
        fro = i * 1000
        to = i * 1000 + 999
        baseURL = “https://gis.ohiodnr.gov/arcgis/rest/services/OCM_Services/LakeErie_Bathymetry/MapServer/1/query”
        #where = “1=1”
        print ‘OBJECTID >= ‘ + str(fro) + ‘ and OBJECTID = ‘ + str(fro) + ‘ and OBJECTID <=' + str(to)
        fields = "*"
        query = "?where={}&outFields={}&returnGeometry=true&f=json".format(where, fields)
        fsURL = baseURL + query
        fs = arcpy.FeatureSet()
        fs.load(fsURL)
        Groupoffeatures = i * 1000

        arcpy.CopyFeatures_management(fs, "C:\Users\Jie Li\Desktop\data\data.gdb" + str(Groupoffeatures)

      • Hi Lee. You have the idea in your script, but you need to save the data from each iteration, then merge it all together. Question, do you need the actual data for analysis, or do you need to just display the data in ArcMap? And maybe a call to Ohio DNR (http://geospatial.ohiodnr.gov/) to get a shapefile might be a better way? But then again, maybe your project requires programming to get the data? Let me know. -mike

  9. Hi everyone, I am actually working on a university project that needs me to get the bathymetric contour data from a map service (https://gis.ohiodnr.gov/arcgis/rest/services/OCM_Services/LakeErie_Bathymetry/MapServer). And there are 6928 records on that.. I tried the script (with loop) and but I didn’t make it. I really know little about coding so could anyone help me on that?

    And here is my script:
    import arcpy
    # (69,28 fstopo quad features)
    arcpy.env.overwriteOutput = True

    for i in range(0, 18):
    print i
    fro = i * 1000
    to = i * 1000 + 999
    baseURL = “https://gis.ohiodnr.gov/arcgis/rest/services/OCM_Services/LakeErie_Bathymetry/MapServer/1/query”
    #where = “1=1”
    print ‘OBJECTID >= ‘ + str(fro) + ‘ and OBJECTID = ‘ + str(fro) + ‘ and OBJECTID <=' + str(to)
    fields = "Z_VALUES"

    query = "?where={}&outFields={}&returnGeometry=true&f=json".format(where, fields)

    fsURL = baseURL + query

    fs = arcpy.FeatureSet()
    fs.load(fsURL)
    Groupoffeatures = i * 1000

    arcpy.CopyFeatures_management(fs, "C:\Users\Jie Li\Desktop\data\data.gdb" + str(Groupoffeatures))

    Reply ↓

  10. I am having a problem with an https://. Anyone able to help?

    At the fs = arcpy.FeatureSet() spot I receive an error

    RecordSetObject:Cannot open table for Load

    Reply ↓

      • Looks like it works for me. To test, I opened a Python command line window and entered the following at the prompts:

        >>> import arcpy
        >>> url = “https://anrmaps.vermont.gov/arcgis/rest/services/Open_Data/OPENDATA_ANR_WATER_SP_NOCACHE_v2/MapServer/186/query?where=1%3D1&returnGeometry=true&f=json”
        >>> fs = arcpy.FeatureSet()
        >>> fs.load(url)
        >>> arcpy.GetCount_management(fs)
        Result ‘3’
        >>>

        Just to check, I got a count of features, which was 3 and that is correct because there are only 3 in that map service. -mike

      • Interesting. I’m able to run it in the python command tool on demand. I am looking to put this code in an overnight geoprocessing python script that will be set off by the computers task scheduler. I am looking to grab an updated copy of the feature each night and update a separate point feature class that fall within these polygons. When I run the task to grab the feature through the task scheduler I receive the RecordSetObject:Cannot open table for Load error.

      • Try doing a Google search on “arcpy.FeatureSet() error scheduled task RecordSetObject Cannot open table for Load” and see what you get. Also, and I am sure you are doing this, run the scheduled task with a login that has administrative rights and access to ArcGIS.

      • I am not sure why but I changed https:// to http:// and I am able to run the script successfully in overnight geoprocessing and extract the features. Thanks for the time spent looking into this!

  11. Hi This is really awesome…i am trying to get some data from a municipality, I am able to view there files by adding the info via GIS server…but i’d like to make shapefiles as well. i am getting this error:
    RuntimeError: RecordSetObject: Cannot open table for Load

    i modified your script, but probably didnt get it right.

    import arcpy

    arcpy.env.overwriteOutput = True

    baseURL = “http://capeims.capecoral.net/arcgis/rest/services/Pipes/MapServer/0/query”
    where = “1=1”
    fields = “PIPE_SIZE,PIPE_TYPE”

    query = “?where={}&outFields={}&returnGeometry=true&f=json”.format(where, fields)

    fsURL = baseURL + query

    fs = arcpy.FeatureSet()
    fs.load(fsURL)

    arcpy.CopyFeatures_management(fs, “G:\Engineering\PROJ-ENG\C\CCST\04DD\GIS\CCST-Ex_SW.gdb\Pipes”)

    I am trying to get the pipe data from here: http://www.arcgis.com/home/webmap/viewer.html?url=http%3A%2F%2Fcapeims.capecoral.net%2Farcgis%2Frest%2Fservices%2FOpenData%2FPublicWorks%2FFeatureServer%2F15&source=sd

    Any recommendations? Thanks!

    Reply ↓

  12. Hi Mike, Thank you so much for writing this post! You just saved me so much time and heart break! Thanks and Happy Friday. -Josh

    Reply ↓

Leave a reply to Alex Dunko Cancel reply