This repository has been archived by the owner on Feb 1, 2024. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 36
/
edwebinar_mar19_ornldaac_tutorial.Rmd
565 lines (432 loc) · 26.2 KB
/
edwebinar_mar19_ornldaac_tutorial.Rmd
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
---
title: "Introduction to Geospatial Analysis in R"
author: "Presented by the ORNL DAAC https://daac.ornl.gov"
date: "March 13, 2019"
output:
html_document:
keep_md: yes
number_sections: yes
toc: yes
html_notebook:
number_sections: yes
toc: yes
pdf_document:
toc: yes
subtitle: NASA Earthdata Webinar
editor_options:
chunk_output_type: inline
---
***
<!--------------------------- LOAD THIS FILE INTO RSTUDIO --------------------------->
# Install and Load Packages
> *Functions featured in this section:*
> **rasterOptions {raster}**
> set global options used by the raster package
In addition to the built-in functionality of R, we will use four packages throughout this exercise. Packages are a collection of documentation, functions, and other items that someone has created and compiled for others to use in R. Install the packages, as well as their dependencies, using the function `install.packages()`.
```{r message=FALSE, warning=FALSE, eval=FALSE}
install.packages("raster", dependencies = TRUE)
install.packages("rgdal", dependencies = TRUE)
install.packages("sf", dependencies = TRUE)
install.packages("tigris", dependencies = TRUE)
```
Most functions we will use are from the raster package or are included upon installation of R. Notice that we can set options for the raster package with `rasterOptions()`. These will help you see how long your code will take to run and help manage large objects.
```{r message=FALSE, warning=FALSE, results="hide"}
library(raster)
rasterOptions(progress = "text") # show the progress of running commands
rasterOptions(maxmemory = 1e+09) # increase memory allowance
rasterOptions(tmpdir = "temp_files") # folder for temporary storage of large objects
library(rgdal)
library(sf)
library(tigris) # provides states() function
```
For package details try `help()` (e.g., `help("raster")`), and to view the necessary arguments of a function try `args()` (e.g., `args(cover)`).
# Load Data
> *Functions featured in this section:*
> **raster {raster}**
> creates a raster object
> **states {tigris}**
> downloads a shapefile of the United States that will be loaded as a SpatialPolygonsDataFrame object
Two GeoTiff files are needed to complete this tutorial, both from the dataset titled "CMS: Forest Carbon Stocks, Emissions, and Net Flux for the Conterminous US: 2005-2010" and freely available through the ORNL DAAC integrated web platform. The dataset provides maps of estimated carbon emissions in forests of the conterminous United States for the years 2006-2010. We will use the maps of carbon emissions caused by fire (GrossEmissions_v101_USA_Fire.tif) and insect damage (GrossEmissions_v101_USA_Insect.tif). These maps are provided at 100 meter spatial resolution in GeoTIFF format using Albers North America projection. Refer to the accompanying "README.md" for instructions on how to download the data.
To begin, be sure to set your working directory using `setwd()` and the filepath to where you saved the data (we use the folder "./data/").
With the `raster()` function, load "GrossEmissions_v101_USA_Fire.tif" and name it *fire* then load "GrossEmissions_v101_USA_Insect.tif" and name it *insect*. The contents of these two files are stored as raster objects. *fire* and *insect* are the primary recipients of our manipulations throughout this exercise.
The function `states()` downloads a shapefile of the United States from the United States Census Bureau. Name the shapefile *myStates*, and it will be stored as a simple feature data frame object.
```{r message=FALSE, warning=FALSE, results="hide"}
fire <- raster("./data/GrossEmissions_v101_USA_Fire.tif")
insect <- raster("./data/GrossEmissions_v101_USA_Insect.tif")
myStates <- states(cb = TRUE) # will download a generalized (1:500k) file
```
# Check the Coordinate Reference System and Plot a Raster
> *Functions featured in this section:*
> **crs {raster}**
> gets the coordinate reference system of a raster object
Use `print()` to view details about the internal data structure of the raster object we named *fire*.
```{r message=FALSE, warning=FALSE}
print(fire)
```
The output lists important attributes of *fire*, like its dimensions, resolution, spatial extent, coordinate reference system, and the minimum and maximum values of the cells (i.e., carbon emissions).
```{r message=FALSE, warning=FALSE}
fire@crs
```
The above command retrieves the coordinate reference system (CRS) of *fire*. Notice the PROJ.4 representation. The first argument of the is "+proj=" and defines the projection. "aea" refers to the NAD83 / Albers NorthAm projection (also shown following "PROJCRS"), and "+units=m" tells us that the resolution of the raster object is in meters. Refer to the attributes of *fire* provided by `print()`. The resolution of the raster is "100, 100 (x, y)" meaning that each cell is 100 meters by 100 meters.
Use the `plot()` function to make a simple image of *fire* and visualize the carbon emissions from fire damage across the forests of the conterminous United States between 2006 and 2010. According to the documentation for the dataset, gross carbon emissions were measured in megagrams of carbon per year per cell.
```{r message=FALSE, warning=FALSE}
plot(fire,
main = "Gross Carbon Emissions from Fire Damage\n across CONUS Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = "Mg C/yr\n", side = 3),
colNA = "black",
box = FALSE)
```
The spatial extent of the raster object is displayed on the x- and y-axes. All NA cells (i.e., cells that have no values) are colored black for better visualization of fire damage. The legend offers the range of cell values and represents them using a default color theme.
Let's examine the raster object we named *insect*. `crs()` retrieves the CRS arguments for *insect* as a Vector object. We use `identical()` to determine if *fire* and *insect* have the same CRS.
```{r message=FALSE, warning=FALSE}
identical(crs(fire), crs(insect))
```
The CRS for the two raster objects are identical.
Plot *insect* but change the content for the argument "main = ", which defines the main title of the plot.
```{r message=FALSE, warning=FALSE}
plot(insect,
main = "Gross Carbon Emissions from Insect Damage\n across CONUS Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = "Mg C/yr\n", side = 3),
colNA = "black",
box = FALSE)
```
You can likely imagine an outline of the United States given the spatial data distribution of the two raster objects.
# Select Data Within a Region of Interest
> *Functions featured in this section:*
> **CRS {rgdal}**
> creates a CRS object using PROJ.4 arguments
> **st_transform {sf}**
> provides re-projection given a CRS
> **st_bbox {sf}**
> provides bounding of a simple feature
> **crop {raster}**
> returns a geographic subset of an object as specified by an Extent object
> **mask {raster}**
> creates a new raster object with the same values as the input object, except for the cells that are NA in the second object
Next, we reduce the size of *fire* and *insect* by choosing a smaller extent of the raster objects. Use `print()` to view details about the internal data structure of the simple feature we named *myStates*.
```{r message=FALSE, warning=FALSE}
print(myStates)
```
*myStates* has 56 rows (features, i.e., polygons) and ten columns (variables or features).
For this exercise, we will focus on carbon emissions for the states Idaho, Montana, and Wyoming. We can use column referencing and indexing to select all column information contained in *myStates*, but for only three rows (polygons). Name the resultant simple feature *threeStates*.
```{r message=FALSE, warning=FALSE}
threeStates <- myStates[myStates$NAME == "Idaho" |
myStates$NAME == "Montana" |
myStates$NAME == "Wyoming", ]
print(threeStates)
```
*threeStates* has only three rows, but the same number of columns as *myStates*.
What does *threeStates* look like plotted? We'll plot the geometry (i.e., the 10th column) of *threeStates* so we only see the outline.
```{r message=FALSE, warning=FALSE}
plot(threeStates$geometry)
```
We can get the *fire* and *insect* data that occurs "within" *threeStates*. First, we must confirm that the three objects share a CRS before we can "match" them on a coordinate plane.
```{r message=FALSE, warning=FALSE}
identical(crs(fire), crs(threeStates))
```
*threeStates* does not have the same CRS as *fire*, so we will make a simple feature object with the projection of *fire* using `st_transform()`. We also use `CRS()` to properly format the projection arguments of *fire*.
```{r message=FALSE, warning=FALSE}
transStates <- st_transform(threeStates, CRS(proj4string(fire)))
plot(transStates$geometry)
```
Plotting the geometry of the new object *transStates* shows that the projection has changed. Notice how the orientation of the polygons has shifted to match the NAD83 / Albers NorthAm projection.
Now that our objects share a CRS, we will compare the extent of *fire* and *transStates*. For the simple feature *transStates* we will use `st_bbox()` to view the bounding (i.e., bounding box) of the object.
```{r message=FALSE, warning=FALSE}
cat("fire extent\n"); fire@extent; cat("transStates extent\n"); st_bbox(transStates)
```
*fire* has a much larger extent than *transStates*.
We will use the `crop()` function to reduce the extent of the two raster objects. Cropping will create a geographic subset of *fire* and *insect* as specified by the extent of *transStates*. We will name the new raster objects to reflect this manipulation.
```{r message=FALSE, warning=FALSE, results="hide"}
# this will take a minute to run
cropFire <- crop(fire, transStates) # crop(raster object, extent object)
cropInsect <- crop(insect, transStates)
```
Now when we plot *cropFire* and *cropInsect*, we will also plot *transStates* "on top" to envision how carbon emissions are distributed across the three states.
```{r message=FALSE, warning=FALSE}
plot(cropFire,
main = "Gross Carbon Emissions from Fire Damage\n across ID, MT, WY Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = "Mg C/yr\n", side = 3),
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
plot(cropInsect,
main = "Gross Carbon Emissions from Insect Damage\n across ID, MT, WY Forests (2005-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = "Mg C/yr\n", side = 3),
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
```
If you look closely at the cells "outside" the boundary of the *transStates* polygons, you can still see cells values. That's because `crop()` changed the extent of the two raster objects to match that of the simple feature object, but the boundary of the *transStates* polygons are rotated to fit the NAD83 / Albers NorthAm projection and does not extend to the entire rectangular extent of the raster objects.
To remove those extraneous cell values, use the `mask()` function to create two new rasters, one for fire damage and one for insect damage. **Note:** You can use `mask()` or `crop()` in either order.
```{r message=FALSE, warning=FALSE, results="hide"}
# this will take a couple of minutes to run
maskFire <- mask(cropFire, transStates) # mask(raster object, mask object)
maskInsect <- mask(cropInsect, transStates)
```
Plot *maskFire* and *maskInsect*.
```{r message=FALSE, warning=FALSE}
plot(maskFire,
main = "Gross Carbon Emissions from Fire Damage\n across ID, MT, WY Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = "Mg C/yr\n", side = 3),
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
plot(maskInsect,
main = "Gross Carbon Emissions from Insect Damage\n across ID, MT, WY Forests (2005-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = "Mg C/yr\n", side = 3),
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
```
These plots demonstrate that the extraneous cells has been removed from outside the boundary of the *transStates* polygons.
# Examine Raster Value Summaries
> *Functions featured in this fection:*
> **extract {raster}**
> extracts values from a raster object at the locations of other spatial data
In this section, we will compare the three states by their carbon emissions from fire damage only.
We will use the `extract()` function to collect the cell values of *maskFire* where the *transStates* simple feature object overlaps the raster object on their shared coordinate reference system. We will use `summary()` to examine the distribution of cell values that we collect.
```{r message=FALSE, warning=FALSE}
# this can take up to an hour to run, so I will load a saved copy for the demonstration
if(file.exists("./data/val_fireStates.Rds")) {
val_fireStates <- readRDS("./data/val_fireStates.rds")
summary(val_fireStates)
}else{
val_fireStates <- extract(maskFire, transStates, df = TRUE) # extract(raster object, extent object)
summary(val_fireStates)
}
```
There are two columns for *val_fireStates*. One is ID, which corresponds with the three states; 1 = Idaho, 2 = Montana, and 3 = Wyoming. The second column is a summary of all cell values across those three states. On average, 56 megagrams of carbon per year are a result of forest destruction by fire damage for all states combined.
To look at the summary for cell values by state, we will use `subset()` to split the data frame into three. In the code below, we subest *val_fireStates* so that only the rows with a "1" for the ID number will be returned. We name the new object with the prefix "temp".
```{r message=FALSE, warning=FALSE}
temp_val_id <- subset(val_fireStates, subset = ID %in% 1)
summary(temp_val_id)
```
The summary demonstrates that there is now only a single value is included in the ID column, and that the distribution of cell values has changed. This resultant data frame object is quite large and has more information than we need. We need only the second column and we don't care for the large number of NA's.
We will use the functions `which()` and `is.na()` to make a new object from the temporary one. We tell R that we want only the second column and the rows of *temp_val_id* that are not NA.
```{r message=FALSE, warning=FALSE}
val_id <- temp_val_id[which(!is.na(temp_val_id$GrossEmissions_v101_USA_Fire)), 2]
summary(val_id)
```
The resultant object, *val_id*, is a vector object (a single column of numbers) with no NA's.
We will do the same with *val_fire* for the states Montana and Wyoming.
```{r message=FALSE, warning=FALSE}
temp_val_mt <- subset(val_fireStates, subset = ID %in% 2)
val_mt <- temp_val_mt[which(!is.na(temp_val_mt$GrossEmissions_v101_USA_Fire)), 2]
temp_val_wy <- subset(val_fireStates, subset = ID %in% 3)
val_wy <- temp_val_wy[which(!is.na(temp_val_wy$GrossEmissions_v101_USA_Fire)), 2]
```
What's the average and range of values for carbon emissions from fire damage within each state for the period 2006 to 2010?
```{r message=FALSE, warning=FALSE}
cat("Idaho\n"); summary(val_id); cat("Montana\n"); summary(val_mt); cat("Wyoming\n"); summary(val_wy)
```
On average, Montana has the highest carbon emissions, but the maximum gross carbon emissions from a single cell occurred in Idaho.
In addition to using `summary()`, we can create graphs to visualize carbon emissions from fire damage within each of the three states. The function `hist()` plots the frequency of cell values. We will set some arguments of the plot so that we can compare carbon emissions across all three states.
```{r message=FALSE, warning=FALSE}
par(mfrow=c(2,2))
hist(val_id,
main = "Idaho",
ylab = "number of cells",
xlab = "megagrams of carbon per year (Mg C/yr)",
ylim = c(0, 120000), # same y-axis limit for all three states
xlim = c(0, 350)) # same x-axis limit for all three states
hist(val_mt,
main = "Montana",
ylab = "number of cells",
xlab = "megagrams of carbon per year (Mg C/yr)",
ylim = c(0, 120000),
xlim = c(0, 350))
hist(val_wy,
main = "Wyoming",
ylab = "number of cells",
xlab = "megagrams of carbon per year (Mg C/yr)",
ylim = c(0, 120000),
xlim = c(0, 350))
```
The histogram shows the number of times (on the y-axis) each unique cell value (on the x-axis) occurs in each state. In other words, it illustrates the variation in carbon emissions from fire damage within the three different states.
# Reclassify Raster Values
> *Functions featured in this section:*
> **reclassify {raster}**
> reclassifies groups of values of a raster object to other values
> **calc {raster}**
> calculates values for a new raster object from another raster object using a formula
Now we are going to change the values of our two raster objects using different methods.
Beginning with *maskFire*, we will use the `calc()` function to code all cells that have fire damage to be two. To use `calc()`, we must define a function that will detect certain cell values and change them to other values.
```{r message=FALSE, warning=FALSE, results="hide"}
reclassFire <- calc(maskFire,
fun = function(x) {
x[x > 0] <- 2
return(x) })
```
The function we defined changed all *maskFire* cell values that were greater than zero to be two.
Check that our reclassification of *maskFire* worked as expected using `summary()`
```{r message=FALSE, warning=FALSE}
summary(reclassFire[])
```
Yes, all values are two or NA.
All cell values of *reclassFire* should be at the same locations as *maskFire* but with a single value.
```{r message=FALSE, warning=FALSE}
plot(reclassFire,
main = "Locations of Forest Disturbance from Fire Damage\n across ID, MT, WY Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend = FALSE,
col = "red",
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
```
The plot of *reclassFire* now illustrates locations where there were carbon emissions owing to fire damaging the forest. Notice that we chose a single color to represent the presence of values using the argument "col = "red"".
Now we will reclassify all values of *maskInsect* that are greater than zero to be one, but instead of using `calc()`, we will use the `reclassify()` function. `reclassify()` uses a matrix to identify the target cell values and to what value those cells will change.
```{r message=FALSE, warning=FALSE, results="hide"}
reclassInsect <- reclassify(maskInsect,
rcl = matrix(data = c(1, 285, 1), # c(from value, to value, becomes)
nrow = 1, ncol = 3))
```
The argument following "rcl =" tells R that values from two to 285 should be reclassified as one. Essentially, we are making the presence of insect damage equal one.
Check the reclassification of *maskInsect* using `summary()`.
```{r message=FALSE, warning=FALSE}
summary(reclassInsect[])
```
All values are one or NA.
Plot *reclassInsect*. All the cell values should be at the same locations as *maskInsect* but will all be the value one.
```{r message=FALSE, warning=FALSE}
plot(reclassInsect,
main = "Locations of Forest Disturbance from Insect Damage\n across ID, MT, WY Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend = FALSE,
col = "dark green",
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
```
The plot illustrates locations where there were carbon emissions owing to insect damaging the forest, so now the information conveyed by the *maskInsect* raster object is presence or absence of insect damage.
# Combine Two Rasters
> *Functions featured in this section:*
> **cover {raster}**
> replaces NA values in the first raster object with the values of the second
Next, we will join *reclassFire* and *reclassInsect* to form a single raster object. According to the documentation for this dataset, there are no overlapping, non-NA cells between the two raster objects. That is, if you were to combine the two rasters object, a cell could take only the value provided by *reclassFire* (i.e., two) or *reclassInsect* (i.e., one), or be NA. This allows us to use the `cover()` function to combine objects. `cover()` is unique because it will replace NA values of *reclassFire* with non-NA values of *reclassInsect*.
```{r message=FALSE, warning=FALSE, results="hide"}
# this will take a couple of minutes to run
fireInsect <- cover(reclassFire, reclassInsect)
```
Check the combination of *reclassFire* and *reclassInsect* using `summary()`.
```{r message=FALSE, warning=FALSE}
summary(fireInsect[])
```
The data distribution of the new raster object shows that the minimum value is now one (i.e., the insect damage value we specified during reclassification) and the maximum value is two (i.e., the fire damage value).
The plotting arguments below now reflect the "breaks" in the values we would like to see illustrated on the plot. Insect damage is displayed as green cells and fire damage as red.
```{r message=FALSE, warning=FALSE}
plot(fireInsect,
main = "Locations of Forest Disturbance\n across ID, MT, WY Forests (2006-2010)",
xlab = "horizontal extent (m)",
ylab = "vertical extent (m)",
legend.args = list(text = " Disturbance\n", side = 3),
breaks = c(0, 1, 2),
col = c("dark green", "red"),
axis.args = list(at = c(0.5, 1.5), labels = c("insect", "fire")),
colNA = "black",
box = FALSE)
plot(transStates$geometry,
border = "white",
add = TRUE)
```
# Reproject and Write a Raster
> *Functions featured in this section:*
> **projectRaster {raster}**
> projects the values of a raster object to a new one with a different projection
> **writeRaster {raster}**
> writes an entire raster object to a file
Reprojecting a raster in R is different than transforming the CRS as we did with the simple feature earlier in the exercise. To reproject a raster we use the `projectRaster()` function and the `CRS()` function to correctly format the projection information.
```{r message=FALSE, warning=FALSE, results="hide"}
# this will take several minutes to run
prjFireInsect <- projectRaster(fireInsect,
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
```
Now, check what the raster object we made using `print()`.
```{r message=FALSE, warning=FALSE}
print(prjFireInsect)
```
It's a new raster object named *prjFireInsect* that has the standard Geographic projection with latitude and longitude expressed in decimal degrees (DD) as its CRS.
We will plot *prjFireInsect* with slightly different arguments than *fireInsect* to "zoom in" to the center of the plot. Also, we will use *threeStates* instead of *transStates* because *threeStates* also uses the Geographic projection.
```{r message=FALSE, warning=FALSE}
plot(prjFireInsect,
main = "Locations of Forest Disturbance\n across ID, MT, WY Forests (2006-2010)",
xlab = "longitude (DD)",
ylab = "latitude (DD)",
legend.args = list(text = " Disturbance\n", side = 3),
las = 1,
ext = prjFireInsect@extent/1.25,
breaks = c(0, 1, 2),
col = c("dark green", "red"),
axis.args = list(at = c(0.5, 1.5), labels = c("insect", "fire")),
box = FALSE)
plot(threeStates$geometry,
border = "black",
add = TRUE)
```
Let's use the `writeRaster()` function to save *prjFireInsect* to the data directory. We will save the file in \*.tif format so that the geographic information of the raster object is retrievable outside of R.
```{r message=FALSE, warning=FALSE, results="hide"}
writeRaster(prjFireInsect, filename = "./data/prjFireInsect.tif", overwrite=TRUE)
```
Use the function `file.exists()`, which tests for the existence of a given file, to ensure that *prjFireInsect* was successfully saved to our working directory.
```{r message=FALSE, warning=FALSE}
file.exists("./data/prjFireInsect.tif")
```
Now we are able to share the raster with others or open it in another program.
# Export a Plot as PNG and Raster as KML
> *Functions featured in this section:*
> **KML {raster}**
> exports raster object data to a KML file
To save the final plot, we use `png()`. This function will open a graphics device that will save the plot we run in \*.png format. We will use the function `dev.off()` to tell R when we are finished plotting and want to close the graphics device.
```{r message=FALSE, warning=FALSE, results="hide"}
png("prjFireInsect.png", width = 800, res = 80)
plot(prjFireInsect,
main = "Locations of Forest Disturbance\n across ID, MT, WY Forests (2006-2010)",
xlab = "longitude (DD)",
ylab = "latitude (DD)",
legend.args = list(text = " Disturbance\n", side = 3),
las = 1,
ext = prjFireInsect@extent/1.25,
breaks = c(0, 1, 2),
col = c("dark green", "red"),
axis.args = list(at = c(0.5, 1.5), labels = c("insect", "fire")),
box = FALSE)
plot(threeStates$geometry,
border = "black",
add = TRUE)
dev.off()
```
Let's also save *prjFireInsect* in \*.kmz format. KML stands for Keyhole Markup Language, and KMZ is the compressed format. These formats were developed for geographic visualization in Google Earth.
```{r message=FALSE, warning=FALSE, results="hide"}
KML(prjFireInsect, "./data/prjFireInsect.kmz", col = c("dark green", "red"), overwrite=TRUE)
```
We successfully saved the raster object as a KML file.
***
This is the end to the tutorial. If you liked this tutorial, please tell us on [EarthData Forum](https://forum.earthdata.nasa.gov/). If you would like to make a suggestion for a new tutorial, please email [email protected].
There is a supplemental document included on GitHub that offers two additional sections, *Perform a Focal Analysis* and *Get Cell Coordinates*.
<!--------------------------------- END OF TUTORIAL --------------------------------->