2

I am a newcomer to ArcGIS and Python.

What I am trying to do is very simple: map 400K business locations using lat/lon coordinates onto a shapefile of California census tracts, and then export the data from the spatial join into a text file so that I have a spreadsheet containing the census tract that each business falls into. The input file of lat/lon data is a csv with 400K records, 5 variables and is ~15MB.

I have successfully completed this task interactively in ArcMap 10.2.1, but I would like to have a record of my process for future replication, hence my foray into ArcPy. In ArcMap, this process takes a few minutes. However, my ArcPy script takes multiple hours to complete. It appears to be the CopyFeatures step that requires significant computer resources and so I am not including ArcMap details or ArcPy code beyond this step. Note also that I have run my full ArcPy script on a subset of my data (the first 20 records of my 15MB file) and it works fine and finishes in just a few seconds. Given that, I know the issue is likely with the size of my file, but I have not been able to determine or successfully implement a work around to decrease the memory usage of my code.

Are there any ways I can modify my script to decrease processing time for this simple task?

I'm hopeful there's a way to make this as efficient with ArcPy as it is with ArcMap.

ArcMap steps:

  1. Import shapefile containing census tracts

    • Add data -> choose shapefile
  2. Import and prepare lat/lon data for spatial join

    • File -> Add Data -> Add XY Data -> select my 15GB csv -> X Field (longitude) -> Y Field (latitude) -> Coordinate System (WGS 1984)

    • Right click on XY data -> Data -> Export data -> All Features -> Output feature class is a .shp saved in my GIS workspace -> upon completion of this step, a pop up displays asking whether want to add to map as layer and I choose Yes (don't think this is significant however)

ArcPy script:

# Import modules
import arcpy, csv
import os
from arcpy import env

# Set environment settings
env.workspace = "C:/Users/ccwright/Desktop/GIS_workspace"
arcpy.env.overwriteOutput = True

# Set local variables:
in_table = "infousa_latlon_small.csv"
x_coords = "longitude"
y_coords = "latitude"
out_layer = "lat_lon_layer"
feature_class = "lat_lon_feature_class.shp"
census_tracts_shp = "census_tracts_CA_only/US_tract_2000_CA_only.shp"
spatial_join = "spatial_join.shp"
spatial_join_output_txt = "spatial_join_data.txt"

# Set coordinate system to WGS-1984:
WKID = 4326 
sr = arcpy.SpatialReference()
sr.factoryCode = WKID
sr.create()
env.outputCoordinateSystem = sr

# Process: Make XY Event Layer
arcpy.MakeXYEventLayer_management(in_table, x_coords, y_coords, out_layer, sr, "")

# Process: Copy Features
arcpy.CopyFeatures_management(out_layer, feature_class, "", "0", "0", "0")
  • You say that your problem occurs during the CopyFeatures step so I think you should edit your question to remove any steps subsequent to that. If you get past CopyFeatures in this Q&A then you can always research/ask about any subsequent problematic steps later. The other thing that I think you should edit in are the precise steps that you are performing in ArcMap so that we can see that your ArcPy script is accurately reproducing them. – PolyGeo Jul 20 '14 at 3:45
  • Thank you for your suggestions PolyGeo. I have removed my ArcPy code subsequent to the CopyFeatures step and added detail on my interactive steps in ArcMap. I look forward to hearing your thoughts! – ccwright Jul 21 '14 at 20:11
  • It looks like what needs to be compared to your ArcPy code is just steps 1 & 2. I'll try and find time to generate similar sized datasets and test the alternative procedures to see if I can reproduce the same performance you are reporting, but hopefully someone else can beat me to it. It looks like you are doing the same thing two ways so you have provided an excellent foundation from which to develop such a test. – PolyGeo Jul 21 '14 at 20:42
  • Thanks much, PolyGeo! I removed the steps beyond 1 & 2 for clarity. I appreciate the feedback and hope we can figure this thing out. I will continue looking for potential solutions in the meantime. Thank you! – ccwright Jul 21 '14 at 20:55
2

Try changing the CopyFeatures_management to FeatureClassToFeatureClass_conversion(). It also takes a feature layer and is essentially going to do the same thing you would when you right-click on your XY Event Layer in the TOC, go to Data and select Export.

Also, you could probably simplify setting your output spatial coordinate system by just directly referencing the spatial reference object. no need to create the object. The env object knows how to do all that stuff already.

# Set output coordinate system to WGS-1984:
env.outputCoordinateSystem = arcpy.SpatialReference(4326)
0

I suggest to change

arcpy.SpatialJoin_analysis(census_tracts_shp, feature_class, spatial_join, "JOIN_ONE_TO_MANY", "KEEP_ALL","", "INTERSECT", "", "")

to

arcpy.SpatialJoin_analysis(feature_class,census_tracts_shp, spatial_join, "JOIN_ONE_TO_ONE", "KEEP_ALL","", "INTERSECT", "", "")

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

  • Thank you for the suggestion. This modification takes about just as long as my original script unfortunately. Also, I specified the census tract shapefile as the target features and the lat/lon data as the join features because I want my output table to contain rows for all census tracts, even those that do not have any points mapped to them. Am I misunderstanding the need to keep my code as I originally had it in order to get this result? – ccwright Jul 21 '14 at 23:44
  • I am just struggling to understand the output you'd like to see, because you'll end up with n polygons (census) each reffering to 1 out of n businesses inside that polygon. If this is what you realy want to see disregard my solution. Anyway, very long paths and especially shapefile names are not looking good to me, because shapefiles (.shp) has associated .dbf file. dBase file is very old format from DOS era, when file name length was limited to 8 characters. Try to reduce "lat_lon_feature_class.shp" to "l_l_fc.shp" and see what happen – FelixIP Jul 22 '14 at 0:24
  • Essentially the only output I need is a table that has attached a census tract code to each of my lat/long records. It seems I need to specify the tract polygons as the target, the lat/long data as the join and choose a one to many join because otherwise, using the method you suggest, I only get a table containing a single lat/long record in each census tract. Do you know if there are any alternative mapping procedures I could do that would circumvent the whole shapefile creation, etc and just give me the underlying attribute table? This should certainly reduce processing time significantly... – ccwright Jul 24 '14 at 20:58
  • And thank you for the suggestions about naming conventions! I'm not sure if it reduced processing time in this instance but it is definitely very valuable to know in general! – ccwright Jul 24 '14 at 21:03
  • "Essentially the only output I need is a table that has attached a census tract code to each of my lat/long records." Great at last, because this is exactly what I realised you were trying to achieve, despite you saying the opposite. And yes, with my method if you have 5000 businesses and 200 tracts you'll get 5000 records, i.e. one for each business with information on parent track, and not 200 records as you think it will.. – FelixIP Jul 24 '14 at 22:37

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.