2

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

enter image description here

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

enter image description here

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.

  • 1
    Use IsNull to classify into data/nodata, this gives you a binary raster with 1 where the original is NoData and 0 where it is data. Set raster properties tool to set nodata to 1 on the isnull raster then when converting the isnull raster to polygon use simplify and it will speed up enormously. First make sure that the nodata value is set properly in your input raster. Are you doing this in a script? I see a lot of pictures (which helps) but no code.. what code do you have so far? – Michael Stimson May 4 '18 at 13:51
  • @MichaelStimson when I execute IsNull I get a binary raster with 0 in it which indicates data as far as I understand. I haven't written a script so far as I want to test the approach for just one image first and then incorporate it in a script by looping. What I noticed is that my data looks kind of strange: The Color Index in NoData areas is 0 and as well a field Color(a,r,g,b) 255, 0,0,0 is indicated when I use the Identify tab... – redshoe May 4 '18 at 14:25
  • Because it's black you haven't set the nodata value, this is important before running isnull, there is also setnull but that uses only one band, to set the null value on the raster use set raster properties and set the nodata value to RGB 0 which would be nodata="1 0;2 0;3 0". – Michael Stimson May 5 '18 at 10:51
3

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

enter image description here

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.