looping through rasters in an mxd

[RosieBell]

Hi all, newbie question here. I'd like to use the Python window to write a piece of code which loops through all the rasters in an mxd, and creates a polygon of the outline of each raster.

The Raster Domain tool in the 3D analyst extension works well for this task. However, when I try and use the code snippet in the Python window, I encounter the following problem. When I go to type in the In Raster ("dtm_grd") in the example below, I get a drop down list of all the rasters in my mxd. If I choose one from this list, I then get an error message that the raster "does not exist or is not supported". However, if I change the environment to the folder where the raster is currently located, and copy the source of the raster from the Properties box, the tool executes successfully.

Is there a way of looping through these layers, without having to manually change the address for each one? I know (very roughly) how to loop through a list of layers in python, but there doesn't seem much point trying if I can't get the Python snippet to work on even one layer.

nb not using the Raster to Polygon tool as most of the rasters are floating point and would need to be converted to integer before this tool would work.

arcpy.CheckOutExtension("3D") env.workspace = "C:/data" arcpy.RasterDomain_3d("dtm_grd", "raster_domain.shp", "POLYGON")

[NathanHeick]

This is my best answer.  I am not an expert, but I try to answer someone's question when someone else helps me.

I would try to first use desc = arcpy.Describe(rasterLayer) to describe each raster layer in the map document.  See below:

http://resources.arcgis.com/en/help/main/10.1/index.html#/Describe/018v00000066000000/

Then, I would try to use the desc.dataElement property of the layer describe object to get the raster dataset describe object.  See here:

http://resources.arcgis.com/en/help/main/10.1/index.html#/Layer_properties/018v00000063000000/

Lastly, I think you could use desc.dataElement.CatalogPath to get the source of the raster.  See here:

http://resources.arcgis.com/en/help/main/10.1/index.html#/Describe_object_properties/018v00000013000...

Then, you could use it directly as an argument to RasterDomain without setting the workspace.  It sounds like you have various source folders for the raster datasets and you care specifically about the rasters in an MXD.  I don't know how you were planning on listing the layers, but I am assuming you would have to use the arcpy.mapping module instead of arcpy.ListRasters().

Hope this helps.

[ModyBuchbinder]

Hi

Try something like this:

mxd = arcpy.mapping.MapDocument(r"C:\temp\test.mxd")
lst = arcpy.mapping.ListLayers(mxd,"*")
for l in lst:
       print l.isRasterLayer

Have fun

[DustinEdge]

Hi there

Once you call the Describe function on each raster, you can get the Extent from the Describe properties, which returns an array of the coordinates defining the extent. You can pass that as an array to create a new polygon geometry.

dsc = arcpy.Describe(raster)

xmin = dsc.Extent.XMin
ymin = dsc.Extent.YMin
xmax = dsc.Extent.XMax
ymax = dsc.Extent.YMax


then use the code in the Polygon documentation to create a new shape for each raster
http://resources.arcgis.com/en/help/main/10.1/index.html#/Polygon/018z00000061000000/

Hope this helps

Cheers
Dustin

[RosieBell]

Thanks to everyone for taking the time to help out. Unfortunately, I'm still not quite clear on how to achieve this - but this is a reflection on my skill level, rather than your help! I'll keep playing with it, and hopefully understand it soon. Cheers.

[T__WayneWhitley · Accepted solution]

This may help you, see the code...

http://forums.arcgis.com/threads/74002-Create-index-outline-of-rasters?p=258761#post258761

When I get around to it, I intend to update it to arcpy and improve it, make a more universal tool out of it....later, when there's more free time.  But the basic scripting is there.  Guess if you're still having trouble I could make a finished tool out of it over the weekend for you to use, since I've been meaning to do it anyway.

Enjoy,
Wayne

[RosieBell]

Hi Wayne, thanks very much for your reply. I actually got it to work - this is the longest piece of code that I've ever used so quite exciting.

I've posted the v10 version of your code below - all credit to Wayne and none to me btw! The only other differences apart from v9/ v10 stuff was that I designed this to be loaded from the Python window, for the shameful reason that I don't know how to do it any other way.

#Script to loop through rasters in an mxd, and print a shape file of their extent.
# Based on script at the following forum link: http://forums.arcgis.com/threads/74002-Create-index-outline-of-rasters?p=258761#post258761


mxd = arcpy.mapping.MapDocument("CURRENT")
GridName = "boundaries"
IMAGEfldLen = 50
arcpy.CreateFeatureclass_management("<file path for directory which will contain shapefile>", "<insert name of shapefile>", "POLYGON")
fields = ['IMAGE', 'XMIN', 'XMAX', 'YMIN', 'YMAX']
arcpy.AddField_management("<insert name of shapefile>", fields[0], "TEXT", '', '', IMAGEfldLen)

