nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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

Trouver plus de publications taguées avec

Accepted answers

T__WayneWhitley

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

All comments

RosieBell

Ah, I'm with you now. I hadn't realised that my polygon (and line) had so many vertices - about 50,000! Simplifying brought it down to under 20, much more manageable. I think I do need the polygons so have gone back to using POLYGON in the Raster Domain tool, and then simplifying.

This is my script to date (just in case anyone is trying something similar). It's a bit clumsy but it works. I included Delete_management to get rid of my large polygons, once I had a simplified version.

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

import arcpy
from arcpy import env

#disable M and Z values in output
arcpy.env.outputZflag = "Disabled"
arcpy.env.outputMflag = "Disabled"

arcpy.CheckOutExtension("3D")

mxd = arcpy.mapping.MapDocument(<insert path to mxd here>)

env.workspace = <insert path to workspace here>(mxd, "_*")

for raster in rasters:
    
    #Generate a unique name for each polygon/ shapefile. Have given them the prefix "del" so I can find them and delete them later
    outPoly = "del" + str(raster) + ".shp"
    #Execute Raster Domain Tool
    arcpy.RasterDomain_3d(raster, outPoly, "POLYGON")

#get the feature classes from this folder
fcList = arcpy.ListFeatureClasses()

for fc in fcList:
    #Generate a unique name for each simplified polygon
    simpPoly = fc.lstrip("del_")+ "simp_" + fc.rstrip(".shp")
    #Execute Simplify Polygon tool
    arcpy.SimplifyPolygon_cartography(fc, simpPoly, "POINT_REMOVE", 1, "", "", "NO_KEEP")

#get feature classes for deletion from same folder
deleteBigPolys = arcpy.ListFeatureClasses("del*")

#delete these features classes (too large)
for deleteBigPoly in deleteBigPolys:
    arcpy.Delete_management(deleteBigPoly)

T__WayneWhitley

Well, depends what you need it for. I would say now that you have it in line form, you can experiment with some of the generalization tools....because if you don't need vertices so 'densely packed', you can make some adjustments now, saving time in the conversion to polygon process - should you decide you need polygons.

Don't know how well this will work, but it may be worth looking into:

Simplify Line (Cartography)
Desktop » Geoprocessing » Tool reference » Cartography toolbox
http://resources.arcgis.com/en/help/main/10.1/index.html#/Simplify_Line/007000000010000000/

Just curious, but did you know you could export that extent raster polygon fc to a table and use that result in ArcMap as a 'legacy' raster catalog? In other words, in table form, ArcMap reads and recognizes the IMAGE, XMIN, YMIN, XMAX, YMAX format as a catalog and will auto-load your image tiles (in the view extent - and this is adjustable). I find that valuable sometimes, esp. when sharing multiple images without preprocessing a mosaic, etc. That is, 10.0 still 'recognizes' it; I don't know about 10.1, but I wouldn't be surprised if it still is supported.

Enjoy,
Wayne

RosieBell

I added this to the beginning of my code to disable M and Z values:

#disable M and Z values in output
arcpy.env.outputZflag = "Disabled"
arcpy.env.outputMflag = "Disabled"

unfortunately it didn't affect the size.

I also tried using the tool manually and disabling everything I could think of in the environments - again, didn't make any difference.

Saving it as LINE rather than POLYGON did the trick though! Thanks Wayne 🙂 Now I wonder if I should bother trying to convert this line to polygon through my script, or if it's not necessary.

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

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

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

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.

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

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

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/018v00000013000000/

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.

Rejoindre la conversation

Les membres connectés peuvent publier, suivre les mises à jour, et plus encore. Nouveau ici ? Inscrivez-vous gratuitement.

Liens rapides

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.