# randomselection.py # 19 May 2011 # kimo@ollivier.co.nz # # Input: parcel polygon featureclass # item to sum (area) # threshold (check min, mean, max values for reasonable figure) # Assumes values always less than threshold or you get a null list # Output: Total sum and layer with selected parcels # set up to allow iteration # uses ArcGIS 10 but could use 9.3 import arcpy import sys import os import random import datetime try: fc = sys.argv[1] item = sys.argv[2] threshold = sys.argv[3] except: # debugging and standalone operation fc = "c:/project/forum/mobile.gdb/parcel" item = "area" threshold = 100000.0 # square metres def getdict(fc,item): ''' Make a dictionary once for the sampler ''' dParcel = {} for row in arcpy.SearchCursor(fc,"featcode = 'parcel'"): if item.lower() == 'area': dParcel[row.objectid] = row.shape.area else : dParcel[row.objectid] = row.GetValue(item) return dParcel def getsample(dParcel,threshold): '''Given a dict of values return a random selection less than the threshold. Shuffles the list and pops off the ids until threshold reached ie sampling without replacement ''' total = 0 lstObj = [] lstKey = dParcel.keys() random.shuffle(lstKey) while total < threshold: nObj = lstKey.pop() nVal = dParcel[nObj] total = nVal if total > threshold: break # do not exceed threshold else : lstObj.append(nObj) total += nVal return (total,lstObj) # --------------- main --------------- print start = datetime.datetime.now() print start # # make a dictionary dParcel = getdict(fc,item) print "records",len(dParcel) # sample the dictionary total,lstObj = getsample(dParcel,threshold) print "Total: ",total print "No of parcels: ",len(lstObj) print "objectid: value" for ob in lstObj: print ob,":",dParcel[ob],", ", print # make up an sql query to select the features sqlQuery = "OBJECTID in (" + ",".join([str(id) for id in lstObj]) + ")" print sqlQuery arcpy.env.overwriteOutput = True arcpy.MakeFeatureLayer_management(fc,"fc_lay",sqlQuery) lyrName = arcpy.CreateUniqueName("xyz.lyr",os.path.dirname(os.path.dirname(fc))) print lyrName arcpy.SaveToLayerFile_management("fc_lay",lyrName) # Or sample the dictionary 100 times for it in range(100): total,lstObj = getsample(dParcel,threshold) print total,len(lstObj) # finished elapsed = datetime.datetime.now() - start print "\n\nElapsed",elapsed, "Well done!" # positive feedback
2011-05-19 20:40:43.640000 records 7170 Total: 145316.242757 No of parcels: 21 objectid: value 726 : 809.726863177 , 4530 : 809.419546424 , 4752 : 1012.70775115 , 6120 : 809.397052955 , 3022 : 809.418725814 , 253 : 986.240771943 , 4445 : 1132.16238039 , 7161 : 928.75695316 , 5747 : 964.13873029 , 2188 : 864.834770195 , 4164 : 20.9192428737 , 187 : 839.573627519 , 2745 : 811.268061104 , 2572 : 1002.32513709 , 1445 : 1381.14243957 , 2674 : 812.661986759 , 6925 : 25126.7099296 , 5479 : 1395.56544711 , 4791 : 810.812432473 , 5406 : 848.317697436 , 6543 : 72658.1213786 , OBJECTID in (726,4530,4752,6120,3022,253,4445,7161,5747,2188,4164,187,2745,2572,1445,2674,6925,5479,4791,5406,6543) c:/project/forum\xx0.lyr 185921.392647 4 145871.629063 27 154201.63183 1 ... 140524.256299 0 ... 190822.946912 9 193668.483773 23 334292.358346 16 141967.626185 23 15109126.1898 56 ... 413536.709991 2 157130.891363 6 Elapsed 0:00:02.917000 Well done!
Angemeldete Mitglieder können Beiträge verfassen, Updates folgen und mehr. Neu hier? Registriere ein kostenloses Konto.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.