Trying to Geoprocess with a Single Variable from a Multidimensional Mosaic Layer

[MarieCline_Delgado]

I have a multidimensional mosaic layer which contains 3 raster datasets.  I am writing a tool that reads one of the raster datasets into a numpy array by applying a definition query on the Variable associated with the raster.  However, I can't seem to find a way to successfully accomplish this. 

Regardless of the definition query I apply, the output ALWAYS uses the same Variable.  
Example:  I have Variables SnowA, SnowB, SnowC
If I set the definition query to Variable = SnowA, the geoprocessing results will reflect those of SnowC

If I set the definition query to Variable = SnowB, the geoprocessing results will reflect those of SnowC

If I set the definition query to Variable = SnowC, the geoprocessing results will reflect those of SnowC

Any help would be greatly appreciated!

import arcpy, datetime, re
now = datetime.datetime.now().strftime("%Y%m%d_%H%M")

# Define input Variables
ptFLS=arcpy.GetParameter(0) # Feature layers 
mosLyr=arcpy.GetParameterAsText(1) # Mosaic layer
minThresh=int(arcpy.GetParameterAsText(2)) # Integer representing inches
outLoc=arcpy.GetParameterAsText(3) # Out Folder
arcpy.env.overwriteOutput = True

# Other variables
aprx = arcpy.mp.ArcGISProject("CURRENT")
symbologyLyr = r"\\\\XXXX\\XXXX\\XXXX\\XXXX\\XXXX\\XXXX\\HighlightedPoints_toolDependency.lyrx"
saveName = str(arcpy.Describe(mosLyr).nameString).replace("\\", "")
groupName = arcpy.Describe(mosLyr).nameString+"\Image"
mosWhere = "{}".format([lyr.definitionQuery for lyr in aprx.activeMap.listLayers() if lyr.supports("LONGNAME") and lyr.longName == groupName][0])
procTemp = re.findall(r"'([^']*)'", mosWhere)[0]

# Get raster information
arcpy.AddMessage(" > Reading raster properties...")
arcpy.MakeMosaicLayer_management (in_mosaic_dataset=mosLyr, out_mosaic_layer="mos_Lyr", where_clause=mosWhere, processing_template=procTemp)
rast = arcpy.Raster("mos_Lyr")
desc = arcpy.Describe(rast)

# Get upper left corner of the raster and spatial reference
ulx = desc.Extent.XMin
uly = desc.Extent.YMax
sr = desc.spatialReference

# Convert raster to numpy array
arcpy.AddMessage(" > Creating raster array...")
mosLyr_npArr = arcpy.RasterToNumPyArray (in_raster=rast)
xMax, yMax = mosLyr_npArr.shape

# Calculate row and col numpy array position based on how far the point is from the upper left corner and pass result into numpy array to extract
for ptFL in ptFLS:
    arcpy.AddMessage(" > Reading {} data...".format(ptFL))
    arcpy.FeatureClassToFeatureClass_conversion(in_features=ptFL, out_path="memory", out_name="ptFC")
    arcpy.AddField_management(in_table=r"memory\\ptFC", field_name="inUnits", field_type="FLOAT", field_scale=1)
    arcpy.AddMessage(" >> Determining impact at client locations...")
    with arcpy.da.UpdateCursor(r"memory\\ptFC", ["SHAPE@","inUnits"]) as cur:
        for row in cur:
            pnt = row[0].projectAs(sr)
            # Assumes every point falls to the right and below the uperleft corner ... some points fall beyond footprint of raster i.e. Hawaii, Carribean, eronious misplaced points
            deltaX = pnt.centroid.X - ulx
            deltaY = uly - pnt.centroid.Y
            arow = round(deltaY/rast.meanCellHeight)
            acol = round(deltaX/rast.meanCellWidth)
            if arow < xMax and acol < yMax:
                row[1] = mosLyr_npArr[arow,acol]
            cur.updateRow(row)

    arcpy.AddMessage(" >> Exporting impacted locations to output folder...")
    arcpy.MakeTableView_management (in_table=r"memory\\ptFC", out_view=r"memory\\impactLocTBL", where_clause="{} > {}".format("inUnits",  minThresh))
    arcpy.TableToExcel_conversion (Input_Table=r"memory\\impactLocTBL", Output_Excel_File=outLoc+"\SnowStormReport_{}_{}_{}.xls".format(saveName,ptFL,now))
    arcpy.MakeFeatureLayer_management(in_features=r"memory\\ptFC", out_layer="Impacted{}_Layer".format(ptFL), where_clause="{} > {}".format("inUnits",  minThresh))
    #arcpy.ApplySymbologyFromLayer_management (in_layer="Impacted{}_Layer".format(ptFL), in_symbology_layer=symbologyLyr)
    arcpy.SaveToLayerFile_management(in_layer="Impacted{}_Layer".format(ptFL), out_layer=outLoc+"\SnowStormReport_{}_{}_{}.lyrx".format(saveName,ptFL,now), version="CURRENT")
    arcpy.CopyFeatures_management("Impacted{}_Layer".format(ptFL), r"C:\Users\mde\Desktop\Projects\scratchSpace.gdb\rasterCalResult_{}".format(now))
    arcpy.Delete_management(in_data="memory")

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

