2

I am trying to run a Monte Carlo simulation on a point feature class. In order to do so, I am trying to take the values from a column of observed data and randomly "shuffle" them between the points over n additional "simulated" columns. In my case, I'd like to add 199 new columns, each with the values of the the observed data randomly redistributed among the points.

My observed data are binary (1 and 0)

I have been able to do this in R previously and bring the table back into ArcGIS, but the tables have become so unwieldly that I am having difficulty bringing them back into ArcGIS, and would like to run the analysis in ArcGIS, rather than moving in and out.

Does anyone have any idea how to do this? I have tried using cursor, but have failed.

  • Is it that the value for the feature goes into one of the 199 columns randomly or the value gets applied to a random feature? I'm not sure what a Monte Carlo simulation is... – Michael Stimson Sep 3 '14 at 22:07
  • Essentially I have a series of points that, in the "observed" column have values of either 1 or 0. What I would like to do is to create a series of new columns where it is completely random whether or not they are 1 or 0, although the total number of points with 1 must be the same (and conversely, the number of points with 0 must be the same) in each column. In this way, it would be taking the values from the observed column, and randomly redistributing them among the points in the "Shuffled" column. I'd then like to do this 199 times. It's called Monte Carlo because it's like shuffling cards! – Xavier Goldie Sep 3 '14 at 22:12
  • Could you do it one column at a time picking a random feature until the number of 1's reaches the threshold? Otherwise it might be easier to do the shuffling in Excel and relate/join the table - hint, dbf files in shapefiles can be opened with Excel but don't save them with Excel. – Michael Stimson Sep 3 '14 at 22:20
  • 1
    If you have the feature class stored in a database I may have a method to do this – MickyT Sep 3 '14 at 22:21
  • I have approximately 700,000 points so it's a bit unwieldly in excel (I have tried endlessly - I can't import the CSV and Excel sheets back into ArcGIS - they keep getting stuck). The way that you mentioned would be good - it would be fine to pick a random feature, update it as 1, and keep doing so until the number of 1s reaches the thresholds, but I don't know how to do that in ArcGIS/ArcPy? – Xavier Goldie Sep 3 '14 at 22:23
5

You can accomplish this task using the following workflow:

Generate a list of all the values (1's and 0's) in a field

vals = [row[0] for row in arcpy.da.SearchCursor(fc, ["binary"])]

Create a loop equal to the number of fields you want to create

for x in range(2):

Within that loop, randomly shuffle the list of 1's and 0's

random.shuffle(vals)

Incorporate the range value into the Field name

fieldName = "Field" + str((x + 1))

Start an UpdateCursor to update the field in range x

with arcpy.da.UpdateCursor(fc, fieldName) as cursor:

Perhaps most important, create a counter and use that to index the shuffled list. Add that indexed value to the attribute table

row[0] = int(vals[count])

import arcpy, os, random

# Your input feature class    
fc = r'C:\path\to\featureclass'

# Get a list of all the values (i.e. 1's and 0's) in the field called "binary"    
vals = [row[0] for row in arcpy.da.SearchCursor(fc, ["binary"])]

# In your case, you probably want to change the range to 199
for x in range(2):
    random.shuffle(vals)   #Randomly shuffle list
    fieldName = "Field" + str((x + 1))

    arcpy.AddField_management(fc, fieldName, "SHORT")
    count = 0
    with arcpy.da.UpdateCursor(fc, fieldName) as cursor:
        for row in cursor:
            row[0] = int(vals[count])
            count = count + 1
            cursor.updateRow(row)

del count
  • 1
    I like that Aaron, it's a bit more versatile than mine and actually uses existing values (doesn't need to be 0 or 1), also more compact. – Michael Stimson Sep 3 '14 at 22:54
2

Here's something that I put together, it woks in mostly the same way as Aarons' script.

import sys, os, arcpy, random

InFC = sys.argv[1]

ThresholdOfOnes = 50 # percent that is!
Columns = 199

ColRange = range(1,Columns + 1) 
InCount = int(arcpy.GetCount_management(InFC).getOutput(0))
ThresholdOfOnes = (InCount * ThresholdOfOnes) / 100 # convert percent into a number

OIDlist = list()
desc = arcpy.Describe(InFC)
OIDfield = desc.OIDFieldName

