Friday, 4 September 2015

Arcpy/REST: how to save raster dataset from URL


I have a script which downloads features from an ArcGIS MapServer via the REST API. Getting vector features is relatively trivial since the features can be retrieved as JSON. However, when retrieving a map as an image, I'm getting stuck at the point where I have retrieved a response.


The relevant portion of my code, which gives me a valid image URL, is this:


# Export the map to JSON.
urlGetImage = '%s/export?bbox=%s&bboxSR=%s&format=png&f=json' % (
theBasicUrl, bbox, bboxSR)
urlGetImageResponse = urllib2.urlopen(urlGetImage)
jsonImage = json.loads(urlGetImageResponse.read())


# Read the image URL from the JSON and save the image.
urllib.urlretrieve(jsonImage['href'], os.path.join(tempFolder, 'test.png'))

Note the parameter where I retrieve the response as JSON. This is not essential. I don't care what I retrieve it as, so long as I can retrieve the generated raster somehow.


A sample response is this:


{"width": 400, "href": "https://[SERVER NAME]:[PORT]/arcgis/rest/directories/arcgisoutput/imagery/DEM_MapServer/_ags_map26932bd6681a4546a951125dd7c60e88.png", "scale": 20326865.193372156, "extent": {"xmin": 2996004.971208999, "ymin": -3232131.5223958646, "ymax": -1080867.3202355732, "xmax": 5147269.17336929, "spatialReference": {"wkid": 102100, "latestWkid": 3857}}, "height": 400}

This gives me a path to download an image from, as well as the relevant spatial reference WKID. I then save it to an image.


My question is this: how do I go from here (an image on the drive) to a (referenced) raster dataset? Or have I already gone wrong in getting to this point?


EDIT:



To clarify, I need to work with a MapService, not an ImageService, as that's how the rasters will be published. When I export the result above as a KMZ (rather than JSON), and then load it in Google Earth, it appears in the correct location. So in theory, getting a raster image from the map server would seem possible. However, trying to convert this into a raster layer via arcpy.KMLToLayer_conversion gives:


  "ERROR 000401: No features were found for processing."

... and even if I extract the image information in the KML, download the image, and replace the path in the KML (which by default points to the generated image on the server), it still complains with this error.



Answer



It turns out there's a far more convenient, less manual, less error-prone way to do this. It was just a matter of getting the MapServer data as a KMZ and converting from KML to raster using all the parameters for the tool. For vector conversion, you only need two parameters for the KMLToLayer_conversion tool, but for raster / image conversion, you must specify all four. This is the crucial step that I missed earlier.


After some tweaking, and with some refactoring probably still required, here's a working approach (tested with non-tiled map):


## imports: arcpy, os, json, glob, urllib

def getMapImage(theBasicUrl, theClipJsonString, theDestinationFolder,

theDestinationGDB):
"""
Retrieve a map image via URL based on an extent.

Parameters:
theBasicUrl: The URL to which queries can be attached.
theDestinationFolder:
The directory containing the destination GDB.
theDestinationGDB: The ultimate output location.
theClipJsonString: Area to be clipped (as JSON).


Returns:
outputFeature: The location of the output feature as a string.
"""

# Create a temp folder to work in.
# Note: 'now' is a global variable containing a datetime string.
tempFolder = os.path.join(theDestinationFolder,'temp_%s' % now)
os.makedirs(tempFolder)


# Set the image format (could also be a parameter).
imageFormat = 'png'

# Create a JSON file out of the clip area JSON
# (required for JSONToFeatures_conversion).
jsonFile = os.path.join(tempFolder, 'clip.json')
with open(jsonFile, 'w') as outFile:
outFile.write(theClipJsonString)

# Turn the clip shape into a feature class.

clipFeatureClass = os.path.join(tempFolder, 'clipFC.shp')
arcpy.JSONToFeatures_conversion(jsonFile, clipFeatureClass)

# Get the clip feature class's extents.
clipExtents = arcpy.Describe(clipFeatureClass).extent

# Construct a string for the request URL's bbox parameter.
bbox = '%s,%s,%s,%s' % (
clipExtents.XMin,
clipExtents.YMin,

clipExtents.XMax,
clipExtents.YMax)
bboxSR = json.loads(theClipJsonString)['spatialReference']['wkid']

# Download as a KMZ.
urlGetImage = '%s/export?bbox=%s&bboxSR=%s&format=%s&f=kmz' % (
theBasicUrl, bbox, bboxSR, imageFormat)
kmzFile = os.path.join(tempFolder, 'test.kmz')
urllib.URLopener().retrieve(urlGetImage, kmzFile)


'''
NOTE: a caveat applies here which concerns the possible downloading of
source rasters at all available scales, as mentioned in the docs for
KMLToLayer_conversion.
'''

# Convert the KMZ to a raster.
rasterizedKMZ = 'mapLayer'
arcpy.KMLToLayer_conversion(kmzFile, tempFolder, rasterizedKMZ, True)


# Place the raster in the GDB.
rasterizedKMZFile = glob.glob(os.path.join(tempFolder, '%s.grd' % \
rasterizedKMZ, '*.%s' % imageFormat))[0]

outRasterName = 'map_export_%s' % now

arcpy.CopyRaster_management(
rasterizedKMZFile, os.path.join(theDestinationGDB, outRasterName),
'#', '#', '256', 'NONE', 'NONE', '#', 'NONE', 'NONE')


# Delete intermediates
arcpy.Delete_management(tempFolder) # this doesn't work yet

return os.path.join(theDestinationGDB, outRasterName)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...