update

Author: isaacullah

test_pit_sim.py

@@ -9,12 +9,13 @@
 #usage: test_pit_sim.py [-h] [--iters [integer]] [--bounds [integer]] [--padding [integer]] [--sampres [integer]] [--sampint [integer]] [--repeats [integer]] [--sampdist [map name]]
 parser = argparse.ArgumentParser(description='This model simulates additive sampling strategies on a known distribution of artifacts.')
 parser.add_argument('--iters', metavar='integer', type=int, nargs='?', const=3, default=3, help='The number of steps in the additive sampling routine (per iteration)')
-parser.add_argument('--bounds', metavar='integer', type=int, nargs='?', const=300, default=300, help='The size of the sampling universe (a square)')
-parser.add_argument('--padding', metavar='integer', type=int, nargs='?', const=100, default=100, help='The number of additional units to pad the edges of the sampling universe (to allow for grid rotation, etc)')
+parser.add_argument('--bounds', metavar='integer', type=int, nargs='?', const=450, default=450, help='The size of the sampling universe (a square)')
+parser.add_argument('--padding', metavar='integer', type=int, nargs='?', const=50, default=50, help='The number of additional units to pad the edges of the sampling universe (to allow for grid rotation, etc)')
 parser.add_argument('--sampres', metavar='integer', type=int, nargs='?', const=3, default=3, help='The size of each sampling units (also a square)')
 parser.add_argument('--sampint', metavar='integer', type=int, nargs='?', const=50, default=50, help='The interval of the initial sampling units')
 parser.add_argument('--repeats', metavar='integer', type=int, nargs='?', const=100, default=100, help='The number of times to shuffle the grid and resample the distribution (number of iterations for the simulation)')
-parser.add_argument('--sampdist', metavar='map name', type=str, nargs='?', const="", default="", help='The artifact distribution to be sampled (a kernel density map)')
+parser.add_argument('--sampdist', metavar='map name', type=str, nargs='?', const="", default="", help='The artifact distribution to be sampled (a binary presence/absence map or a density map)')
+parser.add_argument('--sites', metavar='map name', type=str, nargs='?', const="Site_areas@gridsim", default="Site_areas@gridsim", help='A binary raster map where areas within sites are coded 1, and areas without are coded 0)')
 # Read in the entered values to internal variables
 args = vars(parser.parse_args())
 iters = args['iters']
@@ -26,6 +27,7 @@
 grass.message("Resampling distance is %s squares" % resamp)
 repeats = args['repeats']
 sampdist = args['sampdist']
+sites = args['sites']
 prefix = '%s_%s' % (sampdist,sampint)
 
 def main(iters,bounds,padding,sampres,sampint,repeats,sampdist,prefix):
@@ -38,10 +40,13 @@ def main(iters,bounds,padding,sampres,sampint,repeats,sampdist,prefix):
     grass.mapcalc("${realpresab}=if(${sampdist} > 0, 1, 0)", overwrite = True, quiet = True, sampdist = sampdist, realpresab = realpresab)
     #extract the number of squares that have at least one artifact
     truepres = int(grass.parse_command("r.stats", quiet = True, flags="cn", input=realpresab, separator="=")["1"])
-    #set up statsfile
+    #set up statsfiles
     f = open('%s%s%s_stats.csv' % (os.getcwd(), os.sep, prefix), 'w+')
     f.write("Iteration,Positive Squares,Interpolated Positives,True Positives,Numeric Difference,Percent Difference,R,R2\n")
     f.flush()
+    f2 = open('%s%s%s_sites_stats.csv' % (os.getcwd(), os.sep, prefix), 'w+')
+    f2.write("Iteration,Inside Positives,Inside Negatives,Outside Positives,Outside Negatives\n")
+    f2.flush()
     #initiate the outer loop, which will randomly rotate/move the sampling grid within the sampling universe
     pad = len(str(repeats)) # set the zfill pad to match the number of repeats
     for x in range(repeats):
@@ -84,8 +89,8 @@ def main(iters,bounds,padding,sampres,sampint,repeats,sampdist,prefix):
         grass.run_command("r.fillnulls", quiet=True, overwrite=True, input=outmap, output=dsurf, method="bicubic")
         # turn density surface into presab to compare to real presab map via regression analysis
         outpres = "%s_%s_interp_pres" % (prefix,str(x + 1).zfill(pad))
-        grass.mapcalc("${outpres}=round(${dsurf})", quiet=True, overwrite=True, dsurf=dsurf, outpres=outpres)
-        # pull some stats
+        grass.mapcalc("${outpres}=eval(a=round(${dsurf}), if(a >= 1, 1, 0))", quiet=True, overwrite=True, dsurf=dsurf, outpres=outpres)
+        # pull some general stats
         presab = grass.parse_command("r.stats", quiet = True,flags="cn", input=outmap, separator="=")
         interp_presab = grass.parse_command("r.stats", quiet = True,flags="cn", input=outpres, separator="=")
         try:
@@ -99,8 +104,29 @@ def main(iters,bounds,padding,sampres,sampint,repeats,sampdist,prefix):
         R = grass.parse_command('r.regression.line', flags='g', mapx=realpresab, mapy=outpres)["R"]
         f.write("%s,%s,%s,%s,%s,%s,%s,%s\n" % (x + 1,pres,intpres,truepres,truepres - intpres,(truepres - intpres)/float(truepres),R, float(R) * float(R)))
         f.flush()
+        # pull some site stats
+        sites_stats = grass.read_command("r.stats", quiet = True,flags="cn", input="%s,%s" % (sites,outmap), separator=",")
+        site_negs = 0
+        site_pos = 0
+        out_negs = 0
+        out_pos = 0
+        for lines in sites_stats.split('\n'):
+            line=lines.split(',')
+            if line[0] is '1' and line[1] is '0':
+                site_negs = line[2]
+            if line[0] is '1' and line[1] is '1':
+                site_pos = line[2]
+            if line[0] is '0' and line[1] is '0':
+                out_negs = line[2]
+            if line[0] is '0' and line[1] is '1':
+                out_pos = line[2]
+            else:
+                pass
+        f2.write("%s,%s,%s,%s,%s\n" % (x+1,site_pos,site_negs,out_pos,out_negs))
+        f.flush()
     #clean up
     f.close()
+    f2.close()
     grass.run_command("g.remove", quiet = True, flags = 'f', type = 'raster', pattern = "*%s" % handle)
 
 if __name__ == "__main__":