Download Queried Raster from ArcGIS REST

Question

I am trying to download/export a raster based on a query of USGS's 3DEP Image Service. The endpoint (if that's the correct terminology) is here:

https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage?bbox=-2.00375070672E7,-1689391.823360186,2.0037507532527123E7,1.88090019719E7

Here is what this looks like using the HTML interface. Clicking "Export Image (GET)" initiates a download of the raster as a .tif file, exactly as I need.

I would like to be able to do this using Python/ArcPy exclusively (updating the parameters as needed), as part of a larger script. I have had a shockingly difficult time Googling this or finding any resources about how to do this, so I'm missing some very basic concepts here.

The code below returns... something, but I don't know how to interact with whatever it is. I just need the actual raster data (a tif? a numpy array? anything?) to do some other geoprocessing with.

import requests

url = \
    r"https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage?"

params = {"f": "html",
          "bbox": [-2.00375070672E7, -1689391.823360186, 2.0037507532527123E7, 1.88090019719E7],
          "imageSR": 102005,
          "format": "tiff",
          "pixelType": "S16",
          "noDataInterpretation": "esriNoDataMatchAny",
          "interpolation": "RSP_BilinearInterpolation"}

response = requests.get(url=url, params=params, allow_redirects=True)
print(response.text)

# Something like this...
# queried_DigitalElevationModel = response.getData()?
# arcpy.Raster(queried_DigitalElevationModel).save("localfilepath or memory")

VinceE · Answer

Thanks for the assistance @Anonymous User, I was able to get this working. I have a follow up question that I'm hoping you may have some insight into. My modified code is below. The extent string is computed beforehand, but seems to work as expected, based on the polyline geometry being read via a SearchCursor. def save_raster(name, output_directory, ext_string): """""" url = r"https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage?" params = {"f": "image", "bbox": ext_string, "bboxSR": 102005, "imageSR": 102005, "format": "tiff", "pixelType": "F32", "size": "600,600", "noDataInterpretation": "esriNoDataMatchAny", "interpolation": "RSP_BilinearInterpolation"} response = requests.get(url=url, params=params, allow_redirects=True) name = "".join([char if char in ascii_letters else "" for char in name]) fullpath = pathlib.Path(output_directory, f"{name}.tif") with open(fullpath, mode="wb") as image_file: image_file.write(response.content) arcpy.management.CalculateStatistics(str(fullpath)) The question I have is about specifying pixel size, which oddly enough does not seem to be an option based on the documentation here: https://developers.arcgis.com/rest/services-reference/enterprise/export-image.htm Below my pixels are roughly 10x10 (meters) because I've specified a size of 600,600. Here my pixels are much, much larger because the size I've specified is 30x30, width/height. Using the dimensions of the polyline extent (extent.XMax - extent.XMin) and the intended pixel resulution (30 meters, for example), some math can be performed to get the "size=," that yields approximately the intended pixel resolution, but it's pretty cumbersome compared to just passing the resolution as a parameter. It's very odd to me that the call takes a pixel count (when is this property relevant?) instead of the actual pixel resolution, which is one of the more important properties of a raster dataset...

Anonymous User · Answer

You're close and a couple things for it-  bbox should just be a string, and the return format (f) should be image, and lastly, you'll need to tell it to save it (response.content) to disk.

import requests

url = r"https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage?"

params = {"f": "image",
          "bbox": '-2.00375070672E7,-1689391.823360186,2.0037507532527123E7,1.88090019719E7',
          "imageSR": 102005,
          "format": "tiff",
          "pixelType": "S16",
          "noDataInterpretation": "esriNoDataMatchAny",
          "interpolation": "RSP_BilinearInterpolation"}

response = requests.get(url=url, params=params, allow_redirects=True)

with open(r'C:\...\Documents\ArcGIS\outtest_s16.tif', mode='wb') as localfile:
    localfile.write(response.content)

gives me this: