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

Question

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"\\xyz\xyz\xyz\xyz\xyz\xyz\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 asterCalResult_{}".format(now)) arcpy.Delete_management(in_data="memory") ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ scriptTool.py.zip

MarieCline_Delgado · Answer

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?

DanPatterson_Retired · Answer

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 · Answer

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