-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathr.agropast.semiadaptive.py
executable file
·1113 lines (1098 loc) · 65.6 KB
/
r.agropast.semiadaptive.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/python
#
############################################################################
#
# MODULE: r.agropast.semiadaptive.tenure.py
# AUTHOR(S): Isaac Ullah, University of Pittsburgh, Arizona State University
# PURPOSE: Simulates agricultural and pastoral landuse and tracks yields and environmental
# impacts. Farming and grazing strategies and yield goals are predetermined by the
# researcher, and do not change (adapt) during the simulation. However, catchment sizes
# can be adapted over time to meet these goals. This version implments a land tenuring alogrithm.
# Requires r.landscape.evol.
# ACKNOWLEDGEMENTS: National Science Foundation Grant #BCS0410269, Center for Comparative Archaeology at
# the University of Pittsburgh, Center for Social Dynamics and Complexity at Arizona
# State University
# COPYRIGHT: (C) 2014 by Isaac Ullah, University of Pittsburgh
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
#############################################################################
#%Module
#% description: Simulates agricultural and pastoral landuse and tracks yields and environmental impacts. Basic farming and grazing strategies and yield goals are predetermined by the researcher, and do not change (adapt) during the simulation. However, catchment sizes can be adapted over time to meet these goals. This version implments a land tenuring alogrithm. Requires r.landscape.evol. Note that some stats files will be written to current mapset, and will be appended to if you run the simulation again with the same prefix.
#%END
##################################
#Simulation Control
##################################
#%option
#% key: years
#% type: integer
#% description: Number of iterations ("years") to run
#% answer: 500
#% required: yes
#% guisection: Simulation Control
#%END
#%option
#% key: prfx
#% type: string
#% description: Prefix for all output maps
#% answer: p30f10sim
#% required: yes
#% guisection: Simulation Control
#%END
##################################
#Agent Properties
##################################
#%option
#% key: numpeople
#% type: double
#% description: Number of people in the village(s) (stays constant throughout the simulation)
#% answer: 30
#% guisection: Agent Properties
#%END
#%option
#% key: agentmem
#% type: integer
#% description: Length of the "memory" of the agent in years. The agent will use the mean surplus/defict information from this many of the most recent previous years when making a subsistence plan for the current year.
#% answer: 5
#% guisection: Agent Properties
#%END
#%option
#% key: aglabor
#% type: double
#% description: The amount of agricultural labor an average person of the village can do in a year (in "person-days")
#% answer: 300
#% guisection: Agent Properties
#%END
#%option
#% key: cerealreq
#% type: double
#% description: Amount of cereals that would be required per person per year if cereals were the ONLY food item (Kg)
#% answer: 370
#% guisection: Agent Properties
#%END
#%option
#% key: animals
#% type: double
#% description: Number of herd animals that would be needed per person per year if pastoral products were the ONLY food item
#% answer: 60
#% guisection: Agent Properties
#%END
#%option
#% key: fodderreq
#% type: double
#% description: Amount of fodder required per herd animal per year (Kg)
#% answer: 680
#% guisection: Agent Properties
#%END
#%option
#% key: a_p_ratio
#% type: double
#% description: Actual ratio of agricultural to pastoral foods in the diet (0-1, where 0 = 100% agricultural and 1 = 100% pastoral)
#% answer: 0.20
#% options: 0.0-1.0
#% guisection: Agent Properties
#%END
#%option
#% key: costsurf
#% type: string
#% gisprompt: old,cell,raster
#% description: Map of movement costs from the center of the agricultural/grazing catchments (from r.walk or r.cost).
#% guisection: Agent Properties
#%END
##################################
#Farming Options
##################################
#%option
#% key: agcatch
#% type: string
#% gisprompt: old,cell,raster
#% description: Map of the largest possible agricultural catchment map (From r.catchment or r.buffer, where catchment is a single integer value, and all else is NULL)
#% guisection: Farming
#%END
#%option
#% key: agmix
#% type: double
#% description: The wheat/barley ratio (e.g., 0.0 for all wheat, 1.0 for all barley, 0.5 for an equal mix)
#% answer: 0.25
#% options: 0.0-1.0
#% guisection: Farming
#%END
#%option
#% key: nsfieldsize
#% type: double
#% description: North-South dimension of fields in map units (Large field sizes may lead to some overshoot of catchment boundaries)
#% answer: 20
#% guisection: Farming
#%END
#%option
#% key: ewfieldsize
#% type: double
#% description: East-West dimension of fields in map units (Large field sizes may lead to some overshoot of catchment boundaries)
#% answer: 50
#% guisection: Farming
#%END
#%option
#% key: fieldlabor
#% type: double
#% description: Number of person-days required to till, sow, weed, and harvest one farm field in a year.
#% answer: 50
#% guisection: Farming
#%END
#%option
#% key: farmval
#% type: double
#% description: The landcover value for farmed fields (Should correspond to an appropriate value from your landcover rules file.)
#% answer: 5
#% guisection: Farming
#%END
#%option
#% key: farmimpact
#% type: double
#% description: The mean and standard deviation of the amount to which farming a patch decreases its fertility (in percentage points of maximum fertility, entered as: "mean,std_dev"). Fertility impact values of indvidual farm plots will be randomly chosen from a gaussian distribution that has this mean and std_dev.
#% answer: 3,2
#% multiple: yes
#% guisection: Farming
#%END
#%option
#% key: maxwheat
#% type: double
#% description: Maximum amount of wheat that can be grown (kg/ha)
#% answer: 3500
#% guisection: Farming
#%END
#%option
#% key: maxbarley
#% type: double
#% description: Maximum amount of barley that can be grown (kg/ha)
#% answer: 2500
#% guisection: Farming
#%END
#%option
#% key: tenuredrop
#% type: double
#% description: Threshold for dropping land out of Tenure (with flag -M). If the value is positive it is interpreted as a percentage below the yearly average yield of all farm cells. If the value is negative, it is interpreted as a percentage of the maximum yield of all farm cells. Value will always be interpreted as a proportion between 0 and 1.
#% answer: 0.0
#% options: -1.0-1.0
#% guisection: Farming
#%END
#%flag
#% key: i
#% description: -i Implement Land Tenure with a satisficing strategy (land is never dropped, only added if needed)
#% guisection: Farming
#%end
#%flag
#% key: m
#% description: -m Implement Land Tenure with a maximizing strategy (land is dropped if below the threshold defined by "tenuredrop")
#% guisection: Farming
#%end
##################################
#Grazing Options
##################################
#%option
#% key: grazecatch
#% type: string
#% gisprompt: old,cell,raster
#% description: Map of the largest possible grazing catchment (From r.catchment or r.buffer, where catchment is a single integer value, and all else is NULL).
#% guisection: Grazing
#%END
#%option
#% key: mingraze
#% type: double
#% description: Minimum amount of vegetation on a cell for it to be considered grazeable (in value of succession stages, matching landcover rules file)
#% answer: 2
#% guisection: Grazing
#%END
#%option
#% key: grazespatial
#% type: double
#% description: Spatial dependency of the grazing pattern in map units. This value determines how "clumped" grazing patches will be. A value close to 0 will produce a perfectly randomized grazing pattern with patch size=resolution, and larger values will produce increasingly clumped grazing patterns, with the size of the patches corresponding to the value of grazespatial (in map units).
#% answer: 50
#% guisection: Grazing
#%END
#%option
#% key: grazepatchy
#% type: double
#% description: Coefficient that determines the patchiness of the grazing catchemnt. Value must be non-zero, and usually will be <= 1.0. Values close to 0 will create a patchy grazing pattern, values close to 1 will create a "smooth" grazing pattern. Actual grazing patches will be sized to the resolution of the input landcover map.
#% answer: 1.0
#% guisection: Grazing
#%END
#%option
#% key: maxgrazeimpact
#% type: integer
#% description: Maximum impact of grazing in units of "landcover.succession". Grazing impact values of individual patches will be chosen from a gaussian distribution between 1 and this maximum value (i.e., most values will be between 1 and this max). Value must be >= 1.
#% answer: 3
#% guisection: Grazing
#%END
#%option
#% key: manurerate
#% type: double
#% description: Base rate that animal dung contributes to fertility increase on a grazed patch in units of percentage of maximum fertility regained per increment of grazing impact. Actual fertility regain values are thus calculated as "manure_rate x grazing_impact", so be aware that this variable interacts with the grazing impact settings you have chosen.
#% answer: 0.2
#% options: 0.0-100.0
#% guisection: Grazing
#%END
#%option
#% key: fodder_rules
#% type: string
#% gisprompt: string
#% description: Path to foddering rules file for calculation of fodder gained by grazing
#% answer: /home/mdlpd/Dropbox/Scripts_Working_Dir/rules/fodder_rules.txt
#% guisection: Grazing
#%END
#%flag
#% key: f
#% description: -f Do not graze in unused portions of the agricultural catchment (i.e., do not graze on "fallowed" fields, and thus no "manuring" of those fields will occur)
#% guisection: Grazing
#%end
#%flag
#% key: g
#% description: -g Do not allow "stubble grazing" on harvested fields (and thus no "manuring" of fields)
#% guisection: Grazing
#%end
##################################
#Landcover Dynamics Options
##################################
#%option
#% key: fertilrate
#% type: double
#% description: Comma separated list of the mean and standard deviation of the natural fertility recovery rate (percentage by which soil fertility increases per year if not farmed, entered as: "mean,stdev".) Fertility recovery values of individual landscape patches will be randomly chosen from a gaussian distribution that has this mean and std_dev.
#% answer: 2,0.5
#% multiple: yes
#% guisection: Landcover Dynamics
#%END
#%option
#% key: inlcov
#% type: string
#% gisprompt: old,cell,raster
#% description: Input landcover map (Coded according to landcover rules file below)
#% guisection: Landcover Dynamics
#%END
#%option
#% key: infert
#% type: string
#% gisprompt: old,cell,raster
#% description: Input fertility map (Coded according to percentage retained fertility, with 100 being full fertility)
#% guisection: Landcover Dynamics
#%END
#%option
#% key: maxlcov
#% type: string
#% gisprompt: old,cell,raster
#% description: Maximum value of landcover for the landscape (Can be a single numerical value or a cover map of different maximum values. Shouldn't exceed maximum value in rules file')
#% answer: 50
#% guisection: Landcover Dynamics
#%END
#%option
#% key: maxfert
#% type: string
#% gisprompt: old,cell,raster
#% description: Maximum value of fertility for the landscape (Can be a single numerical value or a cover map of different maximum values. Shouldn't exceed 100)
#% answer: 100
#% guisection: Landcover Dynamics
#%END
#%option
#% key: lc_rules
#% type: string
#% gisprompt: string
#% description: Path to reclass rules file for landcover map
#% answer: /home/mdlpd/Dropbox/Scripts_Working_Dir/rules/luse_reclass_rules.txt
#% guisection: Landcover Dynamics
#%END
#%option
#% key: cfact_rules
#% type: string
#% gisprompt: string
#% description: Path to recode rules file for c-factor map
#% answer: /home/mdlpd/Dropbox/Scripts_Working_Dir/rules/cfactor_recode_rules.txt
#% guisection: Landcover Dynamics
#%END
#%flag
#% key: c
#% description: -c Keep C-factor (and rainfall excess) maps for each year
#% guisection: Landcover Dynamics
#%end
##################################
#Landscape Evolution Options
##################################
#%option
#% key: elev
#% type: string
#% gisprompt: old,cell,raster
#% description: Input elevation map (DEM of surface)
#% guisection: Landscape Evolution
#%end
#%option
#% key: initbdrk
#% type: string
#% gisprompt: old,cell,raster
#% description: Bedrock elevations map (DEM of bedrock)
#% answer:
#% guisection: Landscape Evolution
#%end
#%option
#% key: k
#% type: double
#% gisprompt: old,cell,raster
#% description: Soil erodability index (K factor) map or constant
#% answer: 0.42
#% guisection: Landscape Evolution
#%end
#%option
#% key: sdensity
#% type: double
#% gisprompt: old,cell,raster
#% description: Soil density map or constant [T/m3] for conversion from mass to volume
#% answer: 1.2184
#% guisection: Landscape Evolution
#%end
#%option
#% key: kt
#% type: double
#% description: Stream transport efficiency variable (0.001 for a soft substrate, 0.0001 for a normal substrate, 0.00001 for a hard substrate, 0.000001 for a very hard substrate)
#% answer: 0.0001
#% options: 0.001,0.0001,0.00001,0.000001
#% guisection: Landscape Evolution
#%end
#%option
#% key: loadexp
#% type: double
#% description: Stream transport type variable (1.5 for mainly bedload transport, 2.5 for mainly suspended load transport)
#% answer: 1.5
#% options: 1.5,2.5
#% guisection: Landscape Evolution
#%end
#%option
#% key: manningn
#% type: string
#% gisprompt: old,cell,raster
#% description: Map or constant of the value of Manning's "N" value for channelized flow.
#% answer: 0.05
#% required : no
#% guisection: Landscape Evolution
#%end
#%option
#% key: kappa
#% type: double
#% description: Hillslope diffusion (Kappa) rate map or constant [m/kyr]
#% answer: 1
#% guisection: Landscape Evolution
#%end
#%option
#% key: rain
#% type: string
#% gisprompt: old,cell,raster
#% description: Precip totals for the average storm [mm] (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
#% answer: 20.61
#% guisection: Landscape Evolution
#%end
#%option
#% key: r
#% type: string
#% description: Rainfall (R factor) constant (AVERAGE FOR WHOLE MAP AREA) (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
#% answer: 4.54
#% guisection: Landscape Evolution
#%end
#%option
#% key: storms
#% type: string
#% description: Number of storms per year (integer) (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
#% answer: 25
#% guisection: Landscape Evolution
#%end
#%option
#% key: stormlength
#% type: string
#% description: Average length of the storm [h] (or path to climate file of comma separated values of "rain,R,storms,stormlength", with a new line for each year of the simulation)
#% answer: 24.0
#% guisection: Landscape Evolution
#%end
#%option
#% key: speed
#% type: double
#% description: Average velocity of flowing water in the drainage [m/s]
#% answer: 1.4
#% guisection: Landscape Evolution
#%end
#%option
#% key: cutoff1
#% type: double
#% description: Flow accumulation breakpoint value for shift from diffusion to overland flow
#% answer: 0
#% guisection: Landscape Evolution
#%end
#%option
#% key: cutoff2
#% type: double
#% description: Flow accumulation breakpoint value for shift from overland flow to rill/gully flow (if value is the same as cutoff1, no sheetwash procesess will be modeled)
#% answer: 100
#% guisection: Landscape Evolution
#%end
#%option
#% key: cutoff3
#% type: double
#% description: Flow accumulation breakpoint value for shift from rill/gully flow to stream flow (if value is the same as cutoff2, no rill procesess will be modeled)
#% answer: 100
#% guisection: Landscape Evolution
#%end
#%option
#% key: smoothing
#% type: string
#% description: Amount of additional smoothing (answer "no" unless you notice large spikes in the erdep rate map)
#% answer: no
#% options: no,low,high
#% guisection: Landscape Evolution
#%end
#%flag
#% key: 1
#% description: -1 Calculate streams as 1D difference instead of 2D divergence
#% guisection: Landscape Evolution
#%end
#%flag
#% key: k
#% description: -k Keep ALL temporary maps (overides flags -drst). This will make A LOT of maps!
#% guisection: Landscape Evolution
#%end
#%flag
#% key: d
#% description: -d Don't output yearly soil depth maps
#% guisection: Landscape Evolution
#%end
#%flag
#% key: r
#% description: -r Don't output yearly maps of the erosion/deposition rates ("ED_rate" map, in vertical meters)
#% guisection: Landscape Evolution
#%end
#%flag
#% key: s
#% description: -s Keep all slope maps
#% guisection: Landscape Evolution
#%end
#%flag
#% key: t
#% description: -t Keep yearly maps of the Transport Capacity at each cell ("Qs" maps)
#% guisection: Landscape Evolution
#%end
#%flag
#% key: e
#% description: -e Keep yearly maps of the Excess Transport Capacity (divergence) at each cell ("DeltaQs" maps)
#% guisection: Landscape Evolution
#%end
import sys
import os
import tempfile
import random
import numpy
import grass.script as grass
#main block of code starts here
def main():
grass.message("Setting up Simulation........")
#setting up Land Use variables for use later on
agcatch = options['agcatch']
nsfieldsize = options['nsfieldsize']
ewfieldsize = options['ewfieldsize']
grazecatch = options['grazecatch']
grazespatial = options['grazespatial']
grazepatchy = options['grazepatchy']
maxgrazeimpact = options['maxgrazeimpact']
manurerate = options['manurerate']
inlcov = options['inlcov']
years = int(options['years'])
farmval = options['farmval']
maxlcov = options['maxlcov']
prfx = options['prfx']
lc_rules = options['lc_rules']
cfact_rules = options['cfact_rules']
fodder_rules = options['fodder_rules']
infert = options['infert']
maxfert = options['maxfert']
maxwheat = options['maxwheat']
maxbarley = options['maxbarley']
mingraze = options['mingraze']
tenuredrop = options['tenuredrop']
costsurf = options['costsurf']
agmix = options['agmix']
numpeople = float(options['numpeople'])
agentmem = int(options['agentmem'])
animals = float(options['animals'])
agratio = 1 - float(options['a_p_ratio'])
pratio = float(options['a_p_ratio'])
indcerreq = float(options['cerealreq'])
indfodreq = float(options['fodderreq'])
aglabor = float(options['aglabor'])
fieldlabor = float(options['fieldlabor'])
cereal_pers = numpeople * agratio
fodder_anim = animals * numpeople * pratio
cerealreq = indcerreq * cereal_pers
fodderreq = indfodreq * fodder_anim
totlabor = numpeople * aglabor
maxfields = int(round(totlabor / fieldlabor))
#these are various rates with min and max values entered in the gui that we need to parse
farmimpact = map(float, options['farmimpact'].split(','))
fertilrate = map(float, options['fertilrate'].split(','))
#Setting up Landscape Evol variables to write the r.landscape.evol command later
elev = options["elev"]
initbdrk = options["initbdrk"]
k = options["k"]
sdensity = options["sdensity"]
kappa = options["kappa"]
manningn = options["manningn"]
cutoff1 = options["cutoff1"]
cutoff2 = options["cutoff2"]
cutoff3 = options["cutoff3"]
speed = options["speed"]
kt = options["kt"]
loadexp = options["loadexp"]
smoothing = options["smoothing"]
#these values could be read in from a climate file, so check that, and act accordingly
rain2 = []
try:
rain1 = float(options['rain'])
for year in range(int(years)):
rain2.append(rain1)
except:
with open(options['rain'], 'rU') as f:
for line in f:
rain2.append(line.split(",")[0])
#check for text header and remove if present
try:
float(rain2[0])
except:
del rain2[0]
#throw a fatal error if there aren't enough values in the column
if len(rain2) != int(years):
grass.fatal("Number of rows of rainfall data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
sys.exit(1)
R2 = []
try:
R1 = float(options['r'])
for year in range(int(years)):
R2.append(R1)
except:
with open(options['r'], 'rU') as f:
for line in f:
R2.append(line.split(",")[1])
#check for text header and remove if present
try:
float(R2[0])
except:
del R2[0]
#throw a fatal error if there aren't enough values in the column
if len(R2) != int(years):
grass.fatal("Number of rows of R-Factor data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
sys.exit(1)
storms2 = []
try:
storms1 = float(options['storms'])
for year in range(int(years)):
storms2.append(storms1)
except:
with open(options['storms'], 'rU') as f:
for line in f:
storms2.append(line.split(",")[2])
#check for text header and remove if present
try:
float(storms2[0])
except:
del storms2[0]
#throw a fatal error if there aren't enough values in the column
if len(storms2) != int(years):
grass.fatal("Number of rows of storm frequency data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
sys.exit(1)
stormlength2 = []
try:
stormlength1 = float(options['stormlength'])
for year in range(int(years)):
stormlength2.append(stormlength1)
except:
with open(options['stormlength'], 'rU') as f:
for line in f:
stormlength2.append(line.split(",")[3])
#check for text header and remove if present
try:
float(stormlength2[0])
except:
del stormlength2[0]
#throw a fatal error if there aren't enough values in the column
if len(stormlength2) != int(years):
grass.fatal("Number of rows of storm length data in your climate file\n do not match the number of iterations you wish to run.\n Please ensure that these numbers match and try again")
sys.exit(1)
#get the process id to tag any temporary maps we make for easy clean up in the loop
pid = os.getpid()
#we need to separate out flags used by this script, and those meant to be sent to r.landscape.evol. We will do this by popping them out of the default "flags" dictionary, and making a new dictionary called "use_flags"
use_flags = {}
use_flags.update({'g': flags.pop('g'), 'f': flags.pop('f'), 'c': flags.pop('c'), 'i': flags.pop('i'), 'm': flags.pop('m')})
#check the two tenure flags, and throw a fatal error if both are set.
if use_flags['i'] is True and use_flags['m'] is True:
grass.fatal("You have set both the i and the m flags, and this is not allowed. Please uncheck one and try again.")
#now assemble the flag string for r.landscape.evol'
levol_flags = []
for flag in flags:
if flags[flag] is True:
levol_flags.append(flag)
#check if maxlcov is a map or a number, and grab the actual max value for the stats file
try:
maxval = int(float(maxlcov))
except:
maxlcovdict = grass.parse_command('r.univar', flags = 'ge', map = maxlcov)
maxval = int(float(maxlcovdict['max']))
#check if maxfert is a map or a number, and grab the actual max value for the stats file
try:
maxfertval = int(float(maxfert))
except:
maxfertdict = grass.parse_command('r.univar', flags = 'ge', map = maxfert)
maxfertval = int(float(maxfertdict['max']))
#set up the stats files names
env = grass.gisenv()
statsdir = os.path.join(env['GISDBASE'], env['LOCATION_NAME'], env['MAPSET'])
textout = statsdir + os.sep + prfx + '_landcover_temporal_matrix.txt'
textout2 = statsdir + os.sep + prfx + '_fertility_temporal_matrix.txt'
textout3 = statsdir + os.sep + prfx + '_yields_stats.txt'
textout4 = statsdir + os.sep + prfx + '_landcover_and_fertility_stats.txt'
statsout = statsdir + os.sep + prfx + '_erdep_stats.txt'
# Make color rules for landcover, cfactor, and soil fertilty maps
lccolors = tempfile.NamedTemporaryFile()
lccolors.write('0 grey\n10 red\n20 orange\n30 brown\n40 yellow\n%s green' % maxval)
lccolors.flush()
cfcolors = tempfile.NamedTemporaryFile()
cfcolors.write('0.1 grey\n0.05 red\n0.03 orange\n0.01 brown\n0.008 yellow\n0.005 green')
cfcolors.flush()
fertcolors = tempfile.NamedTemporaryFile()
fertcolors.write('0 white\n20 grey\n40 yellow\n60 orange\n80 brown\n100 black')
fertcolors.flush()
#Figure out the number of cells per hectare and how many square meters per cell to use as conversion factors for yields
region = grass.region()
cellperhectare = 10000 / (float(region['nsres']) * float(region['ewres']))
#sqmeterpercell = (float(region['nsres']) * float(region['ewres']))
#do same for farm field size
fieldsperhectare = 10000 / (float(nsfieldsize) * float(ewfieldsize))
#find conversion from field size to cell size
#cellsperfield = (float(nsfieldsize) * float(ewfieldsize)) / (float(region['nsres']) * float(region['ewres']))
#find out the maximum amount of cereals per field according to the proper mix
#maxyield = (((1-float(agmix))*float(maxwheat))+(float(agmix)*float(maxbarley)))/fieldsperhectare
#find out number of digits in 'years' for zero padding
digits = len(str(abs(years)))
#set up the agent memory
farmingmemory = []
farmyieldmemory = []
grazingmemory = []
grazeyieldmemory = []
grass.message('Simulation will run for %s iterations.\n\n............................STARTING SIMULATION...............................' % years)
#Set up loop
for year in range(int(years)):
now = str(year + 1).zfill(digits)
then = str(year).zfill(digits)
#grab the current climate vars from the lists
rain = rain2[year]
r = R2[year]
storms = storms2[year]
stormlength = stormlength2[year]
#figure out total precip (in meters) for the year for use in the veg growth and farm yields formulae
precip = 0.001 * (float(rain) * float(storms))
grass.message('_____________________________\nSIMULATION YEAR: %s\n--------------------------' % now)
#make some map names
fields = "%s_Year_%s_Farming_Impacts_Map" % (prfx, now)
outlcov = "%s_Year_%s_Landcover_Map" % (prfx, now)
outfert = "%s_Year_%s_Soil_Fertilty_Map" % (prfx, now)
outcfact = "%s_Year_%s_Cfactor_Map" % (prfx, now)
grazeimpacts = "%s_Year_%s_Gazing_Impacts_Map" % (prfx, now)
outxs = "%s_Year_%s_Rainfall_Excess_Map" % (prfx, now)
#check if this is year one, use the starting landcover and soilfertily and calculate soildepths
if (year + 1) == 1:
oldlcov = inlcov
oldfert = infert
oldsdepth = "%s_Year_%s_Soil_Depth_Map" % (prfx, then)
grass.mapcalc("${sdepth}=(${elev}-${bdrk})", quiet ="True", sdepth = oldsdepth, elev = elev, bdrk = initbdrk)
else:
oldlcov = "%s_Year_%s_Landcover_Map" % (prfx, then)
oldfert = "%s_Year_%s_Soil_Fertilty_Map" % (prfx, then)
oldsdepth = "%s_Year_%s_Soil_Depth_Map" % (prfx, then)
#GENERATE FARM IMPACTS
#create some temp map names
tempfields = "%stemporary_fields_map" % pid
tempimpacta = "%stemporary_farming_fertility_impact" % pid
#temporarily change region resolution to align to farm field size
grass.use_temp_region()
grass.run_command('g.region', quiet = 'True',nsres = nsfieldsize, ewres = ewfieldsize)
#generate the yields
grass.message("Calculating potential farming yields.....")
#Calculate the wheat yield map (kg/cell)
tempwheatreturn = "%stemporary_wheat_yields_map" % pid
grass.mapcalc("${tempwheatreturn}=eval(x=if(${precip} > 0, (0.51*log(${precip}))+1.03, 0), y=if(${sfertil} > 0, (0.28*log(${sfertil}))+0.87, 0), z=if(${sdepth} > 0, (0.19*log(${sdepth}))+1, 0), a=if(x <= 0 || z <= 0, 0, ((((x*y*z)/3)*${maxwheat})/${fieldsperhectare})), if(a < 0, 0, a))", quiet = "True", tempwheatreturn = tempwheatreturn, precip = precip, sfertil = oldfert, sdepth = oldsdepth, maxwheat = maxwheat, fieldsperhectare = fieldsperhectare)
#Calculate barley yield map (kg/cell)
tempbarleyreturn = "%stemporary_barley_yields_map" % pid
grass.mapcalc("${tempbarleyreturn}=eval(x=if(${precip} > 0, (0.48*log(${precip}))+1.51, 0), y=if(${sfertil} > 0, (0.34*log(${sfertil}))+1.09, 0), z=if(${sdepth} > 0, (0.18*log(${sdepth}))+0.98, 0), a=if(x <= 0 || z <= 0, 0, ((((x*y*z)/3)*${maxbarley})/${fieldsperhectare})), if(a < 0, 0, a))", quiet = "True", tempbarleyreturn = tempbarleyreturn, precip = precip, sfertil = oldfert, sdepth = oldsdepth, maxbarley = maxbarley, fieldsperhectare = fieldsperhectare)
#Create the desired cereal mix
tempcerealreturn = "%stemporary_cereal_yields_map" %pid
grass.mapcalc("${tempcerealreturn}=if(isnull(${agcatch}), null(), (((1-${agmix})*${tempwheatreturn})+(${agmix}*${tempbarleyreturn})) )", quiet = "True", tempcerealreturn = tempcerealreturn, agmix = agmix, tempwheatreturn = tempwheatreturn, tempbarleyreturn = tempbarleyreturn, agcatch = agcatch)
grass.message("Figuring out the farming plan for this year...")
#gather some stats from yields maps in order to make an estimate of number of farm plots...
cerealstats2 = grass.parse_command('r.univar', flags = 'ge', map = tempcerealreturn)
#"Fuzz up" the agent's memory of past yields and shortfalls. We do this by padding the actual values of these things to a randomly generated percentage that is drawn from a gaussian probability distribution with mu of the mean value and sigma of 0.0333. This means that the absolute max/min pad can only be up to +- %10 of the mean value (eg. at the 3-sigma level of a gaussian distribution with sigma of 0.0333), and that pad values closer to 0% will be more likely than pad values close to +- 10%. This more closely models how good people are at "educated guesses" of central tendencies (i.e., it's how people "guesstimate" the "average" value). This also ensures some variation from year to year, regardless of the "optimum" solution.
if len(farmyieldmemory) is 0: #if it's the first year, then just use the fuzzed average potential yield from all cells in agcatch, and make the padded amount 1
fuzzyyieldmemory = random.gauss(float(cerealstats2["mean"]), (float(cerealstats2['mean']) * 0.0333))
fuzzydeficitmemory = -1
else:
if agentmem > (year -1):
slicer = year - 1
else:
slicer = agentmem
grass.debug("slicer: %s" % slicer)
fuzzyyieldmemory = random.gauss(numpy.mean(farmyieldmemory[slicer:]), (numpy.mean(farmyieldmemory[slicer:]) * 0.0333))
fuzzydeficitmemory = random.gauss(numpy.mean(farmingmemory[slicer:]), (numpy.mean(farmingmemory[slicer:]) * 0.0333))
#Figure out how many fields the agent thinks it needs based on current average yield
numfields = int(round(float(cerealreq) / fuzzyyieldmemory))
grass.debug("total fields should be %s" % numfields)
#discover how many cells to add or subtract based on agent memory of surplus or deficit
fieldpad = int(round(fuzzydeficitmemory / fuzzyyieldmemory))
grass.debug("fieldpad (sign reversed) %s" % fieldpad)
#add or remove the number of extra cells calculated from agent memory surplus or deficit
numfields = numfields - fieldpad
grass.debug("did numfields change? %s" % numfields)
#Then, make sure we don't exceed our catchment size or labor capacity (and reduce number of fields if necessary)
if numfields > int(round((float(cerealstats2['cells']) - float(cerealstats2["null_cells"])))):
numfields = int(round((float(cerealstats2['cells']) - float(cerealstats2["null_cells"]))))
if numfields > maxfields:
numfields = maxfields
grass.debug("did numfields hit the max and be curtailed? %s" % numfields)
#check for tenure, and do the appropriate type of tenure if asked
if use_flags['m'] is True:
grass.message("Land Tenure is ON, with MAXIMIZING strategy")
#check for first year, and zero out tenure if so
if (year + 1) == 1:
grass.run_command('r.random', quiet = 'True', input = agcatch, npoints = numfields, raster = tempfields)
grass.message('First year, so all farm fields randomly assigned')
tenuredcells = 0
droppedcells = 0
newcells = 0
else:
#make some map names
oldfields = "%s_Year_%s_Farming_Impacts_Map" % (prfx, then)
tenuredfields = "%s_Year_%s_Tenured_Fields_Map" % (prfx, now)
tempagcatch = "%stemporary_agricultural_catchment_map" % pid
grass.message('Performing yearly land tenure evaluation')
#grab stats from old fields
oldfarmstats = grass.parse_command('r.univar', flags = 'ge', map = oldfields)
#check the comparison metric and drop underperforming fields compared to average yield...
if tenuredrop >= 0:
grass.mapcalc("${tenuredfields}=if(isnull(${oldfields}), null(), if(${tempcerealreturn} >= (${meanyield}*${tenuredrop}), 1, null()))", quiet = "True", tenuredfields=tenuredfields, oldfields = oldfields, tempcerealreturn = tempcerealreturn, meanyield = oldfarmstats['mean'], tenuredrop = tenuredrop)
else: #or compared to the maximum yeild.
tenuredrop = tenuredrop * -1
grass.mapcalc("${tenuredfields}=if(isnull(${oldfields}), null(), if(${tempcerealreturn} >= (${maxyield}*${tenuredrop}), 1, null()))", quiet = "True", tenuredfields=tenuredfields, oldfields = oldfields, tempcerealreturn = tempcerealreturn, maxyield = oldfarmstats['max'], tenuredrop = tenuredrop)
#temporarily update the agricultural catchment by withholding tenured cells. All other cells in the catchment will be fair game to be chosen for the new fields.
grass.mapcalc("${tempagcatch}=if(isnull(${agcatch}), null(), if(isnull(${tenuredfields}), ${agcatch}, null() ))", quiet = "True", tempagcatch = tempagcatch, agcatch = agcatch, tenuredfields = tenuredfields)
#pull stats out to see what we did (how many old fields kept, how many dropped, and how many new ones we need this year
tenuredstats = grass.parse_command('r.univar', flags = 'ge', map = tenuredfields)
tenuredcells = int(float(tenuredstats['cells']) - float(tenuredstats['null_cells']))
droppedcells = int(float(oldfarmstats['cells']) - float(oldfarmstats['null_cells'])) - tenuredcells
if numfields-tenuredcells <= 0:
newcells = 0
tempfields = tenuredfields
else:
newcells = numfields-tenuredcells
#Now run r.random to get the required number of additional fields
tempfields1 = "%stemporary_extra_fields_map" % pid
grass.run_command('r.random', quiet = 'True', input = tempagcatch, npoints = newcells, raster = tempfields1)
#patch the new fields into the old fields
grass.run_command('r.patch', quiet = "True", input = "%s,%s" % (tempfields1,tenuredfields), output = tempfields)
grass.message("Keeping %s fields in tenure list, dropping %s underperforming fields, adding %s new fields" % (tenuredcells, droppedcells,newcells))
elif use_flags['i'] is True:
grass.message("Land Tenure is ON, with SATSFICING strategy")
#check for first year, and zero out tenure if so
if (year + 1) == 1:
grass.run_command('r.random', quiet = 'True', input = agcatch, npoints = numfields, raster = tempfields)
grass.message('First year, so all farm fields randomly assigned')
tenuredcells = 0
droppedcells = 0
newcells = 0
else:
#make some map names
oldfields = "%s_Year_%s_Farming_Impacts_Map" % (prfx, then)
tenuredfields = "%s_Year_%s_Tenured_Fields_Map" % (prfx, now)
tempcomparefields = "%stemporary_Old_Fields_yield_map" % pid
tempagcatch = "%stemporary_agricultural_catchment_map" % pid
grass.message('Performing yearly land tenure evaluation')
#make a map of the current potential yields from the old fields
grass.mapcalc("${tempcomparefields}=if(isnull(${oldfields}), null(), ${tempcerealreturn} )", quiet = "True", tempcomparefields = tempcomparefields, oldfields = oldfields, tempcerealreturn = tempcerealreturn)
#grab stats from the old fields
oldfarmstats = grass.parse_command('r.univar', flags = 'ge', map = tempcomparefields)
#make map of tenured fields for this year
grass.mapcalc("${tenuredfields}=if(isnull(${oldfields}), null(), 1)", quiet = "True", tenuredfields=tenuredfields, oldfields = oldfields)
#temporarily update the agricultural catchment by withholding tenured cells. All other cells in the catchment will be fair game to be chosen for the new fields.
grass.mapcalc("${tempagcatch}=if(isnull(${agcatch}), null(), if(isnull(${tenuredfields}), ${agcatch}, null() ))", quiet = "True", tempagcatch = tempagcatch, agcatch = agcatch, tenuredfields = tenuredfields)
#pull stats out to see what we did (how many old fields kept, how many dropped, and how many new ones we need this year
tenuredstats = grass.parse_command('r.univar', flags = 'ge', map = tenuredfields)
tenuredcells = int(float(tenuredstats['cells']) - float(tenuredstats['null_cells']))
droppedcells = 0
if numfields-tenuredcells <= 0:
newcells = 0
tempfields = tenuredfields
else:
newcells = numfields-tenuredcells
#Now run r.random to get the required number of additional fields
tempfields1 = "%stemporary_extra_fields_map" % pid
grass.run_command('r.random', quiet = 'True', input = tempagcatch, npoints = newcells, raster = tempfields1)
#patch the new fields into the old fields
grass.run_command('r.patch', quiet = "True", input = "%s,%s" % (tempfields1,tenuredfields), output = tempfields)
grass.message("Adding %s new fields" % (newcells))
else:
grass.message("Land Teure is OFF")
tenuredcells = 0
droppedcells = 0
newcells = 0
#Now run r.random to get the required number of fields
grass.run_command('r.random', quiet = 'True', input = agcatch, npoints = numfields, raster = tempfields)
#use r.surf.gaussian to cacluate fertily impacts in the farmed areas
grass.run_command('r.surf.gauss', quiet = "True", output = tempimpacta, mean = farmimpact[0], sigma = farmimpact[1])
grass.mapcalc("${fields}=if(isnull(${tempfields}), null(), ${tempimpacta})", quiet = "True", fields = fields, tempfields = tempfields, tempimpacta = tempimpacta)
#grab some yieled stats while region is still aligned to desired field size
#first make a temporary "zone" map for the farmed areas in which to run r.univar
tempfarmzone = "%stemporary_farming_zones_map" % pid
grass.mapcalc("${tempfarmzone}=if(isnull(${fields}), null(), ${tempcerealreturn})", quiet = "True", tempfarmzone = tempfarmzone, fields = fields, tempcerealreturn = tempcerealreturn)
#now use this zone map to grab some stats about the actual yields for this year's farmed fields
cerealstats = grass.parse_command('r.univar', flags = 'ge', percentile = '90', map = tempfarmzone)
#calculate some useful stats
cerealdif = float(cerealstats['sum']) - float(cerealreq)
numfarmcells = int(float(cerealstats['cells']) - float(cerealstats['null_cells']))
areafarmed = numfarmcells * float(nsfieldsize) * float(ewfieldsize)
#update agent mempory with the farming surplus or deficit from this year
farmingmemory.append(cerealdif)
farmyieldmemory.append(float(cerealstats['mean']))
#find out the percentage of agcatch that was farmed this year.
basepercent = 100 * (numfarmcells / (float(cerealstats2['cells']) - float(cerealstats2["null_cells"]) ) )
if basepercent > 100:
agpercent = 100.00
else:
agpercent = basepercent
#reset region to original resolution
grass.del_temp_region()
grass.message('We farmed %s fields, using %.2f percent of agcatch...' % (numfarmcells,agpercent))
#GENERATE GRAZING IMPACTS
grass.message("Calculating potential grazing yields")
#generate basic impact values
tempimpactg = "%stemporary_grazing_impact" % pid
grass.run_command("r.random.surface", quiet = "True", output = tempimpactg, distance = grazespatial, exponent = grazepatchy, high = maxgrazeimpact)
#Calculate temporary grazing yield map in kg/ha
tempgrazereturnha = "%stemporary_hectares_grazing_returns_map" % pid
tempgrazereturn = "%stemporary_grazing_returns_map" % pid
grass.run_command('r.recode', quiet = 'True', flags = 'da', input = oldlcov, output = tempgrazereturnha, rules = fodder_rules)
#convert to kg / cell, and adjust to impacts
grass.mapcalc("${tempgrazereturn}=(${tempgrazereturnha}/${cellperhectare}) * ${tempimpactg}", quiet = "True", tempgrazereturn = tempgrazereturn, tempgrazereturnha = tempgrazereturnha, cellperhectare = cellperhectare, tempimpactg = tempimpactg)
grass.message('Figuring out grazing plans for this year....')
#Do we graze on the stubble of agricultural fields? If so, how much fodder to we think we will get?
if use_flags['g'] is False:
#temporarily set region to match the field size again
grass.use_temp_region()
grass.run_command('g.region', quiet = 'True',nsres = nsfieldsize, ewres = ewfieldsize)
#set up a map with the right values of stubble fodder in it, and get them to the right units (fodder per farm field)
stubfod1 = "%stemporary_stubblefodder_1" % pid
stubfod2 = "%stemporary_stubblefodder_2" % pid
stubfod3 = "%stemporary_stubblefodder_3" % pid
#make map of the basic landcover value for fields (grass stubbles)
grass.mapcalc("${stubfod1}=if(${tempfarmzone}, ${farmval}, null())", quiet = "True", stubfod1 = stubfod1, farmval = farmval, tempfarmzone = tempfarmzone)
#turn that map into baseline grazing yields/ha for that landcover value
grass.run_command('r.recode', quiet = 'True', flags = 'da', input = stubfod1, output = stubfod2, rules = fodder_rules)
#Match the variability in stubble yields to that in cereal returns, and convert to yields per field
grass.mapcalc("${stubfod3}=(${stubfod2} / ${fieldsperhectare}) * (${tempcerealreturn}/${maxcereals})", quiet = "True", stubfod3 = stubfod3, stubfod2 = stubfod2, fieldsperhectare = fieldsperhectare, tempcerealreturn = tempcerealreturn, maxcereals = cerealstats['max'])
stubblestats = grass.parse_command('r.univar', flags = 'ge', percentile = '90', map = stubfod3)
#Since we grazed stubble, let's' estimate how much stubbles we got by padding the mean value to a randomly generated percentage that is drawn from a gaussian probability distribution with mu of the mean value and sigma of 0.0333. This means that the absolute max/min pad can only be up to +- %10 of the mean value (eg. at the 3-sigma level of a gaussian distribution with sigma of 0.0333), and that pad values closer to 0% will be more likely than pad values close to +- 10%. This more closely models how good people are at "educated guesses" of central tendencies (i.e., it's how people "guesstimate" the "average" value). This also ensures some variation from year to year, regardless of the "optimum" solution. Once we've done that, do we think there's still some remaining fodder needs to be grazed from wild patches? If so, how much?
if (float(fodderreq) - ( float(stubblestats['mean']) * numfarmcells )) < 0:
remainingfodder = 0
else:
remainingfodder = float(fodderreq) - ( random.gauss(float(stubblestats['mean']), (float(stubblestats['mean']) * 0.0333)) * numfarmcells )
#reset region
grass.del_temp_region()
else:
stubblestats = {"mean": '0', "sum": "0", "cells": "0", "stddev": "0", "min": "0", "first_quartile": "0", "median": "0", "third_quartile": "0", "max": "0"}
remainingfodder = float(fodderreq)
#Do we graze on the "fallowed" portions of the agricultural catchment?
tempgrazecatch = "%stemporary_grazing_catchment_map" % pid
if use_flags['f'] is True:
grass.mapcalc("${tempgrazecatch}=if(isnull(${grazecatch}), null(), if(isnull(${agcatch}), ${tempgrazereturn}, null()))", quiet = "True", tempgrazecatch = tempgrazecatch, grazecatch = grazecatch, agcatch = agcatch, tempgrazereturn = tempgrazereturn)
else:
grass.mapcalc("${tempgrazecatch}=if(isnull(${grazecatch}), null(), if(isnull(${fields}), ${tempgrazereturn}, null()))", quiet = "True", tempgrazecatch = tempgrazecatch, grazecatch = grazecatch, fields = fields, tempgrazereturn = tempgrazereturn)
#Now that we know where we are allowed to graze, how much of the grazing catchment does the agent think it needs to meet its remaining fodder requirements? First grab some general stats from the grazing catchment.
fodderstats = grass.parse_command('r.univar', flags = 'ge', percentile = '90', map = tempgrazecatch)
#"Fuzz up" the agent's memory of past yields and shortfalls. We do this by padding the actual values of these things to a randomly generated percentage that is drawn from a gaussian probability distribution with mu of the mean value and sigma of 0.0333. This means that the absolute max/min pad can only be up to +- %10 of the mean value (eg. at the 3-sigma level of a gaussian distribution with sigma of 0.0333), and that pad values closer to 0% will be more likely than pad values close to +- 10%. This more closely models how good people are at "educated guesses" of central tendencies (i.e., it's how people "guesstimate" the "average" value). This also ensures some variation from year to year, regardless of the "optimum" solution.
if len(grazeyieldmemory) is 0: #if it's the first year, then just use the fuzzed average potential yield from all cells in agcatch, and make the padded amount 1
fuzzygyieldmemory = random.gauss(float(fodderstats['mean']), (float(fodderstats['mean']) * 0.0333))
fuzzygdeficitmemory = -1
else:
fuzzygyieldmemory = random.gauss(numpy.mean(grazeyieldmemory[slicer:]), (numpy.mean(grazeyieldmemory[slicer:]) * 0.0333))
fuzzygdeficitmemory = random.gauss(numpy.mean(grazingmemory[slicer:]), (numpy.mean(grazingmemory[slicer:]) * 0.0333))
#Figure out how many grazing patches the agent thinks it needs based on current average patch yield
numfoddercells = int(round(float(remainingfodder) / fuzzygyieldmemory))
grass.debug("total graze patches should be %s" % numfoddercells)
#discover how many cells to add or subtract based on agent memory of surplus or deficit
fodderpad = int(round(fuzzygdeficitmemory / fuzzygyieldmemory))
grass.debug("fodderpad (sign reversed) %s" % fodderpad)
#add or remove the number of extra cells calculated from agent memory surplus or deficit
numfoddercells = numfoddercells - fodderpad
grass.debug("did numfoddercells change? %s" % numfoddercells)
#check if we need more cells than are in the maximum grazing catchment, and clip to that catchment if so
totgrazecells = int(float(fodderstats['cells'])) - int(float(fodderstats['null_cells']))
grass.debug("did we have to clip numfoddercells to the catchment size? %s" % numfoddercells)
#make the actual grazing impacts map
if numfoddercells > totgrazecells:
celltarget = totgrazecells
grass.mapcalc("${grazeimpacts}=if(isnull(${tempgrazecatch}), null(), if(${oldlcov} <= ${mingraze}, null(), ${tempimpactg}))", quiet = "True", grazeimpacts = grazeimpacts, tempimpactg = tempimpactg, tempgrazecatch = tempgrazecatch, oldlcov = oldlcov, mingraze = mingraze)
elif numfoddercells == 0:
grass.mapcalc("${grazeimpacts}=null()", quiet = "True", grazeimpacts = grazeimpacts)
else:
celltarget = numfoddercells
#now clip the cost surface to the grazable area (including cells with vegetation too low to graze on), and iterate through it to figure out the actual grazing area to be used this year
tempgrazecost = "%stemporary_grazing_cost_map" % pid
grass.mapcalc("${tempgrazecost}=if(isnull(${tempgrazecatch}), null(),if(${oldlcov} <= ${mingraze}, null(), ${costsurf}))", quiet = "True", tempgrazecatch = tempgrazecatch, tempgrazecost = tempgrazecost, costsurf = costsurf, oldlcov = oldlcov, mingraze = mingraze)
catchstat = [float(x) for x in grass.read_command("r.stats", quiet = "True", flags = "1n", input = tempgrazecost, separator="=", nsteps = "10000").splitlines()]
target = 0
cutoff = []
for x in sorted(catchstat):
target = target + 1
cutoff.append(x)
if target >= celltarget:
break
try:
#make the actual grazing impacts map
grass.mapcalc("${grazeimpacts}=if(${tempgrazecost} > ${cutoff}, null(), ${tempimpactg})", quiet = "True", grazeimpacts = grazeimpacts, tempimpactg = tempimpactg, tempgrazecost = tempgrazecost, cutoff = cutoff[-1])
except:
grass.fatal("Uh oh! Somethng wierd happened when figuring out this year\'s grazing catchment! Check your numbers and try again! Sorry!")
sys.exit(1)
#now get some grazing yields stats
#Check if stubble grazing got us everything we wanted, and act appropriately
if numfoddercells == 0:
grazestats = {"mean": '0', "sum": "0", "cells": "0", "stddev": "0", "min": "0", "first_quartile": "0", "median": "0", "third_quartile": "0", "max": "0"}
totalfodder = float(stubblestats['sum'])
grazepercent = 0
fodderdif = totalfodder - float(fodderreq)
areagrazed = 0
else:
#first make a temporary "zone" map for the grazed areas in which to run r.univar
tempgrazezone = "%stemporary_grazing_zones_map" % pid
grass.mapcalc("${tempgrazezone}=if(isnull(${grazeimpacts}), null(), ${tempgrazereturn})", quiet = "True", tempgrazezone = tempgrazezone, grazeimpacts = grazeimpacts, tempgrazereturn = tempgrazereturn)
#now grab the univar stats with this zone file
grazestats = grass.parse_command('r.univar', flags = 'ge', percentile = '90', map = tempgrazezone)
numgrazecells = (float(grazestats['cells']) - float(grazestats["null_cells"]) )
totalfodder = float(grazestats['sum']) + float(stubblestats['sum'])
grazepercent = 100 * (numgrazecells / totgrazecells)
fodderdif = totalfodder - float(fodderreq)
areagrazed = numgrazecells * (float(region['nsres']) * float(region['ewres']))
#update agent mempory with the grazing surplus or deficit from this year
grazingmemory.append(fodderdif)
grazeyieldmemory.append(float(grazestats['mean']))
grass.message('We got %.2f kg of fodder from stubbles, and so we grazed %.2f percent of grazecatch this year.' % (float(stubblestats['sum']),grazepercent))
#figure out how many animals and people were fed
animfed = (totalfodder / indfodreq)
peoplefed = ((float(cerealstats['sum']) / indcerreq) + (animfed / animals))
grass.message('We fed %i herd animals, and supported %i people this year...' %(animfed, peoplefed))
#write the yield stats to the stats file
grass.message('Writing some farming and grazing stats from this year....')