Access layer of raster stored in memory using ArcPy to create NDVI from subset of multi-band raster

Question

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'])

Answer 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.).