forked from nbest937/pad-us_nlcd
-
Notifications
You must be signed in to change notification settings - Fork 0
/
pad-us_nlcd.txt
1791 lines (1327 loc) · 47.2 KB
/
pad-us_nlcd.txt
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
Synthesis of a disaggregated carbon inventory for the conterminous United States
================================================================================
Author: Neil Best, Joshua Elliott <nbest@ci.uchicago.edu, jelliott@ci.uchicago.edu>
Date:
Table of Contents
=================
1 Introduction
1.1 [exact name of biofuels program]
1.2 Motivation
1.3 Analysis Framework
1.4 Summary of Methodology
2 Underlying datasets
2.1 Describe NBCD ALDB and BAWH
2.2 cite NBCD web site and in-review paper
2.3 NLCD land use/land cover and canopy density
2.4 Dependency of NBCD on NLCD
3 Analysis Procedure
3.1 Reproject the $5'$ grid cells
3.2 Download and prepare inputs
3.3 Load the data into GRASS
3.4 Replace null values in NBCD with zeroes
3.5 Aggregate the statistics
4 Plots
5 Appendix
5.1 Software Environment
1 Introduction
===============
1.1 TODO [exact name of biofuels program]
------------------------------------------
some prose about the biofuels program?
1.2 Motivation
---------------
To apply the information contained in the National Biomass and Carbon
Dataset (NBCD) to an analysis of the land use change implications of
expanding biofuels production it is necessary to add qualitative depth
to the data. There are components of the total above-ground live dry
biomass (ALDB) that should be treated differently when modelling
economic decisions about land use change and land management. Rather
than simply adding up the ALDB for each geographical unit in a given
set of spatial boundaries that total biomass can be broken out into
categories by overlaying information from other sources. To add
information about the physical characteristics of the biomass
indicated at a given location and its cultural or ecological context
we use the National Land Cover Database (NLCD). Similarly to
distignuish areas based on their legal protection status we bring in
the Protected Areas Database of the United States (PAD-US). By
assigning attributes from NLCD and PAD-US to the corresponding biomass
inventory indicated by a given pixel in the NBCD data we are able to
disaggregate the total biomass into categories given by the Cartesian
product of the individual datasets' categories. More detail on this in
section ______.
* TODO cite NBCD, NLCD, PAD-US
* TODO cross-reference to section on data set details
1.3 Analysis Framework
-----------------------
A secondary goal of this work is to demonstrate a framework for
reproducible research. Unless you are browsing the source code
repository of this project you are reading a document that is the
result of an automated process that downloads the base data, performs
some necessary transformations of those datasets, performs the
computations from top to bottom, updates the figures, tables, and
maps, and exports the complete document to a PDF file. The particular
framework chosen for this work uses the Emacs editor, Org-mode for
outlining and literate programming, the Geospatial Data Abstraction
Library (GDAL) for manipulation of vector and raster datasets, GRASS
for analysis and tabulation of GIS data, and R for aggregation and
formatting of the data.
* TODO cite Emacs, Org-mode, R, GRASS, GDAL, raster, data.table
Full software version information is available in the appendix.
## library( raster)
library( raster, lib.loc="~/src/R/lib/")
setOptions( progress= "text")
library( plyr)
library( stringr)
library( ascii)
options( asciiType= "org")
library( reshape)
library( Hmisc)
library( data.table)
## overwriteRasters <- TRUE
1.4 Summary of Methodology
---------------------------
The initial goal of this work was to assign protection status to areas
of the various land cover classes. Because PAD-US has coverage
everywhere that NLCD does we extended the analysis to include Alaska,
Hawaii, and Puerto Rico. However, at present NBCD only provides data
for the conterminous United States (cUSA) so the ultimate result of
this process only covers that region. The tabulation of protected
areas by cover clasee for the outlying areas is provided for
completeness leaving open the possibility that more data might become
available in the future.
To satisfy various modelling and reporting requirements we have
adopted three modes of disaggregation of the woody biomass: state,
county, and 5 arc-minute ($5'$, one twelfth of a degree) grid cells.
Clearly the first two are dependent and nested -- the total biomass in
a state is the sum of the biomass of its counties. Apart from the
political boundaries it is common practice to subdivide the globe or
regions thereof into units based on geographic coordinates,
i.e. longitude and latitude. At $5'$ resolution a grid cell is
roughly 10 km on a side at the equator. The variation in area of
these grid cells must be accounted for in any comparison among or
summation across grid cells, but this complication is accepted as a
trade-off against working in familiar map coordinates.
In all cases we chose to reproject the units of disaggregation, being
the states, counties, and grid cells, into the native representation
of the input raster data (NBCD and NLCD), which happen to be in
regional Albers equal-area projections chosen to minimize area
distortions for their respective extents: cUSA, Alaska, Hawaii, and
Purto Rico. The political boundaries and protection status
delineations are provided in vector format making them less
computationally intensive to reproject than their raster counterparts.
After reprojection we convert the political and protection status
boundaries into their raster representation at the same resolution as
the NBCD and NLCD data, thereby assigning additional qualitative
categories to each NBCD/NLCD pixel. From a programming perspective it
was expedient to reproject the grid cells in raster space but they
just as well could have been defined as vector entities as well.
Therefore assigning NBCD pixels to their respective $5'$ grid cells
was perhaps more expensinve in terms of processing time to perform but
simpler from a conceptual and implementation perspective. By bringing
our reporting units into the coordinate space of the NBCD data we
guarantee that the sum total of all given biomass will be conserved.
* DONE tangle out the R code and run it from the Makefile
* TODO Add an emacs-lisp function call to the makefile to evaluate the R code
This has to happen before tangling out the noweb components. However,
it is possible that the need for this has gone away with new versions
of Org since this project began.
2 Underlying datasets
======================
This section describes the input datasets that go into this analysis
in greater detail than what was offerd in the introduction.
Conceptually the analysis employs land use / land cover information
from NLCD, legal protection status from PAD-US, and political
boundaries from the US Census to disaggregate undifferentiated
estimates of woody biomass, i.e. ALDB, which is assumed to be 50%
carbon by mass, taken from NBCD. The analysis also employs another
spatial divisions of disaggregation, the $5'$ grid cells, which are
also defined here for completeness.
2.1 TODO Describe NBCD ALDB and BAWH
-------------------------------------
2.2 TODO cite NBCD web site and in-review paper
------------------------------------------------
Kellndorfer, J., Walker, W., LaPoint, E., Bishop, J., Cormier, T.,
Fiske, G., Hoppus, M., Kirsch, K., and Westfall, J. 2012. NACP
Aboveground Biomass and Carbon Baseline Data (NBCD 2000),
U.S.A., 2000. Data set. Available on-line [http://daac.ornl.gov] from
ORNL DAAC, Oak Ridge, Tennessee,
U.S.A. [http://dx.doi.org/10.3334/ORNLDAAC/1081] or at
[http://www.whrc.org/mapping/nbcd/]
Kellndorfer, J.M., Walker, W.S., Kirsch, K.M., Fiske, G., Bishop, J.,
LaPoint, E., Hoppus, M, & Westfall, J. The National Biomass and Carbon
data set 2000 (NBCD 2000): A high resolution baseline for the
conterminous U.S. Remote Sensing of Environment, Submitted.
2.3 NLCD land use/land cover and canopy density
------------------------------------------------
* TODO Describe NLCD land use/land cover and canopy density
2.4 Dependency of NBCD on NLCD
-------------------------------
One important factor that influenced our operational choices in
carrying out the analysis is the dependence between NLCD and NBCD.
NBCD's ALDB densities are only calculated where NLCD's canopy density
is non-zero and WHRC's radar-derived basal area weighted height
(BAWH) is also non-zero. Because NBCD is modelled based on the first
version of the 2001 NLCD and that is the only version for which
canopy density is provided it is meaningless to disaggregate ALDB
based on subsequent versions of NLCD. The classification accuracy of
NLCD 2001 version 1 is described in _____.
* TODO cite carbon mass fraction of ALDB from WHRC
* TODO cite NLCD 2001 v1 accuracy assessment
* TODO explain which NLCD classes are modelled for ALDB
3 Analysis Procedure
=====================
Conceptually our method reduces to a brute-force tabulation exercise.
We are simply constructing classifications of the following forms:
- grid ID, land use/cover, protection status
- state, county, land use/cover, protection status
By tabulating NBCD pixels on these bases we are able sum up the ALDB
densities for a given set of identifiers in the above schema. For
this method to work it is necessary to convert all of the variables
into rasters of the same projection, extent, and resolution so that
there is a one-to-one relationship between pixels among the variables.
For this analysis we have chosen to work in the Albers equal-area
projection used by NBCD and NLCD to avoid performing any
transformation on them. PAD-US and the political boundaries are
provided in vector form so they are reprojected and converted to
rasters. Similarly the $5'$ grid cells that are regular in
geographic coordinates are reprojected into irregular areas in the
target projection so that each NBCD/NLCD pixel is assigned a grid
cell ID.
3.1 Reproject the $5'$ grid cells
----------------------------------
For convenience we adopt a scheme for identifying the $5'$ grid cells
using sequential integers rather than decimal coordinates. By
stipulating that any grid we wish to consider in our models is either
a global grid or a subset thereof the only parameter that must be
tracked is the resolution. Cell number one has its upper left corner
at 180 degrees west longitude and 90 degrees north latitude (-180,90)
by convention. From there the cell numbers increase to the east and
then resume at the western edge when eastern edge is reached. =R='s
=raster= package provides a convenient abstraction for specifying and
instantiating such simply formulated rasters. The following code
performs these steps:
1. Create a global raster in geogrpahic coordinates by default.
2. Set the resolution to $5'$.
3. Assign values to its cells from a sequence that runs from one to
the number of cells.
4. Write the data to a file.
world <- raster()
res( world) <- 5/60
world[ ] <-
1:ncell( world)
world <-
writeRaster(
world, "data/grid5minWorld.tif",
datatype= "INT4U",
overwrite= overwriteRasters)
The cUSA study area covers only a subset of these cells. Using the
=gdalwarp= utility from GDAL we can specify a target projection,
extent, and resolution that match NBCD/NLCD.
cusaProj='+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 cusaProj'
cusaProj+='+lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs'
gdalwarp -overwrite -of HFA \
-t_srs ${cusaProj} \
-te -2493045 177285 2342655 3310005 \
-tr 30 30 -co "COMPRESSED=YES" \
data/grid5minWorld.tif data/grid5minAeaCUSA.img
3.2 Download and prepare inputs
--------------------------------
3.2.1 NLCD
~~~~~~~~~~~
The first version of the 2001 NLCD is available as a single raster
file so no special preparation is required.
3.2.2 NBCD
~~~~~~~~~~~
The second version of the NBCD is divided according to "mapping zones"
and so must be patched together to form a continuous coverage. We do
this using GDAL's =gdalbuildvrt= to build a virtual raster from the
collection of downloaded mapping zone tiles. The caveat which comes
with this approach is that the individual NBCD rasters include a small
buffer outside the boundary of the respective mapping zone which means
that there will be one or more NBCD rasters that contain estimates for
pixels on either side of the mapping zone boundaries. Because the
model that produces the NBCD estimates has a stochastic component
estimates for a given pixel within these overlapping buffer zones are
not constrained to have the same value. From [the `gdalbuildvrt' page on gdal.org]:
ORG-BLOCKQUOTE-START
If there is some amount of spatial overlapping between files, the
order may depend on the order they are inserted in the VRT file, but
this behaviour should not be relied on.
ORG-BLOCKQUOTE-END
The best action we could take to maintain reproducibility in this
situation is to ensure that the arguments to the =gdalbuildvrt=
utility, the names of the NBCD tiles, are given to in lexicographic
order.
[the `gdalbuildvrt' page on gdal.org]: http://gdal.org/gdalbuildvrt.html
* TODO either copy or include code from NBCD submodule
3.2.3 PAD-US
~~~~~~~~~~~~~
3.2.4 State and county boundaries
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3.3 Load the data into GRASS
-----------------------------
Due to the volume of data involved in computations at continental
extents and 30m resolution processing performance becomes a critical
factor. At the time of this writing the performance of =GRASS='s
=r.stats= utility far outstrips that of Hijman & van Etten's =raster=
package for =R= in computing pixel count tabulations. Therefore it is
worth the trouble to populate as =GRASS= database and convolve our
input data sets there.
* TODO check if =raster= has improved tabulation performance.
3.3.1 start GRASS
~~~~~~~~~~~~~~~~~~
cd ~/Desktop/pad-us_nlcd
export GISRC=./.grassrc6
grass64 -text data/grass/cUSA/PERMANENT
* TODO check that GRASS processing works out of the box
Make sure that the bare GRASS database in the Git repo is enough to
start the project cleanly.
3.3.2 Load the NBCD data
~~~~~~~~~~~~~~~~~~~~~~~~~
# g.gisenv set=LOCATION_NAME=cUSA
# g.gisenv set=MAPSET=PERMANENT
eval $(g.gisenv)
# export GRASS_MESSAGE_FORMAT=plain
r.in.gdal --overwrite -e \
input=nbcd/data/nbcdAldb.vrt output=nbcdAldb
r.in.gdal --overwrite -e \
input=nbcd/data/nbcdBawh.vrt output=nbcdBawh
3.3.3 load the $5'$ grid cell IDs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
r.in.gdal in=data/grid5minAeaCUSA.vrt out=grid_5min
3.3.4 laod the PAD-US data
~~~~~~~~~~~~~~~~~~~~~~~~~~~
r.in.gdal input=pad-us/data/PADUS1_2_cUSA_GAP.tif output=gap
3.3.5 Load the NLCD data
~~~~~~~~~~~~~~~~~~~~~~~~~
r.in.gdal input=nlcd/pr_landcover_wimperv_10-28-08_se5.img output=Nlcd01v1PR
r.in.gdal --overwrite -e \
input=nlcd/data/nlcd2001_canopy_mosaic_1-29-08/nlcd_canopy_mosaic_1-29-08.img \
output=canopy
3.3.6 Load the state and county rasters
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
r.in.gdal --overwrite -e \
input=data/cusaStatesAea.img output=states
r.in.gdal --overwrite -e \
input=data/cusaCountiesAea.img output=counties
3.4 Replace null values in NBCD with zeroes
--------------------------------------------
g.region rast=Nlcd01v1
r.mask -o input=Nlcd01v1 maskcats="1 thru 95"
r.mapcalc nbcdMask='eval( nbcd=canopy > 0 && nbcdBawh > 0, if( isnull( nbcd), 0, nbcd))'
# echo grid_5min,Nlcd01v1,gap,aldb,n > data/statsNbcdNlcd01v1Grid5min.csv && \
# r.stats -Nc input=grid_5min,Nlcd01v1,gap,nbcdAldb fs=, >> data/statsNbcdNlcd01v1Grid5min.csv 2> data/statsNbcdNlcd01v1Grid5min.err &
echo state,county,Nlcd01v1,gap,nbcd,aldb,n > data/statsNbcdNlcd01v1County.csv \
&& r.stats -Nc input=states,counties,Nlcd01v1,gap,nbcdMask,nbcdAldb fs=, \
>> data/statsNbcdNlcd01v1County.csv \
2> data/statsNbcdNlcd01v1County.err &
3.5 Aggregate the statistics
-----------------------------
3.5.1 by 5' grid cells
~~~~~~~~~~~~~~~~~~~~~~~
stats <-
read.csv( "statsNbcd.csv",
na.strings= "*",
col.names= c( "cell", "nlcd", "gap", "nbcd", "n"),
colClasses= c("numeric", "factor", "factor", "numeric"))
stats <-
within(
stats,
{ levels( gap) <- c( levels( gap), "0")
gap[ is.na( gap)] <- "0"
gap <- combine_factor( gap, c(0,1,1,1,0))
levels( gap) <- c( "no", "yes")
nbcd[ is.na( nbcd)] <- 0
})
dt <- data.table( stats)
setkey( dt, cell, nlcd, gap)
wm <- dt[, list( wm= weighted.mean( nbcd, n)), by= "cell,nlcd,gap"]
wmCt <-
cast(
data= wm,
formula= cell ~ gap + nlcd,
## fun.aggregate= sum,
## margins= "grand_col",
value= "wm" )
write.csv(
format.df(
wmCt,
cdec= c( 0, rep( 1, ncol( wmCt) - 1)),
numeric.dollar= FALSE,
na.blank= TRUE),
row.names= FALSE,
file= "nbcdFiaAldb.csv",
quote= FALSE)
zip( "pad-us_nlcd_nbcd.zip", "fracscUSA.csv")
zip( "pad-us_nlcd_nbcd.zip", "nbcdFiaAldb.csv")
+ TODO convert NAs to zeros for \*Fr and \*Ha in CSVs and SHPs
+ TODO trim spaces in char data frames before writing CSVs
3.5.2 load r.stats output
~~~~~~~~~~~~~~~~~~~~~~~~~~
library( reshape)
library( Hmisc)
library( data.table)
library( stringr)
library( ggplot2)
library( foreign)
stateAttrs <-
read.dbf( "shp/tl_2010_us_state10.dbf")
stateNames <-
data.table( stateAttrs[, c( "STATEFP10", "STUSPS10", "NAME10")])
setnames(
stateNames,
names( stateNames),
c( "state", "usps", "name"))
setkey( stateNames, state)
keycols <- c(
"state", "county", "nlcd",
"gap", "nbcd", "aldb")
rawCountyStats <-
read.csv(
"data/statsNbcdNlcd01v1County.csv",
na.strings= "*",
header= TRUE,
col.names= c( keycols, "n"),
colClasses= c(
"character", "character", "character",
"integer", "integer", "integer", "integer"))
rawCountyStats <- data.table( rawCountyStats)
setkeyv( rawCountyStats, keycols)
{
rawCountyStats[ is.na( state), state := "0"]
rawCountyStats[ is.na( county), county := "0"]
rawCountyStats[, state :=
str_pad(
as.character( state),
2, pad= "0")]
rawCountyStats[, county :=
str_pad(
as.character( county),
3, pad= "0")]
rawCountyStats[ is.na( gap), gap := 0L]
rawCountyStats[ gap == 4L, gap := 0L]
rawCountyStats[ gap != 0L, gap := 1L]
rawCountyStats[ is.na( aldb), aldb := 0L]
rawCountyStats[, gap := as.logical( gap)]
rawCountyStats[, nbcd := as.logical( nbcd)]
rawCountyStats[ !nbcd, aldb := NA]
}
setkeyv( rawCountyStats, keycols)
3.5.3 collapse the records for GAP values that are no longer differentiated
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rawCountyStats <-
rawCountyStats[, list( n= sum( n)),
keyby= keycols ]
rawStateStats <-
rawCountyStats[, list( n= sum( n)),
keyby= keycols[ -2] ]
## same thing in functional form
## rawStateStats <- `[`(
## x= rawCountyStats,
## j= list( n= sum( n)),
## keyby= keycols[ -2])
Show the low-density forests that we intend to backfill.
qplot(
data= rawCountyStats[ nlcd %in% as.character( c(41:43, 90))],
x= aldb,
weight= n,
binwidth= 5,
xlim= c( 0, 200),
geom= "histogram",
group= gap,
fill= gap,
position= "dodge",
## position= "stack",
facets= nlcd ~ .)
[file:pad-us_nlcd/images/forestHist.png]
qplot(
data= rawCountyStats[ nlcd %in% as.character( c(41:43, 90))],
x= aldb,
weight= n,
binwidth= 1,
xlim= c( 0, 20),
geom= "histogram",
group= gap,
fill= gap,
position= "dodge",
## position= "stack",
facets= nlcd ~ .)
[file:pad-us_nlcd/images/forestHistDetail.png]
3.5.4 repeat for states
~~~~~~~~~~~~~~~~~~~~~~~~
* calculate the average densities for original and modified pixel counts
+ TODO Figure out how to do self-join in functional form
statsNbcdState <- `[`(
x= rawStateStats,
j= list(
aldb= weighted.mean( aldb, n),
n= sum( n),
ha= sum(n) * 30^2 / 10^4),
keyby= "state,nlcd,gap,nbcd")
stateAreas <- `[`(
x= statsNbcdState ,
j= list( totHa= sum(ha)),
keyby= "state")
statsNbcdState <-
statsNbcdState[ stateAreas][, frac:= ha / totHa][, totHa := NULL]
## statsNbcdState <- `[`(
## x= `[`(
## x= statsNbcdState,
## j= list( stateAreas)) ,
## j= frac:= ha / totHa)
## `[`( x= statsNbcdState, j= stateAreas)
## `[`( x= statsNbcdState,
## j= list(
## names( statsNbcdState),
## frac= ha / sum( .SD[, ha])),
## keyby= "state")
To backfill those NLCD/GAP combinations we must calculate national averages.
nbcdMean <-
statsNbcdState[ n != n2][, list(
aldb= weighted.mean( aldb2, n, na.rm= TRUE)),
keyby= "nlcd,gap"]
ascii(
cast(
nbcdMean,
nlcd ~ gap,
value= "aldb"),
digits = 1,
include.rownames= FALSE)
nlcd FALSE TRUE
------+-------+-------
41 111.4 119.5
42 103.7 120.5
43 110.2 134.0
90 96.4 91.1
* calculate the grand total
statsNbcdState[ I( nbcd), list( Gt= sum( aldb * ha) / 10^9)]
* write out the state data in long form
write.csv(
statsNbcdState[, list( state,nlcd,gap,nbcd,aldb,n,ha,frac)],
row.names= FALSE,
file= "pad-us_nlcd/nbcdStateSerial.csv",
quote= FALSE)
zip( "pad-us_nlcd.zip", "pad-us_nlcd/nbcdStateSerial.csv")
* cross-tabulate the state data using backfilled densities
nbcdStateAldb <-
data.table(
cast(
data= statsNbcdState,
subset= nbcd,
formula= state ~ gap + nlcd,
value= "aldb"),
key= "state")
setnames(
nbcdStateAldb,
colnames( nbcdStateAldb),
str_replace(
str_replace(
colnames(nbcdStateAldb),
"TRUE_", "yes"),
"FALSE_", "no"))
nbcdStateAldbMeans <-
data.table(
cast(
data=
statsNbcdState[
I( nbcd),
list( aldbAve= weighted.mean( aldb, ha)),
by= c( "state", "gap")],
formula= state ~ gap,
value= "aldbAve"),
key= "state")
setnames(
nbcdStateAldbMeans,
c( "FALSE", "TRUE"),
c( "noAll", "yesAll"))
nbcdStateGapFrac <-
data.table(
cast(
data= statsNbcdState,
formula= state ~ gap,
value= "frac",
fun.aggregate= sum,
na.rm = TRUE),
key= "state")
setnames(
nbcdStateGapFrac,
c( "FALSE", "TRUE"),
c( "noAllFr", "yesAllFr"))
nbcdStateGapHa <-
data.table(
cast(
data= statsNbcdState,
formula= state ~ gap,
value= "ha",
fun.aggregate= sum,
na.rm= TRUE),
key= "state")
setnames(
nbcdStateGapHa,
c( "FALSE", "TRUE"),
c( "noAllHa", "yesAllHa"))
nbcdStateGapHa <-
data.table(
cast(
data= statsNbcdState[, list( ha=sum( ha)),
keyby= "state,gap,nbcd"],
formula= state ~ gap + nbcd,
value= "ha"),
key= "state")
setnames(
nbcdStateGapHa,
c( "FALSE_FALSE", "FALSE_TRUE", "TRUE_FALSE", "TRUE_TRUE"),
c( "noAllHaN", "noAllHa", "yesAllHaN", "yesAllHa"))
nbcdStateFrac <-
data.table(
cast(
data= statsNbcdState,
subset= nbcd,
formula= state ~ gap + nlcd,
value= "frac"),
key= "state")
setnames(
nbcdStateFrac,
colnames( nbcdStateFrac)[ -1],
paste(
str_replace(
str_replace(
colnames( nbcdStateFrac)[ -1],
"TRUE_", "yes"),
"FALSE_", "no"),
"Fr", sep= ""))
nbcdNullStateFrac <-
data.table(
cast(
data= statsNbcdState,
subset= !nbcd,
formula= state ~ gap + nlcd,
value= "frac"),
key= "state")
setnames(
nbcdNullStateFrac,
colnames( nbcdNullStateFrac)[ -1],
paste(
str_replace(
str_replace(
colnames( nbcdNullStateFrac)[ -1],
"TRUE_", "yes"),
"FALSE_", "no"),
"FrN", sep= ""))
nbcdStateHa <-
data.table(
cast(
data= statsNbcdState,
subset= nbcd,
formula= state ~ gap + nlcd,
value= "ha"),
key= "state")
setnames(
nbcdStateHa,
colnames( nbcdStateHa)[ -1],
paste(
str_replace(
str_replace(
colnames( nbcdStateHa)[ -1],
"TRUE_", "yes"),
"FALSE_", "no"),
"Ha", sep= ""))
nbcdNullStateHa <-
data.table(
cast(
data= statsNbcdState,
subset= !nbcd,
formula= state ~ gap + nlcd,
value= "ha"),
key= "state")
setnames(
nbcdNullStateHa,
colnames( nbcdNullStateHa)[ -1],
paste(
str_replace(
str_replace(
colnames( nbcdNullStateHa)[ -1],
"TRUE_", "yes"),
"FALSE_", "no"),
"HaN", sep= ""))
nbcdState <-
nbcdStateAldb[ nbcdStateAldbMeans]
nbcdState <-
nbcdState[ nbcdStateGapFrac][ nbcdStateGapHa]
nbcdState <-
nbcdState[ nbcdStateFrac][ nbcdStateHa]
nbcdState <-
nbcdState[ nbcdNullStateFrac][ nbcdNullStateHa]
setnames(
nbcdState,
"state", "fips")
setcolorder(
nbcdState,
order( colnames( nbcdState)))
nbcdStateChar <-
str_trim(
format.df(
nbcdState,
cdec= sapply(
colnames( nbcdState),
function( x) {
ifelse(
x == "fips", 0,
ifelse(
str_detect( x, "HaN?$"), 1,
ifelse(
str_detect( x, "FrN?$"), 3,
1)))
}),
numeric.dollar= FALSE,
na.blank= TRUE))
write.csv(
nbcdStateChar,
row.names= FALSE,
file= "pad-us_nlcd/nbcdState.csv",
quote= FALSE)
zip( "pad-us_nlcd.zip", "pad-us_nlcd/nbcdState.csv")
options(useFancyQuotes = FALSE)
cat(
sapply(
colnames( nbcdState),
function( x) {
dQuote(
ifelse(
x == "fips", "String(2)",
ifelse(
str_detect( x, "HaN?$"),
"Real(10.1)",
ifelse(
str_detect( x, "FrN?$"),
"Real(5.3)",
"Real(5.1)"))))
}),
sep= ",",
file= "pad-us_nlcd/nbcdState.csvt")
options(useFancyQuotes = TRUE)
ogr2ogr <-
paste(
"ogr2ogr -overwrite -progress -sql",
sprintf(
"\"select %s from cusaStatesAea a",
paste( colnames( nbcdState), collapse= ",")),
"left join 'pad-us_nlcd/nbcdState.csv'.nbcdState b",
"on a.GEOID10 = b.fips\"",
"pad-us_nlcd/nbcdState.shp shp/cusaStatesAea.shp")
system( ogr2ogr)
zip(
"pad-us_nlcd.zip",
list.files(
path= "pad-us_nlcd",
pattern= "^nbcdState\\.[a-z]{3}",
full.names= TRUE))
3.5.5 repeat for counties
~~~~~~~~~~~~~~~~~~~~~~~~~~