POST
|
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)
... View more
01-21-2013
08:05 PM
|
0
|
0
|
176
|
POST
|
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.
... View more
01-21-2013
04:11 PM
|
0
|
0
|
176
|
POST
|
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
... View more
01-20-2013
08:45 PM
|
0
|
0
|
660
|
POST
|
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.
... View more
01-20-2013
06:17 PM
|
0
|
0
|
660
|
POST
|
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.
... View more
01-17-2013
06:23 PM
|
0
|
0
|
660
|
POST
|
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")
... View more
01-13-2013
05:16 PM
|
0
|
12
|
2256
|
POST
|
Rosie, I can reproduce this. A fix that works for me is to reload arcpy: >>> reload(arcpy) Thanks Curtis! I tried this, and it works perfectly.
... View more
12-03-2012
03:25 PM
|
0
|
0
|
422
|
POST
|
Bumping this up as it's still an issue for me. Anyone else experience this, or find a solution? I'm sure this is ridiculously simple, but I can't get it to work. I'm using the arcpy.Clip_management tool. When I type the following into the window: arcpy.Clip_management(" A drop down box appears, showing all the raster layers that exist in my mxd. I can select one and then complete the rest of the script. However, if I've started using the python window, then add another layer to my mxd, it doesn't appear in the down down box when I reach this point in the script. It's almost like it needs a refresh to appear. At the moment the only thing I've found to work is restarting ArcGIS. The attached image shows 3 rasters in the table of contents - bunbury, cli7, and perth - but only the original two come as options in the Python window. TIA.[ATTACH=CONFIG]18634[/ATTACH]
... View more
11-29-2012
04:55 PM
|
0
|
0
|
422
|
POST
|
I've been experimenting with different methods of clipping a raster in ArcGIS 10. Extract by Mask seems to be the preferred option - "unless you don't have the spatial analyst extension". I do have the extension, but it seems very slow. Is there any advantage to using this to clip a raster, rather than Clip in Data Management, or exporting the data (after drawing the area with the draw tools). Both the second methods seem quick and effective. The raster I'm clipping is a DEM (GRID file) if that makes a difference. TIA Rosie
... View more
11-25-2012
09:01 PM
|
0
|
1
|
1587
|
POST
|
You've probably checked this - but you are not currently editing this layer, or another in the same location? You can't add fields if you are in an edit session.
... View more
11-22-2012
08:44 PM
|
0
|
0
|
251
|
POST
|
Hi Mathew, I'm running ArcGIS Desktop 10 SP3. I do start a new line after adding. I've also tried clearing all my python, closing and reopening the window - no luck though.
... View more
10-22-2012
10:17 PM
|
0
|
0
|
422
|
POST
|
I'm sure this is ridiculously simple, but I can't get it to work. I'm using the arcpy.Clip_management tool. When I type the following into the window: arcpy.Clip_management(" A drop down box appears, showing all the raster layers that exist in my mxd. I can select one and then complete the rest of the script. However, if I've started using the python window, then add another layer to my mxd, it doesn't appear in the down down box when I reach this point in the script. It's almost like it needs a refresh to appear. At the moment the only thing I've found to work is restarting ArcGIS. The attached image shows 3 rasters in the table of contents - bunbury, cli7, and perth - but only the original two come as options in the Python window. TIA.[ATTACH=CONFIG]18634[/ATTACH]
... View more
10-22-2012
07:16 PM
|
0
|
5
|
857
|
POST
|
I didn't realise there was a button to do this, so I added to my toolbar to see what it would do. It is always greyed out for me, too. I think it's because ArcGIS doesn't recognise an active layer (even if selected in TOC) so the button doesn't know which layer to open. Possibly a redundant feature from earlier versions? The ArcGIS 10.1 help files don't make any mention of it. http://resources.arcgis.com/en/help/main/10.1/index.html#//005s00000002000000 The keyboard shortcut CTRL_T works to open the selected layer in the TOC - this may help.
... View more
10-07-2012
06:48 PM
|
0
|
0
|
1409
|
POST
|
Sorry for the slow reply. I tried both solutions shortly after they were posted, but couldn't get either of them to work. I tried again jm's solution again today and it did the job, many thanks. carrierkh, no doubt it was my inexperience that was the problem, so thankyou for your help.
... View more
09-27-2012
11:12 PM
|
0
|
0
|
249
|
POST
|
I made my first Python tool this week, a Search by Attribute tool to search our cadastre. When I have required parameters it works beautifully. As soon as I make the parameters optional, it doesn't return any results (even if I complete all the fields). What am I doing wrong? I added the 'if' line after reading another thread about optional parameters, but it doesn't make any difference with this problem. Thanks. #define user parameter for road name, and locality SelCondition = arcpy.GetParameterAsText(0) if (not SelCondition) or (SelCondition == "#") or (len(SelCondition.strip()) == 0): SelCondition = "" LocalityName = arcpy.GetParameterAsText(1) if (not LocalityName) or (LocalityName == "#") or (len(LocalityName.strip()) == 0): LocalityName = "" # select cadastral lots on road name and locality, and lot no arcpy.SelectLayerByAttribute_management("Cadastre","NEW_SELECTION", "\"ROAD_NAME\" = "+SelCondition) arcpy.SelectLayerByAttribute_management("Cadastre","SUBSET_SELECTION", "\"LOCALITY\" = "+LocalityName)
... View more
08-20-2012
11:42 PM
|
0
|
3
|
2001
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:23 AM
|