1

I am trying to create a final output of an NDVI raster for a number of Landsat images run through a for loop in Python 3 (for ArcGIS Pro). I am able to produce the NDVI rasters when I save the intermediary Landsat subset from ExtractByMask, but receive the following error when I instead try to store the intermediary subset in memory:

nir_raster = arcpy.Raster(nir)

RuntimeError: ERROR 000732: Input Raster: Dataset in_memory\myMemoryRaster\Layer_4 does not exist or is not supported

Here is the code where I attempt to refer to layers of the raster in memory (in my case layers 4 and 3 refer to NIR and red bands, respectively in my particular multi-band raster).

import arcpy, os

arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
arcpy.env.workspace = 'in_memory'

# Define area of interest for subset
aoi = 'P:/aoi/aoi.shp'
aoi_name = os.path.basename(aoi)

# Select Landsat Images 
images= "P:/Landsat/L7_p39r37_01242016_b2-5_7.img; P:/Landsat/L8_p39r37_02172016_b2-7.img"
image_list = images.split(";") # used here because real script takes multi-value parameter for images

# Select output folder to which NDVI rasters will be saved to
output_workspace = 'P:/NDVI'

# 1. Generate NDVI raster
for i in image_list:

    image_name = os.path.basename(i) 
    image_name_chunks = image_name.split('_')
    ndvi_name = 'ndvi_' + image_name_chunks[2] +'.img'
    ndvi = os.path.join(output_workspace, ndvi_name)

    # 1.1 Subset image
    outExtractByMask = arcpy.sa.ExtractByMask(in_raster = i, in_mask_data = aoi)
    output = "in_memory" + "\\" + "myMemoryRaster"
    outExtractByMask.save(output)

    # 1.2 Calculate NDVI on subset
    nir = output + '\Layer_4'
    nir_raster = arcpy.Raster(nir)

    red = output + '\Layer_3'
    red_raster = arcpy.Raster(red)

    numerator = arcpy.sa.Float(nir_raster - red_raster)
    denominator = arcpy.sa.Float(nir_raster + red_raster)

    ndvi_output = arcpy.sa.Divide(numerator, denominator)
    ndvi_output.save(ndvi)

arcpy.Delete_management('in_memory')

However, when I substitute the following code in sections 1.1 and 1.2 (i.e. not saving subset to memory) my code runs fine:

# 1.1 Subset image
subset_name = image_name[:-4] + '_' + aoi_name[:-8] + '_subset.img'

subset = os.path.join(output_workspace, subset_name)

subset_output = arcpy.sa.ExtractByMask(in_raster = i, in_mask_data = aoi)
subset_output.save(subset) 

# 1.2 Calculate NDVI on subset

nir = subset + '\Layer_4'
nir_raster = arcpy.Raster(nir)
red = subset + '\Layer_3'
red_raster = arcpy.Raster(red)

I've visited: Access individual bands and use them in map algebra and What is the proper syntax and usage for ArcGIS "in_memory" workspace?

Update: adding multi-line version of user2856's Describe note below:

for child in arcpy.da.Describe(outExtractByMask)['children']:
    print(child['catalogPath'])
  • 1
    Could you please show what the following output looks like: print(nir)? – Aaron May 24 at 4:02
  • @Aaron For sure. It reads: in_memory\myMemoryRaster\Layer_4 – Iracambi May 24 at 16:08
1

You may have the band name wrong, check them using arcpy.da.Describe, i.e.

print([child['catalogPath'] for child in arcpy.da.Describe(outExtractByMask)['children']])

I get the band names as 'Band_1' to 'Band_<N>'.

You could create a little function* to automatically pull out the band names, then your code would look something like:

import arcpy, os

arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
arcpy.env.workspace = 'in_memory'

def get_bands(raster):
    bands = [arcpy.Raster(band['catalogPath'] for band in arcpy.da.Describe(raster)['children']]
    return bands

# Define area of interest for subset
etc......

# 1. Generate NDVI raster
for i in image_list:
    etc...

    # 1.2 Calculate NDVI on subset
    bands = get_bands(outExtractByMask)
    nir_raster = bands[3]  # list begins at 0, so bands[3] = Band 4
    red_raster = bands[2]

* Note for any ArcGIS Desktop users who may try to use this function, arcpy.da.Describe is only available in ArcGIS Pro arcpy. This code will not work in ArcGIS Desktop (ArcMap etc.).

  • Very interesting! print([child['catalogPath'] for child in arcpy.da.Describe(outExtractByMask)['children']]) reads: ['in_memory\\myMemoryRaster\\Band_1', 'in_memory\\myMemoryRaster\\Band_2', 'in_memory\\myMemoryRaster\\Band_3', 'in_memory\\myMemoryRaster\\Band_4', 'in_memory\\myMemoryRaster\\Band_5', 'in_memory\\myMemoryRaster\\Band_6'] – Iracambi May 24 at 16:09
  • Thank you! It ran perfectly once I simply changed Layer to Band in section 1.2 as in: nir = output + '\\Band_4' – Iracambi May 24 at 16:14
  • Though it's confusing to me why the result of print(nir) still reads in_memory\myMemoryRaster\Layer_4? Band is nowhere there. – Iracambi May 24 at 16:18
  • @Iracambi I believe you are explicitly defining Layer_4 in your script. – Aaron May 24 at 17:13
  • @Aaron, Ah, yes, of course! Thank you :) – Iracambi May 24 at 17:19

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.