What does the script look like, and the numpy array?

[MarieCline_Delgado]

Hi Dan - I updated the original post with the script.  The numpy array shape is [850, 1500].  Which means it's only reading in a single dimension, right?  It's just reading in the wrong dimension...

[DanPatterson_Retired]

not sure what the script is doing... but

you have your array's UL corner... you have the rows and columns, you need the cell size and/or the raster extent to determine the metric unit that one cell represents. You don't need to do that for every point in the raster.

you want to query the array I presume and/or 'do' something with it.  Queries are easy

a = np.random.randint(0, 10, size=(10,10))  # ---- make some fake array data

np.where(a == 5, 99, a) # ---- if 'a' == 5 assign it 99

array([[ 6,  2,  9,  3,  3,  7,  7,  4,  0,  1],
       [ 8,  1,  8,  3,  8,  7,  6,  1, 99,  6],
       [ 6,  2, 99,  6,  9,  4,  7,  1,  9,  4],
       [99,  3,  4,  9,  0,  6,  4, 99,  7,  1],
       [ 8, 99, 99,  6,  3,  0,  7, 99,  8,  8],
       [ 1, 99, 99,  0,  2,  0,  4,  2,  0, 99],
       [ 1,  9,  7,  2,  8,  8,  9,  3, 99,  0],
       [ 0,  3,  1,  2,  8,  7,  4,  0, 99,  2],
       [ 4,  9,  4,  7, 99,  6,  9,  3,  4,  4],
       [ 4,  1,  4,  9,  2,  9,  2,  4,  9,  9]])

# ---- an indicator raster, similar question, different result

b = (a == 5).astype('int32')

b
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
       [1, 0, 0, 0, 0, 0, 0, 1, 0, 0],
       [0, 1, 1, 0, 0, 0, 0, 1, 0, 0],
       [0, 1, 1, 0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

So the query language and outputs are very much like the spatial analyst, but the terminology and syntax is different... just like any language

[MarieCline_Delgado]

The tool is just suppose to mimic the Extract Values to Points tool (I do not have Spatial Analyst).  

This tool inputs include a point feature class and a mosaic layer.  I want the tool to use a single, user-specified dimension from my multidimensional Mosaic Layer and put the raster values in an array (using RasterToNumPyArray).  The Mosaic Layer is in a map and has a multidimensional filter applied as well as a definition query to display the desired data.  The problem is that when I use this mosaic layer as the input, the methods do not seem to recognize the filter and query on the mosaic layer.  Regardless of which dimension is filtered, the values for the same dimension (call it "SnowC") always wind up being in the array.  

I have tried being more forceful and using arcpy.MakeMosaicLayer_management (line20) with the input mosaic layer to apply the filter and query properties; however, the same result still occured - the raster values were for the SnowC dimension and not the intended SnowA dimension.

Any thoughts on what I may not be doing correctly or a method that could read the desired mosaic layer dimension into the array?