2

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.

2

Are you working with an ImageService? From what I can tell is that you are working with a map service and if that is the case, the map export option will only return a static image of the map service or individual layer.

Looking at what gets returned in your example, we can see there is a png being returned. This is not a true raster but rather just an image of it. Therefore, it will not be georectified and it will not have the pixel values of the raster (they will just be RGB values).

These appear to be services on your own servers? If so, you can just publish the raster by itself as an image service and then you can get the image back as a tif, which will have the spatial reference information as well as the pixel values.

EDIT

After further discussion, we determined that you just needed to just georeference the output map image. I have written a small function based on your code above to do this and it worked for me:

import urllib, urllib2, json, os, arcpy
arcpy.env.overwriteOutput = True

def exportMap(url, bbox, bboxSR, out_png):

    # get bounding box properties (assumes is comma delimited)
    xmin, ymin, xmax, ymax = map(float, bbox.split(','))
    size = ','.join(map(str, [int(abs(xmin-xmax)), int(abs(ymin-ymax))]))

    # make requests         
    urlGetImage = '%s/export?bbox=%s&bboxSR=%s&size=%s&format=png&f=json' %(url, bbox, bboxSR, size)
    urlGetImageResponse = urllib2.urlopen(urlGetImage)
    jsonImage = json.loads(urlGetImageResponse.read())

    # Read the image URL from the JSON and save the image.
    tmp_png = os.path.join(arcpy.env.scratchFolder, 'temp_png.png')
    urllib.urlretrieve(jsonImage['href'], tmp_png)

    # read the raster extent with arcpy and determine shift params
    ext = arcpy.Describe(tmp_png).extent
    ext_json = jsonImage['extent']
    x_shift = ext_json['xmin'] + ext.XMin
    y_shift = ext_json['ymax'] + ext.YMax #might have to adjust these?

    # run shift tool
    arcpy.env.outputCoordinateSystem = ext_json['spatialReference']['wkid']
    arcpy.management.Shift(tmp_png, out_png, x_shift, y_shift)
    arcpy.management.Delete(tmp_png)
    return out_png

if __name__ == '__main__':

    # test it
    url = 'http://arcgis.dnr.state.mn.us/arcgis/rest/services/elevation/mn_dem_1m/MapServer'
    bbox = '230720.3137,4954797.5862,232077.6289,4955591.3378'
    sr = 26915
    out_png = r'C:\TEMP\Map_export.png'
    exportMap(url, bbox, sr, out_png)
    print 'done'

note: if the map service is tiled, it will not align exactly because it is pulling out the actual tiles for the output image.

  • Thanks, but I am constrained to working with a MapService. I will modify the question to reflect this. – Wolfie Inu Oct 14 '15 at 12:29
  • What kind of raster data set are you extracting? – crmackey Oct 14 '15 at 12:31
  • At the moment it's a DEM being served as a MapService, but it could in theory be any map that the client decides to make. As you mention in your answer, the value of such a dataset (even when successfully extracted) is questionable. But I'm given to understand that the client's connection is too slow to just have them use a WMS directly in ArcMap without causing frustration, so they'd prefer to extract an image of it to work with, clipped to a polygon that they provide as a parameter (that I'm approximating with a bounding box). Odd little use-case. – Wolfie Inu Oct 14 '15 at 12:41
  • Gotcha. I struggled with this as well at one point and ended up giving up because I couldn't quite get the image to reference correctly (see post here). Also, an ImageService is different from a WMS service. One thing you can do to is to publish the DEM as an image service and give the client an Add-In that allows them to draw a box then clip out the DEM from your image service. This works very quickly (we use a 1m DEM image service for the entire state of MN and a tiff is returned in ~6 seconds) – crmackey Oct 14 '15 at 13:12
  • 1
    Thanks for all the effort - turns out this didn't work for me due to some CRS-related shenanigans. I ended up going with the KMZ approach after all - please see the accepted answer – Wolfie Inu Oct 16 '15 at 11:12
1

The downloaded image does not contain any information about it's geographical location (at least not to my knowledge), so you need to create a world file manually. The format is trivial, see the ArcGIS Desktop help on World files for raster datasets.

Additionally, you need to set the spatial reference for the image, use arcpy.DefineProjection_management() to do that: http://desktop.arcgis.com/en/desktop/latest/tools/data-management-toolbox/define-projection.htm.

EDIT

I just thought of an alternative: You could create an (un)managed raster catalog, the necessary geographic information then goes into a table. Please consult the help for the Create Raster Catalog toolbox.

  • How would I create a world file (or deduce the necessary parameters to create a blank raster corresponding to the relevant area) without an input raster dataset? And how would I define the projection for a plain image file (which is not a feature layer or geodataset)? – Wolfie Inu Oct 14 '15 at 12:27
  • A world file is simply a text file containing 6 numbers defining the extent of the image, see the help for ArcGIS Desktop. You should be able to create that with Python. Name it test.pngw. Since you downloaded the image you already have a dataset use it's filename in the DefineProjection command. – Berend Oct 14 '15 at 12:43
  • So when the reference for the DefineProjection_management tool says that its input needs to be either "Feature Layer" or "Geodataset," that includes a plain image file? – Wolfie Inu Oct 14 '15 at 13:22
  • @WolfieInu Exactly, an image is a geodataset as far as ArcGIS is concerned. In ArcCatalog, right-click on an image file, and open it's properties window. It will show Data Type: File System Raster – Berend Oct 14 '15 at 14:27
  • 1
    Since the resulting JSON returns the extent info, I think it is easier to just run the Shift tool to put it in the right spot rather than generating a world file. See the edit to my answer, I ran the code on several different map services and the output image is correctly georeferenced (well, it is off by ~1m, may need to play with the math a little bit). – crmackey Oct 15 '15 at 14:55
1

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)
  • 1
    Very cool! Glad you found a working solution. Yeah it's probably safer using this since the coordinate info is automatically tied to the KML. – crmackey Oct 16 '15 at 14:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.