Creating index layer around rasters to rename files using ArcPy with ArcMap?

I have a few hundred high resolution geotiffs (cell size 2m) such as the example below which I am trying to rename according to the location in an index layer (red rectangles in example below). I would like to add the attribute of the index rectangle in front of the filename (in this case "B6_filename.tif") which covers the most of the data pixels in the geotiff (black are NoData values, coloured values are data pixels).

What I tried so far is this:

• Reclassify the raster to a binary mask with data pixels and NoData values using the reclassify tool


• Convert Raster to Polygon (takes over an hour probably due to the high resolution data and is therefore not suitable when trying to automate this for hundreds of geotiffs with arcpy)


• Creating an extent shapefile using the Raster Domain tool available within the 3d Analyst licence

I was thinking of further using this extent shapefile together with an Intersect analysis to calculate the area percentage of each intersected part. However the extent of the raster extracted by the Raster Domain tool results in a rectangle and not around the data pixel values (blue rectangle in example below).

My approaches were not successful.

Is there another way to solve this problem?

In the end I am attempting to rename my sample file "filename.tif" to "B6_filename.tif" with B6 referring to the proportionally largest area covering the index rectangle.

Answer

This is an approach using raster extent, in case you dont get an answer using the actual pixels. It runs very quickly.

An alternative would be to create a fishnet using each raster extent as template. Then use the label points created with fishnet and then run Extract Values to Points. The points holding values will be an approximation of the pixels extents. Select these and convert to polygons with minimum bounding geometry, buffer or point to polygons. I dont have Spatial Analyst extension so I cant try this.

import arcpy, os

rasterfolder =  r"X:\Rasters"
zones = r"X:\Zoneshape.shp"
zonefieldname = 'ZoneID'



arcpy.env.workspace = rasterfolder
rasterlist = arcpy.ListRasters()

#Store extent polygons in a list
polygons = []
for r in rasterlist:
    raster = arcpy.sa.Raster(r)
    polygons.append(raster.extent.polygon)


#Create a feature class of them and add an attribute with Raster filename    
arcpy.env.workspace = 'in_memory'
arcpy.CopyFeatures_management(in_features=polygons, out_feature_class='extents')
arcpy.AddField_management(in_table='extents', field_name="RasterID", field_type="TEXT")
p = iter(rasterlist)
with arcpy.da.UpdateCursor('extents', "RasterID") as cursor:
    for row in cursor:
        row[0]=p.__next__()
        cursor.updateRow(row)


#Intersect extents with zones 
arcpy.Intersect_analysis(in_features=['extents',zones], out_feature_class='Intersect')

#List results
l=[i for i in arcpy.da.SearchCursor('Intersect',['RasterID',zonefieldname,'SHAPE@AREA'])]

#Dictionary of raster name and largest zone name since the largest comes last in list and will overwrite the smaller ones.
d = {k[0]:k[1] for k in sorted(l, key=lambda x: x[0::2])}

#Rename the rasters using dictionary

for k in d:
    os.rename(src=os.path.join(rasterfolder,k),
              dst=os.path.join(rasterfolder,d[k]+'_'+k))