# get the OIDs into a list
with arcpy.da.SearchCursor(InFC,OIDfield) as srch:
    for row in srch:
        OIDlist.append(row[0])

for ColNum in ColRange:
    FieldName = "Col_%d" % ColNum
    SelOIDlist = list()
    while len(SelOIDlist) < ThresholdOfOnes:
        PosIndex = random.randint(0,len(OIDlist)) # pick a random number
        if OIDlist[PosIndex] not in SelOIDlist:
            # if the OID isn't already in the list
            SelOIDlist.append(OIDlist[PosIndex])
    # now to select with the oidlist
    # start building a definiton query
    DefQ = "%s in (" % OIDfield
    isFirst = True
    for ThisOID in SelOIDlist:
        if isFirst:
            # the first one doesn't get a comma
            DefQ = DefQ + str(ThisOID)
            isFirst = False
        else:
            DefQ = DefQ + "," + str(ThisOID)
    DefQ = DefQ + ")"  # close the bracket

    # add the field if it doesn't exist
    fList = arcpy.ListFields(InFC,FieldName)
    if not fList:
        arcpy.AddField_management(InFC,FieldName,"SHORT")

    # blank the field
    arcpy.CalculateField_management(InFC,FieldName,"0","PYTHON")

    # start a cursor with just the random features
    with arcpy.da.UpdateCursor(InFC,FieldName,DefQ) as UpCur:
        for UpRow in UpCur:
            UpRow[0] = 1           # set the value
            UpCur.updateRow(UpRow) # store it

This will work for shapefiles or feature classes in databases.. it gets a list of the oids and then randomly picks from the list (ensuring no duplication) until the threshold is reached then makes the list into a definition query and updates the features. Along the way if the field doesn't exist it will add it, ensuring that previous values are overwritten so if you need to re-random your feature class it's safe to run again.

  • Thank you so much guys - I will run these to see how they work! I really need to work on my cursor knowledge! – Xavier Goldie Sep 3 '14 at 22:58
  • That's some fast coding since the last one we worked on! Nice solution. – Aaron Sep 3 '14 at 23:00
  • Thank you @Aaron, your method of reading the values, using python to shuffle and writing them back is brilliant! Wish I had thought of that... I will be using that later for tasks like randomly determine which 10% gets intensive QA. – Michael Stimson Sep 3 '14 at 23:13
2

Some fancy dictionary thing and a one-pass update cursor would be the fastest/best, but something like this would work too:

import arcpy, random
pntFC = r"C:\temp\test.gdb\test_pnts"
boolField = "ORIG_TRUTH" #original bool field
falseOidList = [r[0] for r in arcpy.da.SearchCursor(pntFC, ["OID@"], boolField + " = 0")]
trueOidList = [r[0] for r in arcpy.da.SearchCursor(pntFC, ["OID@"], boolField + " = 1")]
falseCount = len(falseOidList)
trueCount = len(trueOidList) 
for i in range(1,199 + 1):
   randomFieldName = "RND_" + str(i)
   arcpy.AddField_management(randomFieldName, "SHORT")
   randomTrueOids = random.sample(falseOidList + trueOidList, trueCount)
   arcpy.MakeFeatureLayer_management(pntFC, "fl", "OID in (" + ",".join(str(oid) for oid in (falseOidList + trueOidList)) + ")"
   arcpy.CalculateField_management("fl", randomFieldName, "1", "PYTHON")
   arcpy.MakeFeatureLayer_management(pntFC, "fl", randomFieldName + " IS NULL")
   arcpy.CalculateField_management("fl", randomFieldName, "0", "PYTHON")
  • Your method of joining the OIDs into a where clause is quite cool.. have a look at my answer regarding obtaining the OIDfieldName to work with geodatabase, SDE and shapefile data. Note, if you're working with shapefiles there wont be any null values as shapefiles don't support NULL - but that's ok, they default to '0'. – Michael Stimson Sep 4 '14 at 0:38
  • Yes - I wrote that on the fly, and the Describe(myFC).oidFieldName would surely be necessary. I actually edited the code so it was using a FGDB FC (so my "OID = ..." code " – csny490 Sep 4 '14 at 18:55
  • Argh - can't delete this comment – csny490 Sep 4 '14 at 23:58

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.