for i in range(1,5):
    arcpy.AddField_management(GridName, fields, "DOUBLE", 18, 17)
    print "added field: " + fields

rasters = arcpy.mapping.ListLayers(mxd)
arrayObj = arcpy.CreateObject('Array')
outRows = arcpy.InsertCursor(GridName)

for raster in rasters:
    extent = arcpy.Describe(raster).extent
    arrayObj.add(extent.lowerleft)
    arrayObj.add(extent.lowerright)
    arrayObj.add(extent.upperright)
    arrayObj.add(extent.upperleft)
    arrayObj.add(extent.lowerleft)
    feat = outRows.newRow()
    feat.Shape = arrayObj
    image = '.\\' + str(raster)
    feat.setValue('IMAGE', image)
    feat.setValue('XMIN', extent.xmin)
    feat.setValue('YMIN', extent.ymin)
    feat.setValue('XMAX', extent.xmax)
    feat.setValue('YMAX', extent.ymax)
    outRows.insertRow(feat)
    arrayObj.removeAll()

del outRows

nb in the resulting shapefile, I ended up with an extra row called 'boundaries' (name of my shapefile) which was a polygon which covered the entire area of any rasters. Not sure how I managed that, but I just deleted it anyway.

After all that effort, I've decided I really need to have a polygon that shows the exact boundary of each raster (not just a rectangle) so I'll keep working on it - if I achieve it I'll also post that code on this thread.

Oh and a quick qu for Wayne - why is the range 1 to 5? Cheers.

[T__WayneWhitley]

Hi Rosie, that's great you got it working!

The reason I indexed the list iteration 1-5 is so that the looping for the double type fields (to hold X- and Y- min/max coord vals) starts at index 1 (item 2 in the list) and ends at index 4 (item 5 in the list), meaning index 0 is skipped.
I reserved the index 0 position in the list for the text type field, IMAGE.  Although it could very well have been written differently to handle creating the fields all in the same loop (such as with a dictionary to 'look up' the field type) and maybe I should have done this if for no other reason than to make the code more 'readable', I simply handled the single text field creation in the line prior to entering the loop (fields[0] is IMAGE):

gp.AddField_management(GridName, fields[0], 'TEXT', '', '', IMAGEfldLen)

If this helps, this is how it would be to loop on the entire list:

>>> fields = ['IMAGE', 'XMIN', 'XMAX', 'YMIN', 'YMAX']
>>> for i in range(0, 5):  # or for i in range(5):
 print fields

 
IMAGE
XMIN
XMAX
YMIN
YMAX

...but, then, if accessing the entire list that way it could be done without specifying the index variable i, simply:

>>> for field in fields:
 print field

 
IMAGE
XMIN
XMAX
YMIN
YMAX

Hope that helps.  Again, I'm glad you found that useful!

Enjoy,
Wayne

[RosieBell]

Hi Wayne,

Thanks for that and I think I understand it now 🙂

BTW I did manage to get a script working for getting the exact outline for the rasters

#Script to loop through rasters in an mxd, and print a shape file of their extent (exact boundary, not just polygon).

import arcpy

from arcpy import env

arcpy.CheckOutExtension("3D")

mxd = arcpy.mapping.MapDocument(r"J:\gisprojects\Project\330\80000_89999\3308480_SS_GIS\0051_LiDAR_Index\mxd\ras_to_polygon_test.mxd")

env.workspace = r"J:\gisprojects\Project\330\80000_89999\3308480_SS_GIS\0051_LiDAR_Index\data\shpfile_ras_boundaries"
rasters = arcpy.mapping.ListLayers(mxd)

for raster in rasters:
    
    #Generate a unique name for each polygon
    outPoly = "domain_" + str(raster) + ".shp"
    #Execute Raster Domain Tool
    arcpy.RasterDomain_3d(raster, outPoly, "POLYGON")
print "finished"

Only problem now is that they are very big, I'm guessing because it's looking for 3D values? Each polygon is about 8MB, and I only need the outline. Does anyone know how I can alter this, in the environments perhaps....

Cheers
Rosie

[T__WayneWhitley]

Interesting, what does your output look like?  Yes, the Z component comes into play.  Depending on your rasters, there could be a great many more vertices that are stored...and the possibility of multipart features.  You could try generalizing if that more satisfactorily meets your end purpose.

Oh, also, switch to 'line' out_geometry_type for the output - undoubtedly this will save on space as well.

arcpy.RasterDomain_3d(raster, outPoly, "LINE")