2

I am trying to select and copy shapefiles from one folder to a new location based on if they intersect a polygon feature class. I am getting stuck in somewhere in the for loop.

EDIT: The goal is to copy the whole file if the extent of the shapefile intersects the polygon. I am not wanting to copy a subset of any shape.

# Generate list of shapefiles in the workspace
# Copy shapefiles in a list that intersect a polygon to a new location

import arcpy
import os


set environment settings


arcpy.env.workspace = 'E:\workspace\Contour_Shapefile'


set local variables


out_workspace = 'C:\localvariables\outworkspace\'
select_features = 'C:\localvariables\Project.gdb\AreaOfInterest'


use ListFeatureClasses to generate a list of shapefiles in the workspace shown above.


fc_list = arcpy.ListFeatureClasses()


#Execute CopyFeatures for each input shapefile that INTERSECT Area of interest.
i=1
for shapefile in fc_list:
    out_featureclass = os.path.join(out_workspace, os.path.splitext(shapefile)[0])
    arcpy.SelectLayerByLocation_management(shapefile, 'INTERSECT', select_features)
    arcpy.CopyFeatures_management(shapefile, out_featureclass)
    print(str(i) + ' copied ' + shapefile)
    i += 1
else:
    print( str(i) + ' MISSING')
    i += 1
inpfile.close()
axman2
  • 63
  • 4
  • You state you want to "select and copy shapefiles from one folder to a new location based on if they intersect a polygon feature class". Your logic is selecting a subset IN a shape file and copying those features. You are not copying the whole file. What is it you are actually trying to achieve? – Hornbydd Dec 30 '20 at 00:16
  • Hi, I am wanting to copy the whole file if the extent of the shapefile intersects the polygon. I am not wanting to copy a subset of any shape file. – axman2 Dec 30 '20 at 00:51
  • Here is a similar q&a using rasters. The original poster fell into the trap of different coordinate systems, have you accounted for that? – Hornbydd Dec 30 '20 at 01:02
  • Thanks, I hacked the similar post, linked to in your previous comment, and it works as well. Though it is only copies shapefiles that are completely within the AOI. Coordinate systems is not the issue with my question. – axman2 Dec 30 '20 at 21:08

1 Answers1

4

I see the very last line as error-ing out as inpfile is not defined anywhere. You also didn't store the selection set as a variable or test whether the returned selection layer had any records before copying.

try this.

# Generate list of shapefiles in the workspace
# Copy shapefiles in a list that intersect a polygon to a new location

import arcpy, os


set environment settings


arcpy.env.workspace = 'E:\workspace\Contour_Shapefile'


set local variables


out_workspace = 'C:\localvariables\outworkspace\'
select_features = 'C:\localvariables\Project.gdb\AreaOfInterest'


use ListFeatureClasses to generate a list of shapefiles in the workspace shown above.


fc_list = arcpy.ListFeatureClasses()


#Execute CopyFeatures for each input shapefile that INTERSECT Area of interest.
for shapefile in fc_list:
  #create output feature class
  out_featureclass = os.path.join(out_workspace, os.path.splitext(shapefile)[0])
  #create selection layer in memory
  selection = arcpy.SelectLayerByLocation_management(shapefile, 'INTERSECT', select_features)
  #retreieve record count of selection
  selectionCount = int(arcpy.management.GetCount(selection)[0])
  #print shapefile and record count
  print(f"Selection Shapefile: {shapefile}; Record Count: {selectionCount}")
  #if selection contains records, copy feature
  if selectionCount > 0:
  #on the following line, if you replace 'shapefile' with 'selection' it will copy just the intersecting features.
    arcpy.CopyFeatures_management(shapefile, out_featureclass)
    print(f"Copied Shapefile: {shapefile}")
  else:
    print(f"Skipping shapefile -> {shapefile}")

I stored your selection layer as a variable, tested the selection layer to see if there were records, and if there were more than 0 records in the selection layer, copied the layer into the new directory. This has f-string formatting, so you may need to comment out my print statements if using an older version of python. There are also comments on each line explaining what I did. If you need more help or this doesn't work let me know and let's see if we can't figure it out!

EDIT: this can also be done with the Arcpy Geometry Class. Look at the "disjoint" and "distanceTo" methods. Disjoint returns false if 2 geometries intersect, and distanceTo returns the minimum distance between two objects (returns 0 if they intersect) https://pro.arcgis.com/en/pro-app/latest/arcpy/classes/geometry.htm

John
  • 607
  • 3
  • 10
  • Thanks for the help, John. This method works, but because of the size and number of the shapefiles the script has to process, it takes a few hours to complete (>6hrs). I am curious about using the Arcpy Geometry Class method and how it could potentially reduce the time the script would have to run. Would using the disjoint method replace creating a selection layer? – axman2 Dec 30 '20 at 20:44
  • Nathanael, I am not sure. I haven't used arcpy geometries all that much. If you are looking to speed things up you could try alternatives like running a clip function instead of a search by location or intersection (i believe the clip should run faster.) otherwise I would look into other python libraries like fiona/shapely. Check out these answers: https://stackoverflow.com/questions/25169707/python-shapely-touching-contains and https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python. You could also look into using sql geometry if you are using a SQL database. – John Dec 30 '20 at 21:55