Landfire Data Acquisition and Summarizing
Michael Campbell
8/23/2021
Acquiring the data
- Navigate to https://landfire.gov/
- Select Get Data
- Select Data Distribution site (DDS)
- Click on the map in the approximate area you plan to download data (this step doesn’t really matter – the DDS has a dynamic, navigable map where you will zoom/pan to your area of interest)
- In the DDS, click the Download Tool button
- Select the layer(s) you want to download. If you want to use the same dataset that I used in the PJ paper, then you would use the Existing Vegetation Type layer (under LF 2016 Remap (LF_200))
- You have a few options for defining the spatial extent within which you will be downloading data from the DDS:
- By drawing a box
- By selecting from a “template” (e.g. a county or state boundary)
- By manually entering in bounding box coordinates (in lat/long decimal degrees)
- In this example, I’m going to draw a box around some random location in the Wasatch
- A new tab/window will open up, where you should select Download
- Another new window will open up, which will remain open until your data are ready to download. Large chunks of data can take some time.
- For reference, I’ll be downloading my data into the following directory structure:
- E:\Mickey\TEMP\LF (main directory)
- E:\Mickey\TEMP\LF\a01_compressed (where the ZIP file will be stored)
- E:\Mickey\TEMP\LF\a02_extracted (where the compressed contents will be extracted to)
Extracting the data
Now we’ll get into some R stuff. This could obviously be done outside of R, but why not make things harder than they have to be? :)
# find zip file
in.dir <- "E:\\Mickey\\TEMP\\LF\\a01_compressed"
out.dir <- "E:\\Mickey\\TEMP\\LF\\a02_extracted"
zip.file <- list.files(in.dir, pattern = "*.zip$", full.names = T)
print(zip.file)## [1] "E:\\Mickey\\TEMP\\LF\\a01_compressed/lf12820698_US_200EVT.zip"# extract contents
unzip(zip.file, exdir = out.dir)
new.files <- list.files(out.dir, recursive = T)
print(new.files)## [1] "US_200EVT/info/arc.dir" "US_200EVT/info/arc0000.dat"
## [3] "US_200EVT/info/arc0000.nit" "US_200EVT/info/arc0001.dat"
## [5] "US_200EVT/info/arc0001.nit" "US_200EVT/info/arc0002.dat"
## [7] "US_200EVT/info/arc0002.nit" "US_200EVT/meta1.html"
## [9] "US_200EVT/Metadata.cpg" "US_200EVT/Metadata.dbf"
## [11] "US_200EVT/Metadata.prj" "US_200EVT/Metadata.sbn"
## [13] "US_200EVT/Metadata.sbx" "US_200EVT/Metadata.shp"
## [15] "US_200EVT/Metadata.shx" "US_200EVT/Metadata.xml"
## [17] "US_200EVT/output_parameters.txt" "US_200EVT/US_200EVT.dbf"
## [19] "US_200EVT/us_200evt/dblbnd.adf" "US_200EVT/us_200evt/hdr.adf"
## [21] "US_200EVT/us_200evt/prj.adf" "US_200EVT/us_200evt/sta.adf"
## [23] "US_200EVT/us_200evt/vat.adf" "US_200EVT/us_200evt/w001001.adf"
## [25] "US_200EVT/us_200evt/w001001x.adf"Loading in the LANDFIRE raster data
Now let’s take a look at the data. First, we’ll read in the raster using the raster library.
# load the library
library(raster)## Warning: package 'raster' was built under R version 4.0.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.0.3# read in and plot the raster data
lf <- raster("E:\\Mickey\\TEMP\\LF\\a02_extracted\\US_200EVT\\us_200evt")
plot(lf)
# take a look at the attributes
head(lf@data@attributes[[1]])## ID COUNT EVT_NAME
## 1 7011 38480 Rocky Mountain Aspen Forest and Woodland
## 2 7012 9801 Rocky Mountain Bigtooth Maple Ravine Woodland
## 3 7016 4162 Colorado Plateau Pinyon-Juniper Woodland
## 4 7019 5 Great Basin Pinyon-Juniper Woodland
## 5 7020 15983 Inter-Mountain Basins Subalpine Limber-Bristlecone Pine Woodland
## 6 7049 487 Rocky Mountain Foothill Limber Pine-Juniper Woodland
## LFRDB EVT_FUEL
## 1 7011 2011
## 2 7012 2012
## 3 7016 2016
## 4 7019 2019
## 5 7020 2020
## 6 7049 2049
## EVT_FUEL_N EVT_LF
## 1 Tr Rocky Mountain Aspen Forest and Woodland Tree
## 2 Tr Rocky Mountain Bigtooth Maple Ravine Woodland Tree
## 3 Tr Colorado Plateau Pinyon-Juniper Woodland Tree
## 4 Tr Great Basin Pinyon-Juniper Woodland Tree
## 5 Tr Inter-Mountain Basins Subalpine Limber-Bristlecone Pine Woodland Tree
## 6 Tr Rocky Mountain Foothill Limber Pine-Juniper Woodland Tree
## EVT_PHYS EVT_GP EVT_GP_N
## 1 Hardwood 602 Aspen Forest, Woodland, and Parkland
## 2 Hardwood 605 Bigtooth Maple Woodland
## 3 Conifer 630 Pinyon-Juniper Woodland
## 4 Conifer 630 Pinyon-Juniper Woodland
## 5 Conifer 621 Limber Pine Woodland
## 6 Conifer 621 Limber Pine Woodland
## SAF_SRM EVT_ORDER EVT_CLASS
## 1 SAF 217: Aspen Tree-dominated Open tree canopy
## 2 SRM 418: Bigtooth Maple Tree-dominated Open tree canopy
## 3 SRM 504: Juniper-Pinyon Pine Woodland Tree-dominated Open tree canopy
## 4 SRM 412: Juniper-Pinyon Woodland Tree-dominated Open tree canopy
## 5 SAF 209: Bristlecone Pine Tree-dominated Open tree canopy
## 6 SAF 219: Limber Pine Tree-dominated Open tree canopy
## EVT_SBCLS R G B RED GREEN BLUE
## 1 Deciduous open tree canopy 192 255 138 0.752941 1.000000 0.541176
## 2 Deciduous open tree canopy 171 255 138 0.670588 1.000000 0.541176
## 3 Evergreen open tree canopy 121 255 94 0.474510 1.000000 0.368627
## 4 Evergreen open tree canopy 78 255 59 0.305882 1.000000 0.231373
## 5 Evergreen open tree canopy 63 255 46 0.247059 1.000000 0.180392
## 6 Evergreen open tree canopy 174 245 224 0.682353 0.960784 0.878431Summarizing the data
As you can see, there are several different attributes that may be of interest. The most detailed level of classification is in the EVT_NAME field. But, by the sounds of it, you might be looking for a coarser-level of classification. Whatever level of specificity you’re looking for, the summarization approach will be nearly the same. The other attribute of interest is COUNT. This represents the number of pixels in the raster dataset that each vegetation type possesses. So, let’s say you wanted proportions of the highest level of classification specificity. Here’s how you might go about that:
# convert attribute table to a data frame
df <- as.data.frame(lf@data@attributes[[1]])
head(df)## ID COUNT EVT_NAME
## 1 7011 38480 Rocky Mountain Aspen Forest and Woodland
## 2 7012 9801 Rocky Mountain Bigtooth Maple Ravine Woodland
## 3 7016 4162 Colorado Plateau Pinyon-Juniper Woodland
## 4 7019 5 Great Basin Pinyon-Juniper Woodland
## 5 7020 15983 Inter-Mountain Basins Subalpine Limber-Bristlecone Pine Woodland
## 6 7049 487 Rocky Mountain Foothill Limber Pine-Juniper Woodland
## LFRDB EVT_FUEL
## 1 7011 2011
## 2 7012 2012
## 3 7016 2016
## 4 7019 2019
## 5 7020 2020
## 6 7049 2049
## EVT_FUEL_N EVT_LF
## 1 Tr Rocky Mountain Aspen Forest and Woodland Tree
## 2 Tr Rocky Mountain Bigtooth Maple Ravine Woodland Tree
## 3 Tr Colorado Plateau Pinyon-Juniper Woodland Tree
## 4 Tr Great Basin Pinyon-Juniper Woodland Tree
## 5 Tr Inter-Mountain Basins Subalpine Limber-Bristlecone Pine Woodland Tree
## 6 Tr Rocky Mountain Foothill Limber Pine-Juniper Woodland Tree
## EVT_PHYS EVT_GP EVT_GP_N
## 1 Hardwood 602 Aspen Forest, Woodland, and Parkland
## 2 Hardwood 605 Bigtooth Maple Woodland
## 3 Conifer 630 Pinyon-Juniper Woodland
## 4 Conifer 630 Pinyon-Juniper Woodland
## 5 Conifer 621 Limber Pine Woodland
## 6 Conifer 621 Limber Pine Woodland
## SAF_SRM EVT_ORDER EVT_CLASS
## 1 SAF 217: Aspen Tree-dominated Open tree canopy
## 2 SRM 418: Bigtooth Maple Tree-dominated Open tree canopy
## 3 SRM 504: Juniper-Pinyon Pine Woodland Tree-dominated Open tree canopy
## 4 SRM 412: Juniper-Pinyon Woodland Tree-dominated Open tree canopy
## 5 SAF 209: Bristlecone Pine Tree-dominated Open tree canopy
## 6 SAF 219: Limber Pine Tree-dominated Open tree canopy
## EVT_SBCLS R G B RED GREEN BLUE
## 1 Deciduous open tree canopy 192 255 138 0.752941 1.000000 0.541176
## 2 Deciduous open tree canopy 171 255 138 0.670588 1.000000 0.541176
## 3 Evergreen open tree canopy 121 255 94 0.474510 1.000000 0.368627
## 4 Evergreen open tree canopy 78 255 59 0.305882 1.000000 0.231373
## 5 Evergreen open tree canopy 63 255 46 0.247059 1.000000 0.180392
## 6 Evergreen open tree canopy 174 245 224 0.682353 0.960784 0.878431# remove unnecessary fields
df <- df[,c("EVT_NAME", "COUNT")]
names(df) <- c("veg.type", "num.pix")
# calculate relative abundance (in percent of entire study area)
df$rel.abund <- 100 * df$num.pix / sum(df$num.pix)
# sort by abundance
df <- df[order(df$rel.abund, decreasing = T),]
head(df)## veg.type
## 10 Rocky Mountain Subalpine Dry-Mesic Spruce-Fir Forest and Woodland
## 1 Rocky Mountain Aspen Forest and Woodland
## 7 Southern Rocky Mountain Dry-Mesic Montane Mixed Conifer Forest and Woodland
## 54 Rocky Mountain Cliff Canyon and Massive Bedrock
## 5 Inter-Mountain Basins Subalpine Limber-Bristlecone Pine Woodland
## 8 Southern Rocky Mountain Mesic Montane Mixed Conifer Forest and Woodland
## num.pix rel.abund
## 10 40448 17.019558
## 1 38480 16.191470
## 7 27133 11.416922
## 54 16336 6.873801
## 5 15983 6.725267
## 8 13813 5.812182So, that’s basically it! Then you could use whatever plotting technique you wanted. Here’s a quick pie chart example:
# remove classes with less than 1% abundance
df <- df[df$rel.abund > 1,]
# plot pie chart
pie(df$rel.abund, labels = df$veg.type)
…Pretty ugly, obviously, but there are many ways to pretty it up.
The only other thing I can think of is you possibly wanting to clip the original raster to the extent of some polygon of interest (e.g. a plot area, or a county, or whatever). In that case, you could use the rgdal library, along with the readOGR() function to read in a polygon shapefile of interest. You would then just use the mask() function within the raster library to clip the raster to the extent of your polyogn prior to summarization.
Let me know if you have any questions!