Extract a Map Service Layer to Shapefile using Python

How easy/hard is it to extract a Map Service layer to Shapefile using Python? Not as hard as you would think. Is it legal? Not as easy as you would think! And before I get into this any further, there are also ways of doing this directly in ArcMap - but that wouldn’t be very Pythonic now would it.

UPDATE: Part 2 is now available - Map Service to Shapefile with Python Part 2 – Iteration

Legal Schmegal (and a brief disclaimer)

I’m not a lawyer - but I have watched many TV shows where someone else played one. Oh, and I’ve actually read the Canadian Copyright Act (some of it sunk in), done a few presentations on the topic, and deal with the subject at work on a daily basis. Why do I bring this up? Just because we can do something, doesn’t mean we should. Please make sure you have permission from the owner of the Data/Map Service, or have permission via a License (“Licence” in Canada) before running around downloading shapefiles from the internet.

Oh, and if you are worried about someone downloading from your service: Best practice would be to secure the service, or at least limit the service functionality. Possibly adding a note into the service description on intended use could also help. But if you really want to keep your data secure, please refer to the helpful venn diagram below.

Venn Diagram by Rob Jewitt

Map Service to Shapefile using Python

Alright, now that the legal stuff is out of the way - let’s get down to some coding! We are essentially just going to wrap 3 small steps into a python file:

Query the web service layer for data
Save the results locally
Convert results to a shapefile (or other feature class type)

1 Query the web service layer

For this example we will download US States from Esri’s ArcGIS Server sample server. Once we have the service and layer, we can investigate the attributes and query options using the REST interface or just jump into Python. Either way, we are just going to call the same REST query inside Python using the urllib2 library. The trick here is to make sure the result format type is set to JSON.

myRequest = <URL>
response = urllib2.urlopen(myRequest)
myJSON = response.read()

2 Save the Results

Now that we have the data, we can save it locally as a text file with the JSON extension. Nothing fancy here just create a file object and save the results.

foo = open("jsonOutput.json", "wb")
foo.write(myJSON);
foo.close()

3 Convert the result to a shapefile

With the data successfully downloaded in a JSON file, all we need to do now is convert it to a shapefile. This task is made pretty simple seeing as there is a tool for that, JSON to features.

arcpy.JSONToFeatures_conversion("jsonOutput.json", "finalShapfile.shp")

What could possibly go wrong?!

This approach assumes the service you are connecting to allows query operations, otherwise you will get an error message immediately. Try the query operation directly from your browser first, then copy the parameters/syntax into Python for a quick sanity check.

Another thing with queries is that results are limited. There is a setting in ArcGIS Server to limit the number of features returned from a query operation (1000 by default). You may have noticed in the sample code, the _where _clause was set to where=1%3D1 (the %3D is web speak for equals). This is a simple true query that will return all features to the max results limit.

If more features exist or if you’re not sure how many features there are in the first place, just query it. Make a query and set the Result Option returnCountOnly=true. This returns the count so you can determine if one query will be enough. For the example with US States, there are only 51 features, so the 1=1 clause isn’t a problem. If the number was above the return limit, you would need to use an iterative loop until all features are returned. Don’t worry, guessing ObjectIDs starting at 1 isn’t required since you can send a returnIdsOnly=true query to get the actual list of ObjectIDs - and this query doesn’t have the max return limit applied. Just use the IDs to divide into chunks less than the limit, and loop away.

A final issue that can occur is that the output shapefile isn’t a shapefile - just a DBF. This occurs when a layer in a map service is published but the SHAPE field is hidden from the service. In this case, you still get all the attributes but no spatial component as the SHAPE field must be visible in the service for the geometry to be returned.

Map Service layer to Shapefile - Working Code

Here is the full python script. The URL is split into 3 parts only to help show the various components - a single URL parameter could be used instead. Please note: ArcPy is required for the JSON to Feature tool.