4

The code below takes a huge raster, and a set of point features. For each of them creates a circular buffer (100), then squares it, and then clips that part of the buffer. In doing so, we end up with a png image for each of the point locations.

The script works fine, the main problem, is that it adds each clipped feature to the workspace. These are around 35,000 points, so at about 1100, it just kills arcmap, and I have to restart the process from where it stopped (that is the reason there is an if, to check for 1107).

Is there a way to just do the clipping, without everything being loaded in the environment?

import arcpy
import os
image_name = "MyMosaicDataset2.tif" #has to be loaded in the env
main_path = 'C:\Users\leon\Documents\Data\VCR_Data\Leon'
workspace_name = 'Mars_Tartarus_Colles_Georefrencedg.mxd'
arcpy.env.workspace = os.path.join(main_path, workspace_name)
base_for_buffer = 'VRCs'
buff_name = 'VRC_buff'
output_buffer = os.path.join(main_path, buff_name+'.shp')
output_sq_buffer = os.path.join(main_path, buff_name+'_sq'+'.shp')
radius_buffer = '100 meters'
arcpy.Buffer_analysis(base_for_buffer, output_buffer , radius_buffer, 'FULL', 'ROUND', 'NONE','','PLANAR') #This is to create the buffers
arcpy.FeatureEnvelopeToPolygon_management(output_buffer, output_sq_buffer,"SINGLEPART") #This transforms the circular buffers into rectangular buffers
rows = arcpy.SearchCursor(buff_name+'_sq')#here are the buffer polygons, and stored ina  list
#This is for clipping
cnt=1
for row in rows: #go over the list
    if cnt >= 1107:
        feat = row.Shape
        name = 'vrcp_'+str(cnt)
        output_raster = os.path.join(main_path, name)
        arcpy.Clip_management(image_name, "#", output_raster+'.png',feat, "0", "ClippingGeometry")
    cnt+=1
    enter code here
3
  • is base_for_buffer referring to a layer within your MXD? I think if you change that to refer directly to a feature class, and set your arcpy.env.workspace to not reference your MXD (reference your folder instead), then the outputs shouldn't be added to your MXD.
    – Midavalo
    Jan 26, 2016 at 0:46
  • Just wonder why you are not using Split Raster tool, with split method POLYGON_FEATURES?
    – FelixIP
    Jan 26, 2016 at 1:12
  • You don't need to run the script in ArcMap if you don't interact with the map features - run it from IDE or from a CMD; use arcpy.da.SearchCursor if you are on 10.1 SP1+, 10.0 SearchCursor is much slower; ArcMap menu Geoprocessing > Geoprocessing Options > uncheck "Add results of geoprocessing operations to the display" - no newly created datasets will be added to ArcMap TOC. Jan 26, 2016 at 7:46

2 Answers 2

2

Open mxd document and go to Geoprocessing->Geoprocessing Options-> and uncheck second checkbox from the bottom "Add results of geoprocessing operations to the display". I don't know if this option could be set by arcpy.

I am probably not aware of all circumstances, but considering your code, I think you really don't need mxd as a workspace, exactly as Midavalo suggested in his comment. I mean if you are not manipulating layers symbology, layers order, map layout, or any graphical feature, you can use databases, or folders as a workspaces instead of mxd.

1

The clipped images will not load into the data frame if you run the script as a script tool. This is because you don't create any layers in your script.

I happened to have a raster layer and a point layer handy, so I tried your script to make sure. And it worked as I expected.

I am attaching my version of your script that worked with some notes of explanation. Everywhere it says NOTE means I had to change something to get the script to work as you wanted.

import arcpy
import os
image_name = "MyMosaicDataset2.tif" 
main_path = r'C:\Users\leon\Documents\Data\VCR_Data\Leon' # NOTE path syntax - add r
arcpy.env.workspace = main_path # NOTE an MXD is not a valid workspace
base_for_buffer = 'VRCs' # point layer being buffered
buff_name = 'VRC_buff' # base filename for both shapefiles
output_buffer = os.path.join(main_path, buff_name+'.shp')
output_sq_buffer = os.path.join(main_path, buff_name+'_sq'+'.shp')
radius_buffer = '100 meters'

# Create point buffers as polygon feature class
# NOTE original script had error - too many parameters, path error
# NOTE also that curly braces mean an optional parameter
# don't need to specify dissolve field if not dissolving.
# Buffer_analysis (in_features, out_feature_class, buffer_distance_or_field, {line_side}, {line_end_type}, {dissolve_option}, {dissolve_field})
arcpy.Buffer_analysis(base_for_buffer, output_buffer, radius_buffer, 'FULL', 'ROUND')

# Turn the point buffers into squares by getting their extents
# output to polygon feature class
# FeatureEnvelopeToPolygon_management (in_features, out_feature_class, {single_envelope})
arcpy.FeatureEnvelopeToPolygon_management(output_buffer, output_sq_buffer, "SINGLEPART")

rows = arcpy.SearchCursor(output_sq_buffer) # NOTE this creates a feature class cursor, not a list
                                            # NOTE use entire path when there is not a layer
cnt=1
for row in rows: # loops through squared polygon features
    if cnt >= 1:
        feat = row.Shape # get the square shape
        name = 'vrcp_'+str(cnt) # output raster base name
        output_raster = os.path.join(main_path, name) 
        # raster clip management(clip raster, rectangle, output raster, template dataset for extent, no data value, clipping_geometry)
        # NOTE use "" for NULL parameter when creating script tool
        arcpy.Clip_management(image_name, "", output_raster+'.png', feat, "0", "ClippingGeometry")
    cnt+=1
del rows # make sure you delete your cursor when done with items
# and release your other variables too
del image_name, main_path, base_for_buffer, buff_name
del output_buffer, output_sq_buffer

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.