POST
|
I did finish the tool I was working on. The tool is faster than the individual polygon solution. Originally I produced 3 different solutions for breaking up the polygons and ended up using two of those methods. I have uploaded the tool and documented it. The tool uses Spatial Join and also breaks up the feature class into smaller pieces. http://www.arcgis.com/home/item.html?id=b859b33c616a47d2b99b5e133942db02 Hi Sue, Thanks for sharing your code! I was wondering though, would you mind explaining the logic behind the following sentence taken from your official tool post (your link above) under Description > Solution: "If a feature class is spatially joined to itself, the polygon id in the target and join will match in the top-most layer of polygons. This set of polygons will not overlap any other polygons." I am specifically confused by the phrase "top-most layer of polygons," and how this process produces non-overlapping features. I tried manual testing to better understand it, and when I spatially joined a feature class to itself, the output was identical to the original feature class (with the exception of added fields from the join). Further, joining this back to the original feature class on polygon id did nothing but add more fields - as all the polygons matched/joined because they were all the same as the originals. I have reviewed your code and am still stumped by this process. If you have the time and wouldn't mind going into the logic details, I'd love to understand what's happening here. I have been trying to implement a similar routine (but not using spatial join) that results in grouping points (via a flag field) into non-overlapping subsets based on either their buffer or distance to other points. I'm needing it for the Tabulate Area tool, which has the same overlapping areas problem as the Zonal Stats as Table tool. Any help would be appreciated! Thanks for your time! Abby
... View more
10-18-2012
12:24 PM
|
0
|
0
|
818
|
POST
|
Hi Sue, Thanks for sharing your code! I was wondering though, would you mind explaining the logic behind the following sentence taken from your official tool post under Description > Solution: "If a feature class is spatially joined to itself, the polygon id in the target and join will match in the top-most layer of polygons. This set of polygons will not overlap any other polygons." I am specifically confused by the phrase "top-most layer of polygons," and how this process produces non-overlapping features. I tried manual testing to better understand it, and when I spatially joined a feature class to itself, the output was identical to the original feature class (with the exception of added fields from the join). Further, joining this back to the original feature class on polygon id did nothing but add more fields - as all the polygons matched/joined because they were all the same as the originals. I have reviewed your code and am still stumped by this process. If you have the time and wouldn't mind going into the logic details, I'd love to understand what's happening here. I have been trying to implement a similar routine (but not using spatial join) that results in grouping points (via a flag field) into non-overlapping subsets based on either their buffer or distance to other points. I'm needing it for the Tabulate Area tool, which has the same overlapping areas problem as the Zonal Stats as Table tool. Any help would be appreciated! Thanks for your time! Abby
... View more
10-18-2012
11:59 AM
|
0
|
0
|
818
|
POST
|
Thanks Matthew, Yes, my problem is that nothing is ever removed from my selected features. I have tried changing the selection type (for example, to "WITHIN" the buffer, such that I'm selecting from a selection instead of removing from a selection), and it still does not work. Below is a print out of the first few loops running on a subset of my data. If this doesn't help I can post the data. Meters vs. kilometers should make no difference (10 Kilometers vs. 10000 Meters). The buffers are being created successfully (with proper distance value - proven when I open the shapefile in ArcMap), so the units do not seem to be an issue. I agree that this might be a convoluted work flow, but I'm not sure how else to approach it. My real issue is this: I need to retrieve the area (in square units) of each land-cover type/category within each feature of a polygon feature class. In my case, polygon features = buffers. The "Tabulate Area" tool (Spatial Analyst > Zonal) does this, however there is a problem inherent in the tool: if any polygons overlap, one will get TabArea results for the full polygon while the other(s) will only get results for part of the polygon (the part that does not overlap). I need areas for each full polygon. My original approach (and the approach of my predecessor) was to loop through each feature of interest (i.e. each polygon/buffer) one at a time, and run the TabArea tool individually on each feature. This is quite taxing on implementation time - especially when you're dealing with 35,000+ features. Originally written in VBA, this process took 3-4 seconds per feature. I've gotten Python to mimic the VBA code and processing time (since ESRI is phasing-out VBA), but at 35,000 records, it will take 29-39 hours to complete a single run - and I need to run this process on multiple buffer sizes and multiple land-cover datasets (approx. 35 times). While this code process does work, it's too inefficient for my current needs. My predecessor and I decided that if we could place the features into non-overlapping groups, the tool could run on each group instead of each feature (drastically reducing implementation time). This is the code I included with my original post. Logically, it works (albeit, convoluted), and it can be used with any buffer size of interest. However, as you're aware, I'm having unusual issues with the second selection set. Like I originally posted, this process works in ArcMap - both with the selection tool interfaces and the Python window. Let me know your thoughts, and thanks for your help! BEGINNING STEP 1: GENERAL SETTINGS... Start time: Fri Oct 12 09:42:10 2012 Step 1 COMPLETE! Time elapsed: 14.45 seconds BEGINNING STEP 2: Total (n) Records = 80 Buffer Count = 0 Record Count where flag is zero = 80 Record Count where flag is zero and wells are outside buffers = 80 UnqID = GTC_NO384734 Buffer Count = 1 Record Count where flag is zero = 79 Record Count where flag is zero and wells are outside buffers = 79 UnqID = GTC_NO385488 Buffer Count = 2 Record Count where flag is zero = 78 Record Count where flag is zero and wells are outside buffers = 78 UnqID = GTC_NO323237 Buffer Count = 3 Record Count where flag is zero = 77 Record Count where flag is zero and wells are outside buffers = 77 UnqID = USGS7073 Buffer Count = 4 Record Count where flag is zero = 76 Record Count where flag is zero and wells are outside buffers = 76 UnqID = GTC_NO357868 Buffer Count = 5 Record Count where flag is zero = 75 Record Count where flag is zero and wells are outside buffers = 75 UnqID = PWTS3399 Buffer Count = 6 Record Count where flag is zero = 74 Record Count where flag is zero and wells are outside buffers = 74 UnqID = PWTS753 Buffer Count = 7
... View more
10-12-2012
06:26 AM
|
0
|
0
|
614
|
POST
|
(Note 1: I'm using Arc 10.1 and PyScripter) (Note 2: Thank you, thank you, thank you in advanced for reading the monster explanation below! I am extremely frustrated and would greatly appreciate any help!) I'm not sure how to briefly sum-up my issue, but (very generally) I'm trying to assign values to a flag field in a point shapefile, where the flag field designates groups of points, and each group represents a collection of points that can be buffered to a given buffer size of interest without having their buffers overlap. So basically, I'm sub-setting the point dataset into groups that, when buffered, have buffers that won't overlap. My logic follows a loop whereby: 1) All non-flagged features are selected, 2) Any features falling within "CombinedBuffers.shp" (see definition below) are removed from the selection, 3) First feature in remaining selection is chosen (i.e., basically random), 4) This feature is flagged (for ex., "1" - for group 1), 5) Then this feature is buffered, and the buffer is added (appended) to a "combined buffers" shapefile (each buffer added to this shapefile represents an area under which additional "group 1" points cannot fall; the shapefile is emptied when a new group starts), 6) The code loops back to the top again - selecting all points not flagged, removing those under buffered areas, selecting the first feature, etc. 7) Once all possible points have been designated as group 1, CombinedBuffers.shp is emptied and the loop starts again on the subset of points not flagged (this time, assigning points to group 2, etc.) My problem is introduced in #2 above: removing points from the current selection based on location. (Note: Selection in #1 works fine.) No error is thrown for selection #2, it simply does not remove any points when it should (proven by record count). What I've already tried/checked: 1) "Input" and "CombinedBuffers" datasets are shapefiles, but I've also tried them as feature classes in a personal geodatabase. 2) All data are in the same projected coordinate system; all data are stored in the same location. 3) Selection #1 followed by selection #2 works manually in ArcMap with tool interfaces AND the python window (using the same code!). 4) The appended buffer in CombinedBuffers.shp IS being recognized (feature counts are correct after each added buffer). 5) I have changed the type of selection, and it still will not select properly. 6) I have used similar selection sets in previous codes (and they worked), but I have not tested them using 10.1 yet. Could this be a 10.1 issue? My code is below (Note: sorry for the lengthiness - anyone willing to read through this is highly commended and your help incredibly appreciated!):
#STEP 1: GENERAL SETTINGS-------------------------------------------------------
# Imports...
import time
print "Start time:",time.ctime()
StartSecs = time.clock()
import Functions
import arcpy
import sys
import traceback
import easygui
from arcpy import env
from arcpy.sa import *
# Variables...
BuffDist = 10
BuffUnit = "Kilometers"
Input = "C:/.../TestWells_TabAreaInput.shp"
Output = "C:/.../TestWells_TabAreaOutput.shp"
CmbdBuffs = "C:/.../CmbdBuffs.shp"
TempBuff = "C:/.../TempBuff.shp"
SpatialRef = arcpy.Describe(Input).spatialReference
Scratch_Path = "C:/.../Scratch"
CursorExists = "No"
RowExists = "No"
# Environments...
env.workspace = Scratch_Path
env.scratchWorkspace = Scratch_Path
env.overwriteOutput = "True"
env.extent = "C:/.../FC_Poly_IAState"
env.mask = "C:/.../FC_Poly_IAState"
env.maintainSpatialIndex = False
# Enable error handling...
try:
# STEP 2: CREATE POINTS FEATURE LAYER AND BUFFER FEATURE CLASS------------------
# Make 'wells' feature layer and set flag field = 0...
arcpy.MakeFeatureLayer_management(Input,"InputLyr")
arcpy.CalculateField_management("InputLyr","BuffSelFlg",0)
n = int(arcpy.GetCount_management("InputLyr").getOutput(0))
print "Total (n) Records =", n
# Create empty buffer feature class...
Buffer = str(BuffDist*2) + " " + BuffUnit # For these processes, buffers are defined as 2x buffer size of interest.
if arcpy.Exists(CmbdBuffs):
arcpy.Delete_management(CmbdBuffs)
arcpy.Buffer_analysis("InputLyr",CmbdBuffs,Buffer)
arcpy.DeleteFeatures_management(CmbdBuffs)
BufferCnt = int(arcpy.GetCount_management(CmbdBuffs).getOutput(0))
print "Buffer Count =", BufferCnt
x = 1 # (assigned group value/flag)
# STEP 3: START LOOP------------------------------------------------------------
# Note: The use of "group" below is defined as a group of wells that can be buffered
# to the buffer size of interest without having their buffers overlap.
while x <= n:
# Select wells that have not been added to a group (i.e., where flag = 0)...
arcpy.SelectLayerByAttribute_management("InputLyr","NEW_SELECTION",'"BuffSelFlg" = 0')
RcdCnt = int(arcpy.GetCount_management("InputLyr").getOutput(0))
print "Record Count where flag is zero =", RcdCnt
if RcdCnt > 0: # (i.e., if there are wells that have not been added to any group...)
# Remove (from selection) wells within the current group's buffers...
arcpy.SelectLayerByLocation_management("InputLyr","INTERSECT",CmbdBuffs,"","REMOVE_FROM_SELECTION")
OutsideBuffRcdCnt = int(arcpy.GetCount_management("InputLyr").getOutput(0))
print "Record Count where flag is zero and wells are outside buffers =", RcdCnt
# NOTE: THIS IS THE SELECTION THAT DOES NOT WORK
if OutsideBuffRcdCnt > 0: # (i.e., if there are wells left that are qualified for assignment into this group...)
UpdtCur = arcpy.UpdateCursor("InputLyr","",SpatialRef)
row = UpdtCur.next()
row.BuffSelFlg = x
UpdtCur.updateRow(row)
UnqID = row.UnqRecordI
print "UnqID =",UnqID
arcpy.SelectLayerByAttribute_management("InputLyr","NEW_SELECTION","\"UnqRecordI\" = '" + UnqID + "'")
arcpy.Buffer_analysis("InputLyr",TempBuff,Buffer)
arcpy.Append_management(TempBuff,CmbdBuffs,"NO_TEST")
arcpy.AddSpatialIndex_management(CmbdBuffs)
BufferCnt = int(arcpy.GetCount_management(CmbdBuffs).getOutput(0))
print "Buffer Count =", BufferCnt
arcpy.SelectLayerByAttribute_management("InputLyr","CLEAR_SELECTION")
arcpy.Delete_management(TempBuff)
del UpdtCur
del row
else: # (Start a new group)
arcpy.SelectLayerByAttribute_management("InputLyr","CLEAR_SELECTION")
arcpy.DeleteFeatures_management(CmbdBuffs)
x = x + 1
else: # (All records have been flagged into groups - job done!)
break
# STEP 4: COPY FEATURES AND DELETE FEATURE LAYER----------------------------
arcpy.CopyFeatures_management("InputLyr",Output)
arcpy.Delete_management("InputLyr")
# ERROR HANDLING----------------------------------------------------------------
except arcpy.ExecuteError:
# Get the tool error messages:
msgs = arcpy.GetMessages(2)
# Return tool error messages for use with a script tool:
arcpy.AddError(msgs)
# Print tool error messages for use in Python/Python Window:
print "Geoprocessing tool error:"
print msgs
except:
# Get the traceback object:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
# Concatenate info together concerning the error into a message string:
pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"
# Return python error messages for use in script tool or Python Window:
arcpy.AddError(pymsg)
arcpy.AddError(msgs)
# Print Python error messages for use in Python/Python Window:
print arcpy.GetMessages()
print pymsg
finally:
if CursorExists == "Yes":
del InputCur
if RowExists == "Yes":
del row
if arcpy.Exists("InputLyr"):
arcpy.Delete_management("InputLyr")
print "End time:",time.ctime()
TimeElapsed, Step = (time.clock()-StartSecs), "END"
Functions.PrintTime(TimeElapsed,Step)
... View more
10-11-2012
12:00 PM
|
0
|
3
|
1605
|
POST
|
Huh, yea that's really interesting. The icon by your name means you work for ESRI, right? Can you report this bug or should I? Abby
... View more
08-08-2012
10:58 AM
|
0
|
0
|
633
|
POST
|
Thanks! When I use MakeFeatureLayer and include field restrictions, however, it does not recognize any of my fields. In the code below, all of the "print FieldInfo.getFieldName()" lines return the correct fields/names. When I try to access any of these fields within the cursor though, it throws an error: "Row: Field FID does not exist" (or "Row: Field ___ does not exist" - for any of the other fields it should be reading). Any suggestions?
# Imports...
import time
import Functions
import arcpy
import sys
import traceback
from arcpy import env
from arcpy.sa import *
Input = "C:/AF_WorkingFiles/AgHealth_Iowa/Databases/Shapes/NitrateWells_LandCover/Test_10Locs.shp"
arcpy.MakeFeatureLayer_management(Input,"InputLyr","","","FID; Shape; OBJECTID; UnqRecordI; UnqWellLoc; Flg90_500")
# List land-cover rasters...
L1 = ["LC85","LC90","LC92","LC00","LC02","LC06","LC09"]
for item in L1:
# Environments...
env.workspace = "C:/AF_WorkingFiles/AgHealth_Iowa/Databases/TempData.mdb"
env.overwriteOutput = "True"
env.snapRaster = "'C:/AF_WorkingFiles/AgHealth_Iowa/Databases/AHS_SpatialData.mdb/" + item + "'"
env.cellSize = "'C:/AF_WorkingFiles/AgHealth_Iowa/Databases/AHS_SpatialData.mdb/" + item + "'"
env.extent = "'C:/AF_WorkingFiles/AgHealth_Iowa/Databases/AHS_SpatialData.mdb/" + item + "'"
env.maintainSpatialIndex = False
env.scratchWorkspace = "C:/AF_WorkingFiles/AgHealth_Iowa/Databases/Scratch"
LC_Raster = item
# List buffer sizes...
L2 = ["500","1000","2000","10000"]
for item in L2:
# For each raster/buffer combo, do the following:
BuffSize = item
Dist = item + " Meters"
a = 0
cur1 = arcpy.UpdateCursor("InputLyr")
rowCount = arcpy.GetCount_management("InputLyr")
print "Row Count =",rowCount
desc = arcpy.Describe("InputLyr")
FieldInfo = desc.fieldInfo
print "Field Count =",FieldInfo.count
print "Field 1 =",FieldInfo.getFieldName(0)
print "Field 2 =",FieldInfo.getFieldName(1)
print "Field 3 =",FieldInfo.getFieldName(2)
print "Field 4 =",FieldInfo.getFieldName(3)
print "Field 5 =",FieldInfo.getFieldName(4)
print "Field 6 =",FieldInfo.getFieldName(5)
for row in cur1:
a = a+1
print row.FID
del cur1, row
... View more
08-08-2012
06:23 AM
|
0
|
0
|
633
|
POST
|
Awesome! Thanks so much! I'm interested in limiting the fields in order to make my code run more efficiently. Is this better to do in MakeFeatureLayer versus UpdateCursor? Also, if I decided to set fields in the UpdateCursor and not in the MakeFeatureLayer, is there a way to check that it's only using/reading the fields I've specified? For instance, the fieldInfo.count reads fields in a Layer, but is there a similar feature within a cursor so I can make sure that it's only reading my specified fields?
... View more
08-07-2012
12:13 PM
|
0
|
0
|
633
|
POST
|
Thanks Jake! Two additional questions, if you don't mind... 1) Is there a way to check the number of fields the cursor is using? When I use the following code after setting the cursor, it returns the entire feature class field count instead of being limited to only the fields I listed in the cursor: desc = arcpy.Describe("InputLyr")
FieldInfo = desc.fieldInfo
print "Field Count =",FieldInfo.count 2) Is there a way to include a changing variable in the fields list for the cursor? For example, cur1 = arcpy.UpdateCursor("InputLyr","","","UnqRecordI; UnqWellLoc; ChangingVariableHere") For instance, would this be correct(?): a = 1
ChangingField = "CmpltFlg_" + str(a)
cur1 = arcpy.UpdateCursor("InputLyr","","","UnqRecordI; UnqWellLoc;",ChangingField,"") Thanks for your help!!
... View more
08-07-2012
07:12 AM
|
0
|
0
|
633
|
POST
|
Hi, I'm using Python code in PyScripter with ArcGIS 10.0 to step through rows using an UpdateCursor, however I cannot get it to skip over optional parameters. UpdateCursor takes: UpdateCursor (dataset, {where_clause}, {spatial_reference}, {fields}, {sort_fields}). I am interested in setting the 'dataset' (required) and 'fields' (optional) only, and I have tried the following three methods as suggested on the ArcGIS 10.0 help site: cur1 = arcpy.UpdateCursor("InputLyr","","",["UnqRecordI","UnqWellLoc"])
cur1 = arcpy.UpdateCursor("InputLyr","#","#",["UnqRecordI","UnqWellLoc"])
cur1 = arcpy.UpdateCursor("InputLyr",fields=["UnqRecordI","UnqWellLoc"]) All three result in errors: "Object: Error in parsing arguments for UpdateCursor" for the first two, and "UpdateCursor() got an unexpected keyword argument 'fields'" for the third. I cannot for the life of me figure out what could be wrong with the parsing. Your help/suggestions would be greatly appreciated! Thanks in advance!
... View more
08-06-2012
01:38 PM
|
0
|
9
|
1461
|
POST
|
Hey Chris, Thanks for your input (and the links)! Most of your comments have already been addressed in previous posts, for instance I've tried using "COMPLETELY_WITHIN" (buffers) instead of "WITHIN_A_DISTANCE" (of points), all data are in the same projected coord. system (and stored in a personal GDB), and the "correct" tool seems to be SelectByLocation (based on its return of all points visually within the buffer feature). (I realize "correct" tool is somewhat arbitrary since you cannot compare each tool to a "true" value, and instead can only compare them to each other. My "correct" tool assignment was based on visual inclusion as viewed in ArcMap, where the data frame coord. system is the same as that of my data.) I think you're right in that there could definitely be something going on with differences in geodesic dists vs. euclidean dists. Is possible for this difference to be >80 meters though? For example, I've found cases where SelectLayerByLocation included a feature that was >80m beyond the perimeter of the buffer feature (as measured in ArcMap with the distance tool)... Abby
... View more
05-18-2012
10:53 AM
|
0
|
0
|
441
|
POST
|
Duncan: Thanks! I have tried creating buffers for a subset of the data and had the same problem. The frustrating thing is, it's not really an issue of choosing the 'better' tool; I need to automate this process, and Python only supports the SelectLayerByLocation tool. When I run each of the two tools by hand in ArcMap, I get the same problem that's reflected when I run SelectLayerByLocation in code, so I think it may be true that these two tools just work differently. SelectByLocation seems to produce the 'true' results, however I cannot implement SelectByLocation by hand for thousands of records and manually imput the count from each selection into a table - it really has to be automated, which leaves me no choice but to go with SelectLayerByLocation 😕 I wonder why the two tools work so differently... and I wonder why the tool worked for you. Maybe I'll test out a polyline dataset since that's what you used. Anyhow, many thanks, and happy Friday!! Abby
... View more
05-18-2012
08:52 AM
|
0
|
0
|
441
|
POST
|
Oh, and my geoprocessing environments are good/cleared of any unwanted settings...
... View more
05-17-2012
08:39 AM
|
0
|
0
|
441
|
POST
|
Duncan, Huh, strange you received the same counts with both methods. Did you use the two "different" tools (i.e., SelectByLocation and SelectLayerByLocation)? How large was your spoof dataset? About 2/3 of my data is in agreement between the two tools, and when disagreement occurs it's seemingly random. Could there be a chance your dataset was too small?? The storage format is the same for all the data I'm using: a personal geodatabase created in ArcCatalog 10.0 - so everything's a feature class. Regarding spatial index and geometries, I've actually "re-created" my well points dataset a couple of times when my clients have made changes, including one time recently, and the count problem still occurred. The ancillary data was downloaded from a state website. I would love to upload the data for your testing, but unfortunately it's private and can't be shared. I do appreciate you being willing to look at it though! The link you provided (re: people w/ a similar problem) didn't connect to the appropriate page; it took me to a Developer Summit login page for voting on something. It's definitely relieving to know others have experienced similar issues though! Maybe try sending the link again? And on a random, new-forum-user note... is there a way to get email notification when someone replies to your thread? I'm "subscribed" to this thread but never receive alerts on when posts occur. Again, many thanks! Abby
... View more
05-17-2012
08:36 AM
|
0
|
0
|
441
|
POST
|
Duncan: I tried running the code on a small data subset using "new selection," and I received the same results 😕 I've noticed, however, that running "Select Layer By Location" in the ArcMap Python window and running "Select By Location" manually from the Selection drop-down menu create the count differences I'm seeing. So, I think there is some underlying difference between these two tools regarding how they draw circles. For a single feature I tested, "Select Layer By Location" highlighted two features up to 80 meters outside of my specified distance, and "Select By Location" did not highlight these features. I plan to inquire about it at the ESRI UC this summer... maybe someone there will have an answer?! Thanks for your persistence with trying to help me figure this one out!
... View more
05-10-2012
01:05 PM
|
0
|
0
|
441
|
POST
|
Duncan: Re. comment 1: The ArcGIS10 Resource Center gives an example without units (http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000072000000 - see bottom of page), and I read somewhere (maybe on that page??) that if no unit is given it uses the linear unit of the input coordinate system (meters in my case). All data I'm using is in the same projection, so everything is in meters. Re. comment 2: For example, lets says you are in ArcMap and you select records from "Schools" where "MascotColor" = red. If you would then like to find "Streets" within 10km of the selected "red" features, you would have to use "add to current selection." Even though the streets data is in a different dataset, Arc sees all selections - i.e., from all datasets - and therefore you have to add to your current selection (which is within "Schools") even though you are making a selection from "Streets." If you use "new selection" it will clear all selections in all datasets to make the new selection, and it will no longer be based on "red" records only (i.e., it will be based on all mascot colors). Therefore, to keep my selection of the single row/feature, I have to "add to" when selecting ancillary data. Does this make any sense? The odd bit is that my code works part of the time, and therefore is incorrect only part of the time. It just seems strange to me that it would work sometimes and not others. I feel like if it were something in the way I use a tool, none of my results would be right. Thoughts?
... View more
05-09-2012
10:31 AM
|
0
|
0
|
441
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|