Spatial Data in R
Lab: Land Use Land Cover (LULC)
Michael France Nelson
Spring 2022
Lab Objectives and Overview
The goals and key concepts of this lab are:
- Practice working with raster data
- Gain experience dealing with the internal structure of complex R objects.
- Familiarize yourself with the National Land Cover Database data and retrieving its data
- Patch delineation and quantification of patch-level characteristics
Lab Set-Up
We’ll follow parts of the R walkthrough in F+F (with some modifications as noted below) starting with the Examples in R, section 3.3.
Data and Packages
Note that you’ll have to install some packages to follow along with the code.
Check the updated and verified Chapter 3 code in Moodle for the complete list of packages needed.
First, your life will be easier if you create a variable to hold the path pointing to the book’s data:
require(here)
data_dir = here(
"data",
"Fletcher_Fortin-2018-Supporting_Files",
"data")Now you should be able to read in the NLCD example data:
require(raster)
nlcd = raster(file.path(
data_dir, "nlcd2011gv2sr"
))Fletcher and Fortin Section 3.3.3
Now you’re ready to proceed with section 3.3.3 of the book. Consult the updated chapter 3 code in Moodle, as well as the code in the book.
When you read and factorize the data, you should be able to create the following plot with levelplot() from package rasterVis:
Forest Mask
Follow the instructions to create the forest layer raster mask:
Patch Characteristics
You can do the patch delineation using clump as in the book and plot your results:
## Loading required namespace: igraph
The packages SDMTools, which the book uses, is no longer available. Check the updated code script for chapter 3 in Moodle for alternative options.
We’ll be using functions in the landscapemetrics package, so make sure you install it!
landscapemetrics has a lot of functions available for characterizing patch-level features. You can get a list of the available functions with the list_lsm() function:
list_lsm(level = "patch")## metric name type level function_name
## 1 area patch area area and edge metric patch lsm_p_area
## 2 cai core area index core area metric patch lsm_p_cai
## 3 circle related circumscribing circle shape metric patch lsm_p_circle
## 4 contig contiguity index shape metric patch lsm_p_contig
## 5 core core area core area metric patch lsm_p_core
## 6 enn euclidean nearest neighbor distance aggregation metric patch lsm_p_enn
## 7 frac fractal dimension index shape metric patch lsm_p_frac
## 8 gyrate radius of gyration area and edge metric patch lsm_p_gyrate
## 9 ncore number of core areas core area metric patch lsm_p_ncore
## 10 para perimeter-area ratio shape metric patch lsm_p_para
## 11 perim patch perimeter area and edge metric patch lsm_p_perim
## 12 shape shape index shape metric patch lsm_p_shapeLet’s try out the patch area stats using the 4- and 8-neighbor rules:
patches_area_4 = lsm_p_area(nlcd.forest, directions = 4)
patches_area_8 = lsm_p_area(nlcd.forest, directions = 8)You can plot histograms, noting that the area is held in the value column:
par(mfrow = c(1, 2))
hist(
patches_area_8$value,
main = "8-neighbor rule", xlab = "Patch Area")
hist(
patches_area_4$value,
main = "4-neighbor rule", xlab = "Patch Area")
Those look pretty awful, but perhaps doing a log transformation on the areas would be more informative. Code not shown, you should try to re-create these on your own.
It seems that with the 4-neighbor rule, you get many more smaller patches. Can you reason out why this is?
Patch-Level Characteristics: Agricultural Landscapes
Cranberry Bog Agricultural Landscapes
You’ll use the ideas from the F+F walkthrough to contrast agricultural landscapes in two cranberry-producing regions in Massachusetts and Wisconsin.
Navigate to https://www.mrlc.gov/viewer/ to view and download NLCD data.
It may also be useful to look at some aerial imagery of the two regions. Google Earth is a convenient tool (Google Maps is difficult to access on campus).
You should explore the regions to get a feel for what the cranberry bogs look like from above.
Next, you’ll need to retrieve some NLCD data:
For both regions, make sure to download the 2019 Land Cover data
Download NLCD Data
Tomah, WI
In the interactive map, locate the Tomah, WI region. It’s centered around 44.01 degrees North and 90.44 degrees West.
Trace out a square-ish rectangle of about 100km2.
Carver, MA
Repeat the procedure for the region west of Carver, MA., centered near 41.89 degrees North and 70.82 degrees West.
When the data are ready for download, you’ll get an email. Be sure to unzip them into your data subfolder and rename the directories with more descriptive names! I used tomah_bogs and carver_bogs.
You’ll need to look inside of the unzipped directories to get the full names of raster files.
They will be distributed as .tiff files (tagged image format file), a.k.a. GeoTIFF. Note that there will be two different tiff files (with associated metadata files):
- one will be a change index file (you don’t want this one)
- the other will be the “Land_Cover_L48” file. This is the one you want
Note also that the filenames will be long and horrible (that’s why you renamed the containing directories).
To read the Carver, I used the following code on my computer:
require(raster)
carver = raster(here(
"data", "carver_bogs",
"NLCD_2019_Land_Cover_L48_20210604_ahKJ7ZCAqgiHxH5xKCUr.tiff"))NOTE: the exact filenames will be different for your downloaded data.
NLCD Factor Levels
We’ll need to factorize our tomah and carver rasters (see the book example in the walkthrough for how to do this using as.factor())
I’ve created new objects called tomah_f and carver_f to hold the factorized versions of the rasters (code not shown).
Once we’ve done that, we need some R trickery to set the factor level names.
Factor Level Metadata
You may have noticed the NLCD_landcover_legen... csv files. You should open them up and take a look at the contents. They contain the names of the factor levels and the corresponding numeric codes.
We would like to associate the landcover classes with the cell values in the raster, just like F+F did in their walkthrough.
To do so, we can read in the levels and use a trick to keep only the entries that are present in the data.
First, read the file and notice that there are many cell codes (the Value column) that don’t have corresponding level names (the Legend column):
carver_levels = read.csv(
here(
"data", "carver_bogs",
"NLCD_landcover_legend_2018_12_17_ahKJ7ZCAqgiHxH5xKCUr.csv"),
na.strings = "")
head(carver_levels)## Value Legend
## 1 0 Unclassified
## 2 1 <NA>
## 3 2 <NA>
## 4 3 <NA>
## 5 4 <NA>
## 6 5 <NA>When you looked in the csv file, you probably noticed that in the Legend column, there were lots of missing entries. We’d like these to be coded as missing data by R, i.e. NA values
The na.strings = "" argument is the key: it specifies that any empty cells should be read as NA rather than empty strings.
Here’s the trick: we can use the complete.cases() function to remove all the rows that have a NA in the second column.
complete.cases() is quite a handy function for filtering out unused rows in a data.frame.
# Returns a TRUE for each row that has entries in both columns:
complete.cases(carver_levels)## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE
## [18] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [35] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [52] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [69] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE
## [86] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE# Use the Boolean vector to do a logical subsetting operation via an intermediate vector of row indices (created using the which() function:
carver_lvl_rows =
which(complete.cases(carver_levels))Now we’re ready to set the factor levels, just like they did in the book!
levels(carver_f) = carver_levels[carver_lvl_rows, ]## Error in .checkLevels(levels(x)[[1]], value): the first column name of the raster attributes (factors) data.frame should be "ID"Or not…
Take at the column names of carver_levels
names(carver_levels)## [1] "Value" "Legend"You’ll have to change the second column name (code not shown).
Once you’ve made the change, you can run this again:
levels(carver_f) = carver_levels[carver_lvl_rows, ]## Warning in .checkLevels(levels(x)[[1]], value): the number of rows in the raster attributes (factors)
## data.frame is unexpected## Warning in sort(newv[, 1]) == sort(old[, 1]): longer object length is not a multiple of shorter object
## length## Warning in .checkLevels(levels(x)[[1]], value): the values in the "ID" column in the raster attributes
## (factors) data.frame have changedWhen you’ve got your code correct, you should see the correct table of levels:
levels(carver_f)## [[1]]
## ID Legend
## 1 0 Unclassified
## 12 11 Open Water
## 13 12 Perennial Snow/Ice
## 22 21 Developed, Open Space
## 23 22 Developed, Low Intensity
## 24 23 Developed, Medium Intensity
## 25 24 Developed, High Intensity
## 32 31 Barren Land
## 42 41 Deciduous Forest
## 43 42 Evergreen Forest
## 44 43 Mixed Forest
## 53 52 Shrub/Scrub
## 72 71 Herbaceuous
## 82 81 Hay/Pasture
## 83 82 Cultivated Crops
## 91 90 Woody Wetlands
## 96 95 Emergent Herbaceuous WetlandsLevelplot Your Cranberry Bogs
Creating a nice-looking levelplot of the cranberry bog data isn’t as easy as for the nlcd data in the F+F examples.
I’ve included the complete code I used to create this plot of the Tomah region landcover classes.
You should use this code as a template to create a similar plot for the Carver MA data.
Click to show/hide code.
# Tomah
tomah = raster(here(
"data", "tomah_bogs",
"NLCD_2019_Land_Cover_L48_20210604_ZOZXa4B4SSSqU6Qj6lBf.tiff"))
# Factorizing removes the color information from the "levels" attribute.
tomah_f = as.factor(tomah)
# Read the level legend data
tomah_levels = read.csv(
here(
"data", "tomah_bogs",
"NLCD_landcover_legend_2018_12_17_ZOZXa4B4SSSqU6Qj6lBf.csv"),
na.strings = "")
# Set the level names
names(tomah_levels) = c("ID", "Landcover")
# Get the row numbers of the levels used in the data:
tomah_lvl_rows = which(complete.cases(tomah_levels))
# The colortable() function can retrieve the color information distributed with the raster metadata.
# Use the row indices to extract only the colors you need.
col_tomah = colortable(tomah_f)[tomah_lvl_rows]
# set the landcover type levels for Tomah
# Note that the labels must be the second column for levelplot() to find them.
levels(tomah_f) = tomah_levels[tomah_lvl_rows, ]# Create a plot with the original colors and correctly labeled landcover types
require(rasterVis)
levelplot(
tomah_f,
col.regions = col_tomah,
main = "Landcover Classes in Tomah, WI",
scales = list(draw = F))Cranberry Landscapes: Patch Delineation
We want to reclassify the Carver and Tomah landscapes into binary rasters showing agriculture/non-agriculture.
As a first step, examine the levels table again:
levels(carver_f)[[1]]## ID Legend
## 1 0 Unclassified
## 12 11 Open Water
## 13 12 Perennial Snow/Ice
## 22 21 Developed, Open Space
## 23 22 Developed, Low Intensity
## 24 23 Developed, Medium Intensity
## 25 24 Developed, High Intensity
## 32 31 Barren Land
## 42 41 Deciduous Forest
## 43 42 Evergreen Forest
## 44 43 Mixed Forest
## 53 52 Shrub/Scrub
## 72 71 Herbaceuous
## 82 81 Hay/Pasture
## 83 82 Cultivated Crops
## 91 90 Woody Wetlands
## 96 95 Emergent Herbaceuous WetlandsNote that the double brackets are needed because the object returned by levels(carver_f) is a list:
class(levels(carver_f))## [1] "list"and we want the dataframe that it contains:
class(levels(carver_f)[[1]])## [1] "data.frame"It looks like the level of interest is “Cultivated Crops”.
We can use the reclassify function to create our binary raster, we just need a reclassification raster first.
Examine the help entry for reclassify (it’s not too difficult compared to many R help entries). Based on the help entry, ti seems like we want to use the 2-column reclassification matrix option. This might seem a little bit like the raster reclass procedure in Arc!
Per the help entry:
A
2-column matrixrepresents (“is”, “becomes”) which can be useful for integer values. In that case, the right argument is automatically set to NA
In plain English, that means we want the “becomes” column to hold values of zero in all positions corresponding to classes we don’t care about.
This sounds complicated, but the syntax to build the reclassification matrix is easy:
cbind(
levels(carver_f)[[1]][, 1],
levels(carver_f)[[1]][, 2] == "Cultivated Crops")## [,1] [,2]
## [1,] 0 0
## [2,] 11 0
## [3,] 12 0
## [4,] 21 0
## [5,] 22 0
## [6,] 23 0
## [7,] 24 0
## [8,] 31 0
## [9,] 41 0
## [10,] 42 0
## [11,] 43 0
## [12,] 52 0
## [13,] 71 0
## [14,] 81 0
## [15,] 82 1
## [16,] 90 0
## [17,] 95 0You should save the result in a matrix called reclass_cultivated
Now you can use it to reclassify both the Tomah and Carver rasters. You should factor them too and give the levels descriptive names, so that they will work with levelplot()
# tomah_reclass = reclassify(tomah_f, reclass_cultivated)
# carver_reclass = reclassify(carver_f, reclass_cultivated)
tomah_reclass = as.factor(reclassify(tomah_f, reclass_cultivated))
carver_reclass = as.factor(reclassify(carver_f, reclass_cultivated))
levels(carver_reclass)[[1]]$Cover = c("Other", "Cultivated")
levels(tomah_reclass)[[1]]$Cover = c("Other", "Cultivated")When you are successful, your rasters should look something like these:
levelplot(
tomah_reclass,
main = "Mike's plot of Tomah, WI",
col.regions = c("white", "forestgreen"),
# This suppresses the axes:
scales = list(draw = F))
levelplot(
carver_reclass,
main = "Mike's plot of Carver, MA",
col.regions = c("white", "forestgreen"),
# This suppresses the axes:
scales = list(draw = F))
I called my rasters tomah_reclass and carver_reclass
That was a lot of work to get our rasters ready for analysis! The good news is that now we are ready to proceed.
Finally, let’s calculate some patch-level statistics and see how the agricultural patches in the two regions differ.
Remember our friend from landscapemetrics:
list_lsm(level = "patch")## # A tibble: 12 x 5
## metric name type level function_name
## <chr> <chr> <chr> <chr> <chr>
## 1 area patch area area and edge metric patch lsm_p_area
## 2 cai core area index core area metric patch lsm_p_cai
## 3 circle related circumscribing circle shape metric patch lsm_p_circle
## 4 contig contiguity index shape metric patch lsm_p_contig
## 5 core core area core area metric patch lsm_p_core
## 6 enn euclidean nearest neighbor distance aggregation metric patch lsm_p_enn
## 7 frac fractal dimension index shape metric patch lsm_p_frac
## 8 gyrate radius of gyration area and edge metric patch lsm_p_gyrate
## 9 ncore number of core areas core area metric patch lsm_p_ncore
## 10 para perimeter-area ratio shape metric patch lsm_p_para
## 11 perim patch perimeter area and edge metric patch lsm_p_perim
## 12 shape shape index shape metric patch lsm_p_shapeWe can use these functions to characterize the differences between agricultural patches in the two regions
Report
Q1: Carver Levelplot
- Q1 (3 pts.): Include a levelplot map of your Carver landscape. Your plot must show all of the NLCD categories in different colors, and include a legend. Check out the code I shared above for the Tomah landscape for some hints.
Q2: Agricultural Maps
- Q2 (2 pts.): Include maps of your binary agricultural rasters for Tomah and Carver.
Q3: Qualitative Differences
Examine the maps of your Tomah and Carver binary agricultural landscapes.
- Q3 (3 pts.): Qualitatively describe any differences you see. Some things to consider:
- Shape of the individual agricultural patches
- Areas of the individual agricultural patches
- Configuration (spatial arrangement) of the patches and the matrix.
Q4-6: Quantitative Differences
Explore some of the patch statistics for your Tomah and Carver binary agricultural landscapes.
Choose two patch-level metrics.
Hint: use list_lsm(level = "patch") to print a list of available metrics.
- Q4 (2 pts.): List and briefly describe the two landscape metrics you chose.
- Q5 (2 pts.): Include histograms of the two metrics for both landscapes.
- Q6 (4 pts.): Interpret the differences (or similarities) you see in the histograms in the context of the definition of your two chosen landscape metrics.