-
Notifications
You must be signed in to change notification settings - Fork 7
/
write_scaled_raw_scenarios.r
1669 lines (1455 loc) · 87.3 KB
/
write_scaled_raw_scenarios.r
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
# write_scaled_raw_scenarios.r
####
# California Natural and Working Lands Carbon and Greenhouse Gas
# Model (CALAND) Copyright (c) 2020, The Regents of the University of
# California, through Lawrence Berkeley National Laboratory (subject to
# receipt of any required approvals from the U.S. Dept. of Energy). All
# rights reserved.
# If you have questions about your rights to use or distribute this software,
# please contact Berkeley Lab's Intellectual Property Office at
# IPO@lbl.gov.
#
# NOTICE. This Software was developed under funding from the U.S. Department
# of Energy and the U.S. Government consequently retains certain rights. As
# such, the U.S. Government has been granted for itself and others acting on
# its behalf a paid-up, nonexclusive, irrevocable, worldwide license in the
# Software to reproduce, distribute copies to the public, prepare derivative
# works, and perform publicly and display publicly, and to permit others to do so.
####
# This software and its associated input data are licensed under a modified BSD open source license
# Please see license.txt for details
# Converts an area specific (e.g., county) raw scenario file into a scaled version for generating CALAND input files
# Also saves the conversion info for use by scale_caland_outputs.r, which converts CALAND outputs to the area-specific values
# This is designed to generate files for estimating county level effects of management
# So only the management areas and the county name are needed
# The input units will be preserved in the scaled output
# This is not needed for Seagrass as it resides in the ocean
# The input and output raw scenario files are assumed to be in the 'raw_data' folder
############# Arguments to write_scaled_raw_scenarios.r ##############
# 1. 'scen_file': name of raw scenario file for a specific area, such as a county; default is amador_example.xls
# This file must be in the raw_data folder
# If you want this and the scaled file to be in a folder within raw_data, prefix the file name with the existing folder
# For counties, the Region(s) can be specified explicitly,
# OR the entries can have "All" for the Region - in this case the approprieate regions will be determined
# 2. 'county': name of the county; default is Amador; this must match the county name in the county area file
# 3. 'units': the area units for the input scenario file; default is 'ac'; could also be 'ha'
# 4. 'county_category_areas_file': the file containing the breakdown of county areas for the land categories;
# default is area_lab_sp9_own9_2010lt15_cnty_sqm_stats.csv and this should not be changed
############### Outputs ##################
# Scaled raw scenario file for input to write_caland_inputs()
# The county name and the units will be appended to the input file name
# A file containing expected CALAND output scalars
# this is an excel file with one sheet for each input scenario sheet
# these sheets will be matched to the output files for scaling
# each sheet contains scalars only for the county land categories
# there is a scalar, an expected county area, and an expected region area for each year for each land cat
# "_scalars" will be appended to the name of the scaled scenario file name to create this file name
################# Limitations ###################
# No baseline LULCC can be happening - need to set this using write_caland_inputs()
# this is because the land category distributions are different between region and county,
# which means that the land changes cannot be scaled properly
# the scaling must be able to estimate the changes in area, otherwise the densities will not be valid
# No wildfire can be happening if restoration is happening - need to set this using write_caland_inputs()
# this is because it changes the carbon densities based on the regional fire area distribution, which cannot be mapped to the county
# annual activities and restoration activities cannot be prescribed in the same scenario; they must be in separate scenarios
# this is because carbon densities cannot be changing independently of area changes
# Land protection cannot be estimated because no baseline LULCC can be used
# No non-regeneration because this causes LULCC that cannot be mapped to the county
# this is set when CALAND is run
# If interactions between restoration practices cause the prescriptions to not be met, the carbon densities may be off
# This just throws a warnin if this may be the case, but lets it run
# This is because I can't guarantee that the estimated areas here are exactly the same as in CALAND (even though they should be)
###################### Notes ######################
# Each scenario must at least have the three urban management practices defined in order to run, with all regions and ownerships defined (can use "All"):
# Developed_all Dead_removal: must be defined as 1
# Developed_all Growth: must be defined as 1
# Developed_all Urban_forest: baseline is defined as the statewide value of 0.15, but there are also reigonal values
# The working directory needs to be the main CALAND folder where this file resides
# The output scenarios file name will be the same as the intput file name but with county name and the units ('ac' or 'ha') appended to it
# The output scalars file name will additionally have "scalars" appended to it
####################### Start script #####################
# this enables java to use up to 4GB of memory for reading and writing excel files
options(java.parameters = "-Xmx4g" )
# Load all the required packages
libs <- c( "XLConnect" )
for( i in libs ) {
if( !require( i, character.only=T ) ) {
cat( "Couldn't load", i, "\n" )
stop( "Use install.packages() to download this library\nOr use the GUI Package Installer\nInclude dependencies,
and install it for local user if you do not have root access\n" )
}
library( i, character.only=T )
}
write_scaled_raw_scenarios <- function(scen_file = "amador_example_ac.xls", county = "Amador", units = "ac",
county_category_areas_file = "area_lab_sp9_own9_2010lt15_cnty_sqm_stats.csv") {
### these have been pulled from the arguments because the writing of output scalars is incorrect and has not been fixed for non-county applications
ISCOUNTY = TRUE
avail_area = NA
# 4. 'ISCOUNTY': flag to denote whether this scenario is for a county (TRUE), or for a specified area (FALSE)
# Counties are directly supported for checking areas and scaling automatically
# Other specified areas (ISCOUNTY=FALSE) require additional information about total available area for the project,
# and are restricted to a single specific land category, e.g. Klamath Private Forest
# this is for exploratory research only and has not been tested or validated
# there may be some inconsistencies between project and region for restoration source areas in this mode
# default is TRUE, and any other value will return an error
# 5. 'avail_area': the total available area of the project within a single land category; default is NA
# this applies only to project-level, non-county simulations
# for an annual area practice, this is the total area of the available project area within a land category
# for a restoration practice, this is the initial area of the restored land type within the available project area
# this must be > 0
cat("Start write_scaled_raw_scenarios at", date(), "\n")
in_dir = "./raw_data/"
out_dir = in_dir
xltag = ".xls"
csvtag = ".csv"
# remove all whitespace from the county input string for use in output files names
# this works on a vector as well
ctag = gsub("//s", "", county)
out_file = paste0(out_dir, substr(scen_file, 1, regexpr(".xls", scen_file)-1), "_", ctag, "_", units, xltag)
out_scalar_file = paste0(out_dir, substr(scen_file, 1, regexpr(".xls", scen_file)-1), "_", ctag, "_", units, "_scalars", xltag)
# the column headers are on line 12, for both input and output files
start_row = 12
# the header goes from row 1 to row 10, and is only the first column
last_head_row = 10
# regions
reg_names = c("Central_Coast", "Central_Valley", "Delta", "Deserts", "Eastside", "Klamath", "North_Coast", "Sierra_Cascades", "South_Coast")
num_reg = length(reg_names)
# project-level estimation is only for exploratory research
if (!ISCOUNTY) {
stop("Project-level estimation has not been tested or validated and is for exploratory research only!\n")
}
# check that county is a region name for project-level and that avail_area > 0
if (!ISCOUNTY) {
reg_def = which(reg_names == county)
if (length(reg_def) == 0) {
stop("Arguemnt 'county' must be a region for non-county cases\n")
}
if (is.na(avail_area) | avail_area <= 0) {
stop("Arguemnt 'avail_area' must be > 0\n")
}
}
# land types
lt_names = c("Water", "Ice", "Barren", "Sparse", "Desert", "Shrubland", "Grassland", "Savanna", "Woodland", "Forest", "Meadow",
"Coastal_marsh", "Fresh_marsh", "Cultivated", "Developed_all", "Seagrass")
num_lt = length(lt_names)
# ownerships
own_names = c("BLM", "DoD", "Easement", "Local_gov", "NPS", "Other_fed", "Private", "State_gov", "USFS_nonwild")
num_own = length(own_names)
# urban
urban_man = c("Dead_removal", "Urban_forest", "Growth")
num_urban_man = length(urban_man)
# restoration
restoration_man = c("Restoration", "Afforestation", "Reforestation")
num_restoration_man = length(restoration_man)
# restoration land types
# the first forest type is for afforestation
# the second forest land type is for reforestation
restoration_lt = c("Meadow", "Fresh_marsh", "Coastal_marsh", "Woodland", "Forest", "Forest")
num_restoration_lt = length(restoration_lt)
# restoration sources
restoration_sources = array(dim=c(num_restoration_lt, num_restoration_lt))
num_restoration_sources = c(4, 1, 1, 2, 2, 1)
restoration_sources[,] = NA
restoration_sources[1,1:num_restoration_sources[1]] = c("Shrubland", "Grassland", "Savanna", "Woodland")
restoration_sources[2,1:num_restoration_sources[2]] = c("Cultivated")
restoration_sources[3,1:num_restoration_sources[3]] = c("Cultivated")
restoration_sources[4,1:num_restoration_sources[4]] = c("Grassland", "Cultivated")
restoration_sources[5,1:num_restoration_sources[5]] = c("Shrubland", "Grassland")
restoration_sources[6,1:num_restoration_sources[6]] = c("Shrubland")
meadow_rest_index = 1
source_totals_names = c("Shrubland", "Grassland", "Savanna", "Woodland", "Cultivated")
woodland_st_index = 4
num_sources = length(source_totals_names)
st_sum = array(dim=num_sources)
st_sum_scaled = array(dim=num_sources)
# conversion factors
sqm2ha = 1/10000
sqm2ac = 1/4046.8564
# read in the area file (square meters)
# then convert it to the desired units
areas = read.csv(paste0(in_dir, county_category_areas_file), header=FALSE, stringsAsFactors=FALSE)
colnames(areas) <- c("reg_code", "Region", "lt_code", "Land_Type", "own_code", "Ownership", "cnty_code", "County", "area_sqm")
# drop the code columns
areas$reg_code = NULL
areas$lt_code = NULL
areas$own_code = NULL
areas$cnty_code = NULL
# filter out 'no data' rows
areas = areas[areas$Region != "no data" & areas$Land_Type != "no data" & areas$Ownership != "no data" & areas$County != "no data",]
# convert units
if (units == "ac") {
areas$area = areas$area_sqm * sqm2ac
} else {
areas$area = areas$area_sqm * sqm2ha
}
areas$area_sqm = NULL
scenin_wrkbk = loadWorkbook(paste0(in_dir,scen_file))
# worksheet/table names
scenin_sheets = getSheets(scenin_wrkbk)
num_scenin_sheets = length(scenin_sheets)
# NA values need to be converted to numeric
# the warnings thrown by readWorksheet below are ok because they just state that the NA string can't be converted a number so it is
# converted to NA value
scenin_df_list <- list()
for (i in 1:num_scenin_sheets) {
scenin_df_list[[i]] <- readWorksheet(scenin_wrkbk, i, startRow = start_row, colTypes = c(rep("character",4), rep("numeric",50)),
forceConversion = TRUE)
} # end for i loop to read in scenarios
###### read the scenario headers for the outputs
scen_head_df_list <- list()
for (i in 1:num_scenin_sheets) {
scen_head_df_list[[i]] <- readWorksheet(scenin_wrkbk, i, startRow = 1, endRow = last_head_row, header = FALSE)
}
# get the column headers in order
col_order = colnames(scenin_df_list[[1]])
# get the column headers in order for the scalar output
col_order_scalar = c(col_order, "County", "cnty_lc_area", "lc_area", "man_scalar")
# first copy the input worksheets to the outputs to set up the lists
out_scen_df_list = scenin_df_list
out_scen_sheets = scenin_sheets
# this one is for the table with the scalar values
out_man_df_list = scenin_df_list
# this is the output scalar file
out_scalar_df_list = scenin_df_list
# get the total land category areas
# there are no NA values
lc_areas = aggregate(area ~ Region + Land_Type + Ownership, areas, FUN = sum)
names(lc_areas)[ncol(lc_areas)] <- "lc_area"
# determine total region-land type areas
reg_lt_areas = aggregate(area ~ Region + Land_Type, areas, FUN = sum)
names(reg_lt_areas)[ncol(reg_lt_areas)] <- "reg_lt_area"
# get the county area values
# there will always be a county name even if it isn't used for project-level
# subset county land cat areas
cnty_lc_areas = areas[areas$County == county,]
names(cnty_lc_areas)[ncol(cnty_lc_areas)] <- "cnty_lc_area"
num_cnty_lc = nrow(cnty_lc_areas)
# determine county region-land type areas
cnty_reg_lt_areas = aggregate(cnty_lc_area ~ Region + Land_Type, cnty_lc_areas, FUN = sum)
names(cnty_reg_lt_areas)[ncol(cnty_reg_lt_areas)] <- "cnty_reg_lt_area"
num_cnty_reg_lt = nrow(cnty_reg_lt_areas)
# determine county ownership-land type areas
cnty_own_lt_areas = aggregate(cnty_lc_area ~ Ownership + Land_Type, cnty_lc_areas, FUN = sum)
names(cnty_own_lt_areas)[ncol(cnty_own_lt_areas)] <- "cnty_own_lt_area"
num_cnty_own_lt = nrow(cnty_own_lt_areas)
# determine county land type areas
cnty_lt_areas = aggregate(cnty_lc_area ~ Land_Type, cnty_lc_areas, FUN = sum)
names(cnty_lt_areas)[ncol(cnty_lt_areas)] <- "cnty_lt_area"
num_cnty_lt = nrow(cnty_lt_areas)
# merge these county data into a single data frame
cnty_areas = merge(cnty_lc_areas, cnty_reg_lt_areas, by = c("Region", "Land_Type"), all.x = TRUE)
cnty_areas = merge(cnty_areas, cnty_own_lt_areas, by = c("Ownership", "Land_Type"), all.x = TRUE)
cnty_areas = merge(cnty_areas, cnty_lt_areas, by = c("Land_Type"), all.x = TRUE)
# get the names of the regions in this county
if (!ISCOUNTY) {
county_regions = county
} else {
county_regions = unique(cnty_lc_areas$Region)
}
num_county_regions = length(county_regions)
# store the project-level land cats here
project_lcs = NULL
# loop over the scenario sheets
for (s in 1: num_scenin_sheets) {
# need to start with a fresh output scalar df for each scenario
# the scalar values for the output files are needed for all land cats in the county
# for project-level only the affected land cats are scaled and used - subset and adjust this at the end as necessary
output_scalars = merge(cnty_lc_areas, lc_areas, by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
output_scalars$initial_scalar = output_scalars$cnty_lc_area / output_scalars$lc_area
# Split annual, restoration, and urban records here to check that project-level operates on only annual or restoration
# not sure how to check for urban-only for project-level
annual = out_scen_df_list[[s]][out_scen_df_list[[s]]$Management!="Restoration" & out_scen_df_list[[s]]$Management!="Afforestation" &
out_scen_df_list[[s]]$Management!="Reforestation" & out_scen_df_list[[s]]$Land_Type!="Developed_all",]
cumulative = out_scen_df_list[[s]][out_scen_df_list[[s]]$Management=="Restoration" | out_scen_df_list[[s]]$Management=="Afforestation" |
out_scen_df_list[[s]]$Management=="Reforestation",]
urban = out_scen_df_list[[s]][out_scen_df_list[[s]]$Land_Type == "Developed_all",]
# first check to make sure that annual and restoration are not happening at the same time
if (nrow(annual) > 0 & nrow(cumulative) > 0) {
stop("Specify only annual practices or only restoration practices in a single scenario\n")
}
if (!ISCOUNTY) {
if (nrow(annual) > 0 & nrow(cumulative) > 0) {
stop("Specify only annual practices or only a restoration practice, for a single land category, for project-level estimation\n")
}
}
###### annual area practices
# these are entered as the annual value during each year, including the start and end years
# these require area scaling based on:
# new man area = orig man area * land category area / land category county(or avail-project) area
# need to save the inverse of these ratios in a table by land category for output scaling
# orig man area is checked against available area and is adjusted if necessary, with warnings output to ?terminal/log?
# counties are distributed proportionally among relevant regions if 'All' is entered in Region column
# for projects, only a single, specific land catergory (Region, own, land type) is valid, and the avail area is checked and adjusted
if(nrow(annual) > 0) {
if (!ISCOUNTY) {
# full land category must be specified
if ( !(all(annual$Region %in% reg_names) & all(annual$Ownership %in% own_names) & all(annual$Land_Type %in% lt_names)) ) {
stop("For project-level application the land category has to be completely specified for Region, Ownership, and Land Type\n")
}
# make sure only one land category is annually managed
check_df = unique(annual[,1:3])
if (nrow(check_df) > 1) {
stop("Only one land category can be managed for project-level application\n")
}
# make sure the total project available area value does not exceed the total physical area
check_df = merge(annual, lc_areas, by = c("Region", "Land_Type", "Ownership", ), all.x = TRUE)
if (avail_area > check_df$lc_area[1]) {
stop("The avail_area for project-level application must be <= the total initial physical land category area of ", check_df$lc_area[1], units, "\n")
}
# make sure that the input region matches the record region
if (annual$Region[1] != county) {
stop("The management record region must match the county input argument\n")
}
# put available project area here in the same column name as for counties below
# one input area will define the available area because project-level must be specified by a single land category
# cnty_lc_area is the total available project area for this land type in this ownership and region
# this will be available source area for a restoration land type
annual_reg_own = annual
annual_reg_own$County = county
annual_reg_own$cnty_lc_area = avail_area
} else { # county
# separate any "All" region entries so that the appropriate regions can be specified
# ownerships can be mixed
# the records will be expanded to full land categories
#annual_all_reg = annual[annual$Region == "All",]
#annual_reg = annual[annual$Region != "All",]
#if(nrow(annual_all_reg) > 0 & nrow(annual_reg) > 0){
# stop("Regions must EITHER be specified OR all Region entries must be set to 'All'\n")
#}
# now separate into the four reg-own cases
annual_all = annual[annual$Region == "All" & annual$Ownership == "All",]
annual_all_reg = annual[annual$Region == "All" & annual$Ownership != "All",]
annual_all_own = annual[annual$Region != "All" & annual$Ownership == "All",]
annual_reg_own_orig = annual[annual$Region != "All" & annual$Ownership != "All",]
# distribute records across regions if necessary
# expand the records and adjust management areas
# then strip the added columns to bind them together
# then merge the lc area for scalar calc
annual_reg_own = NULL
# start with the allreg-allown case and append these to the region-allown case below
if(nrow(annual_all) > 0){
# distribute the records across regions
num_recs = nrow(annual_all)
new_records = NULL
for (r in 1:num_recs) {
temp_df = annual_all[r,]
temp_df = merge(temp_df,
unique(cnty_areas[,c("Region", "Land_Type", "cnty_reg_lt_area", "cnty_lt_area")]),
by = c("Land_Type"), all.x = TRUE)
temp_df$Region.x = temp_df$Region.y
temp_df$Region.y = NULL
colnames(temp_df)[colnames(temp_df) == "Region.x"] = "Region"
temp_df$start_area = temp_df$start_area * temp_df$cnty_reg_lt_area / temp_df$cnty_lt_area
temp_df$end_area = temp_df$end_area * temp_df$cnty_reg_lt_area / temp_df$cnty_lt_area
temp_df$cnty_reg_lt_area = NULL
temp_df$cnty_lt_area = NULL
new_records = rbind(new_records, temp_df)
}
# make sure the new table has the same column order
new_records = new_records[,col_order]
if (nrow(annual_all_own) > 0) {
annual_all_own = rbind(annual_all_own, new_records)
} else {
annual_all_own = new_records
}
} # end if all-all case
# do the reg-allown case
if(nrow(annual_all_own) > 0){
# distribute the records across ownerships
num_recs = nrow(annual_all_own)
new_records = NULL
for (r in 1:num_recs) {
# first make sure that the record regions match the county
cr_ind = which(county_regions == annual_all_own$Region[r])
if (length(cr_ind) == 0) {
stop("Record region ", annual_all_own$Region[r], " is not in county ", county, " region list: ", county_regions, "\n")
}
temp_df = annual_all_own[r,]
temp_df = merge(temp_df,
unique(cnty_areas[,c("Region", "Land_Type", "Ownership", "cnty_reg_lt_area", "cnty_lc_area")]),
by = c("Region", "Land_Type"), all.x = TRUE)
temp_df$Ownership.x = temp_df$Ownership.y
temp_df$Ownership.y = NULL
colnames(temp_df)[colnames(temp_df) == "Ownership.x"] = "Ownership"
temp_df$start_area = temp_df$start_area * temp_df$cnty_lc_area / temp_df$cnty_reg_lt_area
temp_df$end_area = temp_df$end_area * temp_df$cnty_lc_area / temp_df$cnty_reg_lt_area
temp_df$cnty_reg_lt_area = NULL
temp_df$cnty_lc_area = NULL
new_records = rbind(new_records, temp_df)
}
# make sure the new table has the same column order
new_records = new_records[,col_order]
annual_reg_own = rbind(annual_reg_own, new_records)
} # end if region-allown case
# do the allreg-own case
if(nrow(annual_all_reg) > 0){
# distribute the records across regions
num_recs = nrow(annual_all_reg)
new_records = NULL
for (r in 1:num_recs) {
temp_df = annual_all_reg[r,]
temp_df = merge(temp_df,
unique(cnty_areas[,c("Region", "Land_Type", "Ownership", "cnty_own_lt_area", "cnty_lc_area")]),
by = c("Ownership", "Land_Type"), all.x = TRUE)
temp_df$Region.x = temp_df$Region.y
temp_df$Region.y = NULL
colnames(temp_df)[colnames(temp_df) == "Region.x"] = "Region"
temp_df$start_area = temp_df$start_area * temp_df$cnty_lc_area / temp_df$cnty_own_lt_area
temp_df$end_area = temp_df$end_area * temp_df$cnty_lc_area / temp_df$cnty_own_lt_area
temp_df$cnty_own_lt_area = NULL
temp_df$cnty_lc_area = NULL
new_records = rbind(new_records, temp_df)
}
# make sure the new table has the same column order
new_records = new_records[,col_order]
annual_reg_own = rbind(annual_reg_own, new_records)
} # end if allregion-own case
# bind the newly expanded records with any fully specified records (reg-own case)
if (nrow(annual_reg_own_orig) > 0) {
for (r in 1:nrow(annual_reg_own_orig)) {
# first make sure that the record regions match the county
cr_ind = which(county_regions == annual_reg_own_orig$Region[r])
if (length(cr_ind) == 0) {
stop("Record region ", annual_reg_own_orig$Region[r], " is not in county ", county, " region list: ", county_regions, "\n")
}
}
annual_reg_own = rbind(annual_reg_own, annual_reg_own_orig)
}
# now all records should be fully expanded to land category
# add the county and county lc area cols
annual_reg_own = merge(annual_reg_own, cnty_lc_areas, by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
# need to make sure that all management land cats are in the county
bad_inds = which(is.na(annual_reg_own$County))
if (length(bad_inds) > 0) {
for (b in 1:length(bad_inds)) {
cat("County land category does not exist for management: ", annual_reg_own$Region[bad_inds[b]], ",", annual_reg_own$Land_Type[bad_inds[b]], ",", annual_reg_own$Ownership[bad_inds[b]], "\n")
stop()
}
}
} # end else county
# now check to make sure that max annual managed areas do not exceed initial available areas
# just check the max managed area over time against the initial available area
# this won't capture limitations due to LULCC
# need to sum simultaneous practices within each land category
unique_lc = unique(annual_reg_own[,c(1:3)])
for (lc in 1:nrow(unique_lc)) {
# extract the records for this land cat
check_df = annual_reg_own[annual_reg_own$Region == unique_lc$Region[lc] & annual_reg_own$Ownership == unique_lc$Ownership[lc] &
annual_reg_own$Land_Type == unique_lc$Land_Type[lc],]
first_year = min(check_df$start_year)
last_year = max(check_df$end_year)
# loop over years
for (y in first_year:last_year) {
check_df$yearcol = 0
check_df$yearcol[check_df$start_year <= y & check_df$end_year >= y] =
apply(check_df[check_df$start_year <= y & check_df$end_year >= y,c("start_area", "end_area")], 1, max)
# sum the col and check it
total = sum(check_df$yearcol)
if (total > check_df$cnty_lc_area[1]) {
cat("Summed managed annual area ", total, " exceeds initial available area", check_df$cnty_lc_area[1], " for:\n")
stop("\t", check_df$Region[1], ", ", check_df$Ownership[1], ", ", check_df$Land_Type[1], "\n")
}
} # end for y loop over years
} # end for lc loop to check managed areas within available area
# now scale the management areas to representative regional areas
# calculate the management scaling ratio region:county for upscaling
# first get the region level land cat area
annual_reg_own = merge(annual_reg_own, lc_areas, by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
annual_reg_own$man_scalar = annual_reg_own$lc_area / annual_reg_own$cnty_lc_area
# scale the values
annual_reg_own$start_area = annual_reg_own$start_area * annual_reg_own$man_scalar
annual_reg_own$end_area = annual_reg_own$end_area * annual_reg_own$man_scalar
if (!ISCOUNTY) {
# save the info for this scenario so that scalars can be subset and adjusted
# add the source columns for binding with other scenarios that may have restoration
annual_reg_own$src_cnty_lc_area = NA
annual_reg_own$src_lc_area = NA
project_lcs = rbind(project_lcs, annual_reg_own)
}
} # end if annual practices
###### restoration practices
# these are entered as cumulative totals by the end of the target end year, so start_area should always be zero
# these require area scaling based on the expected net restored land type area to ensure proper carbon density calculations:
# base scaling factor = expected area of land category in region / expected area of land category in county
# the restored area constitutes a carbon density change, but the converted area maintains its carbon density
# (because there are not other drivers of density change based on the stated limitations)
# this expected area is used because the source distributions are different between county and region
# the inverse of this base scaling factor should be sufficient for output scaling, but is needed for each year
# this base scaling factor needs to be adjusted for management area scaling to make sure that:
# the ratio of region change:total is the same as the ratio of county change:total
# it is important to get this right for both restoration types and sources as it determines the density values
# main challenge is to figure out the woodland and meadow scaling because woodland is a source of meadow
# caounty target cum area is checked against available county source area and an error is thrown if not enough available area
# this is based on the county land cat distribution
# also, the scaled source area needs to be checked to ensure scalability; an error is thrown if there is not enough
# this is based on the region land cat distribution
# this means that in some cases the target restored area for the county may not be able to be simulated
# this can occur if a considerable proprtion of the source area is in the county compared to the region as a whole
# counties are distributed proportionally among relevant regions if 'All' is entered in Region column
# for projects, the sources are for only a single, specific land catergory (Region, own, land type), and the avail area is checked and adjusted
if(nrow(cumulative) > 0) {
# Both Reforestation and Afforestion cannot be done at the same time
# technically, this is just within a particular land category, but throw error regardless of this
# this isn't caught in write_caland_inputs(), but it can crash CALAND
aff_df = cumulative[cumulative$Management == "Afforestation",]
ref_df = cumulative[cumulative$Management == "Reforestation",]
if (nrow(aff_df) > 0 & nrow(ref_df) > 0) {
cat("Afforestation and Reforestation cannot be prescribed in the same scenario!\n")
stop("It is recommended to use Reforestation for all forest area expansion activities (which converts Shrubland), unless you are sure you want to afforest Grassland and Shrubland.\n")
}
# make sure start area is zero
if (sum(cumulative$start_area) > 0) {
stop("Start area for all restoration practices must be zero because this is a cumulative definition\n")
}
# make sure that any delta fresh marsh records are fully specified
# and add them to output_scalars if so
fm_df = unique(cumulative[cumulative$Land_Type == "Fresh_marsh", c("Region", "Land_Type", "Ownership")])
if ( nrow(fm_df) > 0 ) {
for (r in 1:nrow(fm_df)) {
if ( !(fm_df$Region[r] != "All" & fm_df$Ownership[r] != "All") ) {
stop("Delta fresh marsh restoration must be fully specified by region and ownership\n")
} else {
# if delta fresh marsh is included it needs to be added to output_scalars
new_row = output_scalars[1,]
new_row$Region = fm_df$Region[r]
new_row$Land_Type = "Fresh_marsh"
new_row$Ownership = fm_df$Ownership[r]
new_row$County = county
new_row$cnty_lc_area = 0
new_row$lc_area = 0
new_row$initial_scalar = 1
output_scalars = rbind(output_scalars, new_row)
}
}
}
if (!ISCOUNTY) {
# full land category must be specified
if ( !(all(cumulative$Region %in% reg_names) & all(cumulative$Ownership %in% own_names) & all(cumulative$Land_Type %in% lt_names)) ) {
stop("For project-level application the land category has to be completely specified for Region, Ownership, and Land Type\n")
}
# make sure only one land category is cumulatively managed
check_df = unique(cumulative[,1:3])
if (nrow(check_df) > 1) {
stop("Only one land category can be managed for project-level application\n")
}
# make sure that the input region matches the record region
if (cumulative$Region[1] != county) {
stop("The management record region must match the county input argument\n")
}
# check that this land cat exists
# if delta fresh marsh it is ok if associated cultivated exists
cum_lc_area = lc_areas$lc_area[lc_areas$Region == cumulative$Region[1] & lc_areas$Ownership == cumulative$Ownership[1] & lc_areas$Land_Type == cumulative$Land_Type[1]]
cult_lc_area = lc_areas$lc_area[lc_areas$Region == "Delta" & lc_areas$Ownership == cumulative$Ownership[1] & lc_areas$Land_Type == "Cultivated"]
if ( !(length(cum_lc_area) > 0) ) {
if ( !(cumulative$Region[1] == "Delta" & cumulative$Land_Type[1] == "Fresh _marsh" & length(cult_lc_area) > 0) ) {
cat("Restoration type does not exist, so it cannot be restored:\n")
cat("\t", cumulative$Region[1], ", ", cumulative$Ownership[1], ", ", cumulative$Land_Type[1], ",",
cumulative$Management[1], "\n")
stop("\tCannot estimate this land category restoration\n")
}
}
# put available project area here in the same column name as for counties below
# one input area will define the available area because project-level must be specified by a single land category
# cnty_lc_area is the total initial area for this restored land type in this ownership and region
# lc_area will be added below with the region-level total initial area of the restored land type
# source columns are added and filled below
cum_reg_own = cumulative
cum_reg_own$County = county
cum_reg_own$cnty_lc_area = avail_area
} else { # county
# this expansion could probably be done together with the annual
# separate any "All" region entries so that the appropriate regions can be specified
# ownerships can be mixed
# the records will be expanded to full land categories
#cum_all_reg = cumulative[cumulative $Region == "All",]
#cum_reg = cumulative[cumulative $Region != "All",]
#if(nrow(cum_all_reg) > 0 & nrow(cum_reg) > 0){
# stop("Regions must EITHER be specified OR all Region entries must be set to 'All'\n")
#}
# now separate into the four reg-own cases
cum_all = cumulative[cumulative$Region == "All" & cumulative$Ownership == "All",]
cum_all_reg = cumulative[cumulative$Region == "All" & cumulative$Ownership != "All",]
cum_all_own = cumulative[cumulative$Region != "All" & cumulative$Ownership == "All",]
cum_reg_own_orig = cumulative[cumulative$Region != "All" & cumulative$Ownership != "All",]
# distribute records across regions if necessary
# expand the records and adjust management areas
# then strip the added columns to bind them together
# then merge the lc area for scalar calc
cum_reg_own = NULL
# start with the allreg-allown case and append these to the region-allown case below
if(nrow(cum_all) > 0){
# distribute the records across regions
num_recs = nrow(cum_all)
new_records = NULL
for (r in 1:num_recs) {
temp_df = cum_all[r,]
temp_df = merge(temp_df,
unique(cnty_areas[,c("Region", "Land_Type", "cnty_reg_lt_area", "cnty_lt_area")]),
by = c("Land_Type"), all.x = TRUE)
temp_df$Region.x = temp_df$Region.y
temp_df$Region.y = NULL
colnames(temp_df)[colnames(temp_df) == "Region.x"] = "Region"
temp_df$start_area = temp_df$start_area * temp_df$cnty_reg_lt_area / temp_df$cnty_lt_area
temp_df$end_area = temp_df$end_area * temp_df$cnty_reg_lt_area / temp_df$cnty_lt_area
temp_df$cnty_reg_lt_area = NULL
temp_df$cnty_lt_area = NULL
new_records = rbind(new_records, temp_df)
}
# make sure the new table has the same column order
new_records = new_records[,col_order]
if (nrow(cum_all_own) > 0) {
cum_all_own = rbind(cum_all_own, new_records)
} else {
cum_all_own = new_records
}
} # end if all-all case
# do the reg-allown case
if(nrow(cum_all_own) > 0){
# distribute the records across ownerships
num_recs = nrow(cum_all_own)
new_records = NULL
for (r in 1:num_recs) {
# first make sure that the record regions match the county
cr_ind = which(county_regions == cum_all_own$Region[r])
if (length(cr_ind) == 0) {
stop("Record region ", cum_all_own$Region[r], " is not in county ", county, " region list: ", county_regions, "\n")
}
temp_df = cum_all_own[r,]
temp_df = merge(temp_df,
unique(cnty_areas[,c("Region", "Land_Type", "Ownership", "cnty_reg_lt_area", "cnty_lc_area")]),
by = c("Region", "Land_Type"), all.x = TRUE)
temp_df$Ownership.x = temp_df$Ownership.y
temp_df$Ownership.y = NULL
colnames(temp_df)[colnames(temp_df) == "Ownership.x"] = "Ownership"
temp_df$start_area = temp_df$start_area * temp_df$cnty_lc_area / temp_df$cnty_reg_lt_area
temp_df$end_area = temp_df$end_area * temp_df$cnty_lc_area / temp_df$cnty_reg_lt_area
temp_df$cnty_reg_lt_area = NULL
temp_df$cnty_lc_area = NULL
new_records = rbind(new_records, temp_df)
}
# make sure the new table has the same column order
new_records = new_records[,col_order]
cum_reg_own = rbind(cum_reg_own, new_records)
} # end if region-allown case
# do the allreg-own case
if(nrow(cum_all_reg) > 0){
# distribute the records across regions
num_recs = nrow(cum_all_reg)
new_records = NULL
for (r in 1:num_recs) {
temp_df = cum_all_reg[r,]
temp_df = merge(temp_df,
unique(cnty_areas[,c("Region", "Land_Type", "Ownership", "cnty_own_lt_area", "cnty_lc_area")]),
by = c("Ownership", "Land_Type"), all.x = TRUE)
temp_df$Region.x = temp_df$Region.y
temp_df$Region.y = NULL
colnames(temp_df)[colnames(temp_df) == "Region.x"] = "Region"
temp_df$start_area = temp_df$start_area * temp_df$cnty_lc_area / temp_df$cnty_own_lt_area
temp_df$end_area = temp_df$end_area * temp_df$cnty_lc_area / temp_df$cnty_own_lt_area
temp_df$cnty_own_lt_area = NULL
temp_df$cnty_lc_area = NULL
new_records = rbind(new_records, temp_df)
}
# make sure the new table has the same column order
new_records = new_records[,col_order]
cum_reg_own = rbind(cum_reg_own, new_records)
} # end if allregion-own case
# bind the newly expanded records with any fully specified records (reg-own case)
if (nrow(cum_reg_own_orig) > 0) {
# first make sure that the record regions match the county
for (r in 1:nrow(cum_reg_own_orig)) {
cr_ind = which(county_regions == cum_reg_own_orig$Region[r])
if (length(cr_ind) == 0) {
stop("Record region ", cum_reg_own_orig$Region[r], " is not in county ", county, " region list: ", county_regions, "\n")
}
}
cum_reg_own = rbind(cum_reg_own, cum_reg_own_orig)
}
# now all records should be fully expanded to land category
# add the county and county lc area cols
cum_reg_own = merge(cum_reg_own, cnty_lc_areas, by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
# need to make sure that all management land cats are in the county
# except for delta fresh marsh is ok here; source availability is checked below
bad_inds = which(is.na(cum_reg_own$County))
if (length(bad_inds) > 0) {
for (b in 1:length(bad_inds)) {
if ( cum_reg_own$Region[bad_inds[b]] == "Delta" & cum_reg_own$Land_Type[bad_inds[b]] == "Fresh_marsh" ) {
cum_reg_own$County[bad_inds[b]] = county
cum_reg_own$cnty_lc_area[bad_inds[b]] = 0
} else {
cat("County land category does not exist for management: ", cum_reg_own$Region[bad_inds[b]], ",", cum_reg_own$Land_Type[bad_inds[b]], ",", cum_reg_own$Ownership[bad_inds[b]], "\n")
stop()
}
}
}
} # end else county
# first get the region level land cat area
cum_reg_own = merge(cum_reg_own, lc_areas, by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
# this is the managment scalar for upscaling region:county
# calculate the management scaling ratio
# if restoration type is delta fresh marsh need to set the lc area to zero and the management scalar to 1
cum_reg_own$man_scalar = cum_reg_own$lc_area / cum_reg_own$cnty_lc_area
cum_reg_own$lc_area[cum_reg_own$Region == "Delta" & cum_reg_own$Land_Type == "Fresh_marsh"] = 0
cum_reg_own$man_scalar[cum_reg_own$Region == "Delta" & cum_reg_own$Land_Type == "Fresh_marsh"] = 1
# add the source cnty_lc_area and lc_area cols as NA because they will be filled below with total source area
# the county level source area is based on county land cat distribution
# the region level source area is based on region land cat distribution
cum_reg_own$src_cnty_lc_area = NA
cum_reg_own$src_lc_area = NA
# now check to make sure that total prescribed managed areas do not exceed initial available areas
# just check the end managed area against the initial available area of the source categories
# this won't capture limitations due to LULCC - there shouldn't be baseline LULCC!!!!
# but multiple restoration activities may compete for source area
# need to sum like source category needs first, by region-ownership, using the county distributions
# sum the scaled like source category needs, using the regional distributions
# also reset the src_cnty_lc_area to the resepctive source area totals
# also sum up the region-level land cat source areas to check for availability
unique_lc = unique(cum_reg_own[,c(1:3)])
num_lc = nrow(unique_lc)
unique_ro = unique(unique_lc[,c(1,3)])
num_ro = nrow(unique_ro)
source_totals = array(dim=c(num_ro, num_sources))
scaled_source_totals = array(dim=c(num_ro, num_sources))
# store the total source areas needed for restoration, by reg-own and source land type, for county and scaled
# the county areas are base on the county distribution
# the scaled areas are based on the reg-own distribution
source_totals[,] = 0
scaled_source_totals[,] = 0
# store the total prescribed restoration area for each restoration land cat
rest_totals = array(dim=num_lc)
rest_totals[] = 0
# store the reg-own land category areas of sources by restoration land cat
reg_own_source_lc_areas = array(dim=num_lc)
reg_own_source_lc_areas[] = 0
for (lc in 1:num_lc) {
# extract the records for this land cat
check_df = cum_reg_own[cum_reg_own$Region == unique_lc$Region[lc] & cum_reg_own$Ownership == unique_lc$Ownership[lc] &
cum_reg_own$Land_Type == unique_lc$Land_Type[lc],]
# determine the region-ownership index for source totals
ro_reg_ind = which(unique_ro$Region == unique_lc$Region[lc])
ro_own_ind = which(unique_ro$Ownership == unique_lc$Ownership[lc])
ro_ind = intersect(ro_reg_ind, ro_own_ind)
# determine restoration type index and the management type
# Forest has two types of restoration, so two rt_ind values will be returned
# parse this and deal with the different practice types
# afforestation and reforestation are not allowed in the same scenario - this is checked above
rt_ind = which(restoration_lt == unique_lc$Land_Type[lc])
if (length(rt_ind) == 2) {
if (check_df$Management[1] == "Afforestation") {
rt_ind = 5
practice = "Afforestation"
} else {
rt_ind = 6
practice = "Reforestation"
}
} else if (length(rt_ind) == 1) {
practice = "Restoration"
} else if (length(rt_ind) == 0) {
stop("Incorrect restoration land type ", unique_lc$Land_Type[lc], "\n")
} # end if-else restoration type
# sum the restoration area
rest_totals[lc] = sum(check_df$end_area)
# sum the source areas
st_den = 0 # this is to store the county-level total source area for comparison below
st_sum[] = 0
# sum the scaled source areas
st_den2 = 0 # this is to store the region-level total source area for comparison below
st_sum_scaled[] = 0
for (src in 1:num_restoration_sources[rt_ind]) {
st_ind = which(source_totals_names == restoration_sources[rt_ind,src])
st_num = cnty_lc_areas$cnty_lc_area[cnty_lc_areas$Region == unique_lc$Region[lc] & cnty_lc_areas$Ownership == unique_lc$Ownership[lc] &
cnty_lc_areas$Land_Type == source_totals_names[st_ind]]
if(length(st_num) == 0) { st_num = 0 }
st_den = st_den + st_num
st_sum[st_ind] = st_sum[st_ind] + rest_totals[lc] * st_num
# sum up the region-level total source area
temp = lc_areas$lc_area[lc_areas$Region == unique_lc$Region[lc] & lc_areas$Ownership == unique_lc$Ownership[lc] &
lc_areas$Land_Type == source_totals_names[st_ind]]
if(length(temp) == 0) { temp = 0 }
# calculate the scaled source needs
st_den2 = st_den2 + temp
st_sum_scaled[st_ind] = st_sum_scaled[st_ind] + rest_totals[lc] * check_df$man_scalar[1] * temp
} # end for summing up individual sources for non-forest restoration
# keep track of sources in this region-ownership
source_totals[ro_ind,] = source_totals[ro_ind,] + st_sum / st_den
scaled_source_totals[ro_ind,] = scaled_source_totals[ro_ind,] + st_sum_scaled / st_den2
# store the reg-own land category area of source for this restoration land cat
reg_own_source_lc_areas[lc] = st_den2
# set to zero if no source area
source_totals[ro_ind,] <- replace(source_totals[ro_ind,], is.nan(source_totals[ro_ind,]), 0.0)
scaled_source_totals[ro_ind,] <- replace(scaled_source_totals[ro_ind,], is.nan(scaled_source_totals[ro_ind,]), 0.0)
cum_reg_own$src_cnty_lc_area[cum_reg_own$Region == unique_lc$Region[lc] & cum_reg_own$Ownership == unique_lc$Ownership[lc] &
cum_reg_own$Land_Type == unique_lc$Land_Type[lc]] = st_den
cum_reg_own$src_lc_area[cum_reg_own$Region == unique_lc$Region[lc] & cum_reg_own$Ownership == unique_lc$Ownership[lc] &
cum_reg_own$Land_Type == unique_lc$Land_Type[lc]] = st_den2
# loop over the source types and check against the actual county level source availability
# for project-level the check is the total source against the total avail area input
if (!ISCOUNTY) {
# check available initial area against physical initial area
if (cum_reg_own$cnty_lc_area[1] > cum_reg_own$lc_area[1]) {
cat("Input initial restoration type area ", cum_reg_own$cnty_lc_area[1], " exceeds physical regional initial restoration type area ",
cum_reg_own$lc_area[1], " for restoration type:\n")
cat("\t", cum_reg_own$Region[1], ", ", cum_reg_own$Ownership[1], ", cum_reg_own$lc_area[1] ", cum_reg_own$Land_Type[1], "\n")
stop()
}
# check that there is enough local source area
if (sum(source_totals[1,]) > cum_reg_own$src_cnty_lc_area[1]) {
cat("Summed restoration source area ", sum(source_totals[1,]), " exceeds initial available source area ", cum_reg_own$src_cnty_lc_area[1],
" for restoration type:\n")
cat("\t", cum_reg_own$Region[1], ", ", cum_reg_own$Ownership[1], ", ", cum_reg_own$Land_Type[1], "\n")
stop()
}
# check that there is enough total source area for scaling
if (sum(scaled_source_totals[1,]) > cum_reg_own$src_lc_area[1]) {
cat("Summed scaled restoration source area", sum(scaled_source_totals[1,]), " exceeds physical regional initial source area",
cum_reg_own$lc_area[1], " for restoration type:\n")
cat("\t", cum_reg_own$Region[1], ", ", cum_reg_own$Ownership[1], ", cum_reg_own$lc_area[1]", cum_reg_own$Land_Type[1], "\n")
stop()
}
} else {
# county level
# check that initial area for any land cat restoration type is > 0
# if delta fresh marsh it is ok for it to be zero (as set above), as the scaling is just equal to 1 because it doesn't exist anywhere
z_ind = which(cum_reg_own$cnty_lc_area == 0)
if ( length(z_ind) > 0 ) {
for ( z in 1:z_ind) {
z_df = cum_reg_own[z_ind[z],]
fm_df = z_df[z_df$Region == "Delta" & z_df$Land_Type == "Fresh_marsh",]
if ( !(nrow(fm_df) > 0) ) {
cat("Initial restoration type county area is zero for land category and management:\n")
cat("\t", cum_reg_own$Region[z_ind[z]], ", ", cum_reg_own$Ownership[z_ind[z]], ", ", cum_reg_own$Land_Type[z_ind[z]], ",",
cum_reg_own$Management[z_ind[z]], "\n")
stop("\tCannot estimate this land category restoration for this county\n")
}
}
}
# loop over source types to check for sufficent local source area
for (src in 1:length(source_totals[ro_ind,])) {
source_avail = cnty_lc_areas$cnty_lc_area[cnty_lc_areas$Region == unique_lc$Region[lc] & cnty_lc_areas$Ownership == unique_lc$Ownership[lc] &
cnty_lc_areas$Land_Type == source_totals_names[src]]
if(length(source_avail) == 0) { source_avail = 0 }
if (source_totals[ro_ind, src] > source_avail) {
cat("Summed restoration source area ", source_totals[ro_ind, src], " exceeds initial available area ", source_avail, " for source type:\n")
cat("\t", unique_lc$Region[lc], ", ", unique_lc$Ownership[lc], ", ", source_totals_names[src], "\n")
stop("\tHalted at restoration type ", unique_lc$Land_Type[lc], "\n")
}
} # end for loop to check source totals against available area
# loop over source types to check for sufficent scaled areas
for (src in 1:length(scaled_source_totals[ro_ind,])) {
source_avail = lc_areas$lc_area[lc_areas$Region == unique_lc$Region[lc] & lc_areas$Ownership == unique_lc$Ownership[lc] &
lc_areas$Land_Type == source_totals_names[src]]
if(length(source_avail) == 0) { source_avail = 0 }
if (scaled_source_totals[ro_ind, src] > source_avail) {
cat("Summed scaled restoration source area ", scaled_source_totals[ro_ind, src], " exceeds initial physical area ", source_avail, " for source type:\n")
cat("\t", unique_lc$Region[lc], ", ", unique_lc$Ownership[lc], ", ", source_totals_names[src], "\n")
stop("\tHalted at restoration type ", unique_lc$Land_Type[lc], "\n")
}
} # end for loop to check source totals against available area
} # end if project else county for checking available source area
} # end for lc loop to check for managed areas within available area
# calculate modified management scalars for woodland and meadow restoration if necessary
# this happens only if both woodland and meadow restoration exist in the same reg-own
# this is to maintain carbon density fidelity between county and region by:
# ensuring that the ratio of scaled restoration area to initial reg-own area is the same as the corresponding county-level ratio
wl_inds = which(unique_lc$Land_Type == "Woodland")
md_inds = which(unique_lc$Land_Type == "Meadow")
if (length(wl_inds) > 0 & length(md_inds) > 0) {
# check for matching region
reg_matches = intersect(unique_lc$Region[wl_inds], unique_lc$Region[md_inds])
if (length(reg_matches) > 0) {
# check for matching ownership by region
for ( r in 1:length(reg_matches) ) {
reg_own_matches = intersect(unique_lc$Ownership[unique_lc[wl_inds, "Region"] == reg_matches[r]],
unique_lc$Ownership[unique_lc[md_inds, "Region"] == reg_matches[r]])
if (length(reg_own_matches) > 0) {
for ( o in 1:length(reg_own_matches) ) {
# subset this reg-own for woodland and meadow
# there may be multiple rows for each if different time periods are prescribed
ro_df = cum_reg_own[cum_reg_own$Region == reg_matches[r] & cum_reg_own$Ownership == reg_own_matches[o] &
(cum_reg_own$Land_Type == "Woodland" | cum_reg_own$Land_Type == "Meadow"),]
# calculate and a management adjustment scalar and apply it to man_scalar
# this is done per region-ownership
# this should be true at this point, but check anyway
if (nrow(ro_df) > 0) {
# determine the region-ownership index for source totals
ro_reg_ind = which(unique_ro$Region == reg_matches[r])