generated from jtr13/quarto-edav-template
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ch2_a.qmd
502 lines (372 loc) · 26.8 KB
/
ch2_a.qmd
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
---
title-block-banner: true
bibliography: references.bib
---
## Overview
In this chapter we will learn about different types of spatial datasets (`raster` and `vector`). We will visualize these spatial datasets using static and interactive plotting options available in R. We will also explore different color palettes options available for generating spatial maps for scientific/technical reporting of spatial datasets. We will also learn briefly explore coordinate reference systems and map projections for spatial data representation.
### Sample dataset
We will familiarize ourselves with several open-source global datasets and use them to practice spatial mapping and computing in R.
+:-----------------------------------------------------------------------------------------------------------------------------------------+:------------------------------------------:+
| 1. Global surface soil moisture from NASA's Soil Moisture Active Passive (SMAP) satellite [@entekhabi2009] | {width="242"} |
| | |
| <https://smap.jpl.nasa.gov/> | |
+------------------------------------------------------------------------------------------------------------------------------------------+--------------------------------------------+
| 2. Normalized Difference Vegetation Index (NDVI) from Moderate Resolution Imaging Spectroradiometer (MODIS) [@huete1999modis] | {width="250"} |
| | |
| <https://modis.gsfc.nasa.gov/data/dataprod/mod13.php> | |
+------------------------------------------------------------------------------------------------------------------------------------------+--------------------------------------------+
| 3. Global climate reference land regions from Coupled Model Intercomparison Project (CMIP) project [@iturbide2020] | {width="843"} |
| | |
| <https://essd.copernicus.org/articles/12/2959/2020/> | |
+------------------------------------------------------------------------------------------------------------------------------------------+--------------------------------------------+
| 4. Climate classification (Hyper-arid, arid, semi-arid, sub-humid and humid) based on Global Aridity Index Database [@zomer2022version] | {width="497"} |
| | |
| <https://csidotinfo.wordpress.com/2019/01/24/global-aridity-index-and-potential-evapotranspiration-climate-database-v3/> | |
+------------------------------------------------------------------------------------------------------------------------------------------+--------------------------------------------+
### Data download
The sample dataset for this resource is uploaded to `GitHub` for easy access. Download the sample data manually as a zip file from: `https://github.com/Vinit-Sehgal/SampleData` . Once downloaded, extract the zip folder to the current working directory.
Alternatively, use the following script to programmatically download and extract the sample data from the `GitHub` repository.
```{r}
###############################################################
#~~~ Import sample data from GitHub repository
if (dir.exists("SampleData-master")==FALSE){
# First we check if the folder already exists. If not, the download begins
download.file(url = "https://github.com/Vinit-Sehgal/SampleData/archive/master.zip",
destfile = "SampleData-master.zip") # Download ".Zip"
# Unzip the downloaded .zip file
unzip(zipfile = "SampleData-master.zip")
}
# getwd() # Current working directory
list.files("./SampleData-master") # List folder contents. Do you see sample datasets?
```
## Using R as GIS
A `geographic information system`, or `GIS` refers to a platform which can map, analyzes and manipulate **geographically referenced** dataset. A geographically referenced data (or geo-referenced data) is a spatial dataset which can be related to a point on Earth with the help of geographic coordinates. Types of geo-referenced spatial data include: rasters (grids of regularly sized pixels) and vectors (polygons, lines, points).
<br>
<center></center>
<br>
A quick and helpful review of spatial data can be found here: <https://spatialvision.com.au/blog-raster-and-vector-data-in-gis/>
### Plotting raster data: with `terra` and `tmap`
In this section, we plot global raster data of surface (\~5 cm) soil moisture from SMAP. Let's start by first importing the global soil moisture raster.
```{r, message = FALSE, warning = FALSE,results='asis'}
# Import package for raster operations
library(terra)
# Import SMAP soil moisture raster from the downloaded folder
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
```
Once we have imported the `SpatRaster` (short for "spatial raster") using `rast()` function from `terra` package, let's note its attributes. Notice the `dimensions`, `resolution`, `extent`, `crs` i.e. coordinate reference system and `values`. Note that the cell of one raster layer can only hold a single numerical value.
```{r, message = FALSE, warning = FALSE,results='asis'}
# Print raster attributes
print(sm)
# Try:
# dim(sm) # Dimension (nrow, ncol, nlyr) of the raster
# terra::res(sm) # X-Y resolution of the raster
# terra::ext(sm) # Spatial extent of the raster
# terra::crs(sm) # Coordinate reference system
```
Now let's plot the raster using `terra::plot`.
```{r, message = FALSE, warning = FALSE,results='asis'}
# Basic Raster plot
terra::plot(sm, main = "Soil Moisture")
```
### Color palettes
Using a good color palette is an important aspect of spatial mapping. Choice of a good colormap can help the readers understand the key aspects of the map. The selected colors must adequately represent the key features and their differences, wherever applicable, with the least distortion, ambiguity or effort. There are several libraries available in R specifically dedicated to generating color pelettes for scientific mapping. We will also learn the the usage of `cetcolor` and `scico` packages to generate perceptually uniform and color-blindness friendly palettes.
<br>
Some key packages for generating color palettes for scientific mapping are:
1. `RColorBrewer`: <https://colorbrewer2.org/#type=sequential&scheme=BuGn&n=3>
2. `cetcolors` (Perceptually Uniform Colour Maps): <https://cran.r-project.org/web/packages/cetcolor/vignettes/cet_color_schemes.html>
3. `scico` (Scientific colormaps): <https://github.com/thomasp85/scico>
```{r, message = FALSE, warning = FALSE}
# Libraries for generating Colour palettes
library(RColorBrewer)
library(cetcolor)
library(scico)
# To view color palette
library(unikn)
#~~ A) User defined color palette using brewer.pal
mypal1 = RColorBrewer::brewer.pal(10, "Spectral")
# Brewerpal outputs a max of 9-11 colors. So,the pelette may needs expansion.
mypal2= colorRampPalette(mypal1)(20) # Expand pelette to 20 colors
unikn::seecol(mypal2)
# Some more advanced options include opacity and interpolation method
# mypal2= colorRampPalette(mypal1,
# interpolate = c("linear"), # Choose btw linear/spline interpolation
# alpha = 0.8)(20) # Generate 20 colors, opacity set to 0.8
# unikn::seecol(mypal2)
#~~ B) User defined color palette using scico package
mypal1 = scico::scico(20, alpha = 1.0, direction = -1, palette = "vik")
unikn::seecol(mypal1)
#Check: scico_palette_names() for available palettes!
# Try combinations of alpha=0.5, direction =1, and various different color palette
#~~ C) User defined color palette using cetcolor package
mypal2 = cetcolor::cet_pal(20, name = "r2")
unikn::seecol(mypal2)
# Or reverse color pal
mypal2 = rev(cetcolor::cet_pal(20, name = "r2") )
unikn::seecol(mypal2)
```
### Customizing `terra` plot options
There is a long list of customization operations available while plotting rasters in R.
We will start with basic plot from `terra`, and then explore the function by changing different customization options , such as: Try `horizontal=TRUE`, `interpolate=FALSE`, change `xlim=c(-180, 180)` with `asp=1`, or try `legend.shrink=0.4`.
```{r, message = FALSE, warning = FALSE}
sm=rast("./SampleData-master/raster_files/SMAP_SM.tif") # SMAP soil moisture data
terra::plot(sm,
main = "Scientific Plot of Raster",
#Color options
col = mypal2, # User Defined Color Palette
breaks = seq(0, 1, by=0.1), # Sequence from 0-1 with 0.1 increment
colNA = "lightgray", # Color of cells with NA values
# Axis options
axes=TRUE, # Plot axes: TRUE/ FALSE
xlim=c(-180, 180), # X-axis limit
ylim=c(-90, 90), # Y-axis limit
xlab="Longitude", # X-axis label
ylab="Latitue", # Y-axis label
# Legend options
legend=TRUE, # Plot legend: TRUE/ FALSE
# Miscellaneous
mar = c(3.1, 3.1, 2.1, 7.1), # Margins
grid = FALSE # Add grid lines
)
```
### Spatial plotting with `tmap`
```{r, message = FALSE, warning = FALSE}
library(tmap)
# Set tmap mode: Static plot="plot", Interactive plots="view"
tmap_mode("plot")
tmap_SM = tm_shape(sm)+
tm_grid(alpha = 0.2)+ # Transparency of grid
tm_raster(alpha = 0.7, # Transparency of raster plot
palette = mypal2, # Color pellete
style = "pretty", # Select style
title = "Volumetric Soil Moisture")+ # Plot main title
tm_layout(legend.position =
c("left", "bottom"))+ # Placement of legend
tm_xlab("Longitude")+ # x-lab
tm_ylab("Latitude") # y-lab
tmap_SM
```
### Interactive raster visualization aster data with `mapview`
Functionality of `terra` is largely similar to the legacy package `raster` (created by the same developer, Robert Hijmans). The development of `terra` is inspired by computational efficiency in geospatial operations. However, since `terra` is relatively new, and is continually developed, several other packages require conversion of the `SpatRasters` to `rasterLayer` for backwards compatibility.
To convert a `SpatRaster` to `RasterLayer`, use: `sm2=as(sm, "Raster")`
```{r, message = FALSE, warning = FALSE}
library(mapview)
library(raster)
# Convert SpatRaster to raster (from package raster)
sm2=as(sm, "Raster")
mapview(sm2, # RasterLayer
col.regions = mypal2, # Color palette
at=seq(0, 0.8, 0.1) # Breaks
)
```
### Plotting raster data using `tidyterra`
`tidyterra` is a package that add common methods from the tidyverse for SpatRaster and SpatVectors objects created with the `terra` package. It also adds specific `geom_spat*()` functions for plotting rasters with `ggplot2`.
Note on Performance: `tidyterra` is conceived as a user-friendly wrapper of `terra` using the `tidyverse`} methods and verbs. This approach therefore has a cost in terms of performance.
```{r, message = FALSE, warning = FALSE}
library(tidyterra)
library(ggplot2)
ggplot() +
geom_spatraster(data = sm) +
scale_fill_gradientn(colors=mypal2, # Use user-defined colormap
name = "SM", # Name of the colorbar
na.value = "transparent", # transparent NA cells
labels=(c("0", "0.2", "0.4", "0.6", "0.8")), # Labels of colorbar
breaks=seq(0,0.8,by=0.2), # Set breaks of colorbar
limits=c(0,0.8))+
theme_void() # Try other themes: theme_bw(), theme_gray(), theme_minimal()
```
What if we are interested in a particular region and not the entire globe? We can plot the map for a specific extent (CONUS, in this case) by changing the range of `coord_sf`option. We will also use a different theme: `theme_bw`. Try `xlim = c(114,153)` and `ylim = c(-43,-11)`! We will also add state boundaries and coastline to the plot.
```{r, message = FALSE, warning = FALSE}
sm_conus= ggplot() +
geom_spatraster(data = sm) +
scale_fill_gradientn(colors=mypal2, # User-defined colormap
name = "SM", # Name of the colorbar
na.value = "transparent", # transparent NA cells
labels=(c("0", "0.2", "0.4", "0.6", "0.8")), # Labels of colorbar
breaks=seq(0,0.8,by=0.2), # Set breaks of colorbar
limits=c(0,0.8)) +
coord_sf(xlim = c(-125,-67), # Add extent for CONUS
ylim = c(24,50))+
borders("world", # Add global landmass boundaries
colour="gray43", # Fill light-gray color to the landmass
fill="transparent")+ # Transparent background
borders("state", # Add US state borders
colour="gray43", # Use light-gray color
fill="transparent")+ # Use transparent background
theme_bw() # Black & white theme
print(sm_conus)
```
### Plotting vector data
Importing and plotting shapefiles is equally easy in R. We will import the shapefile of the updated global `IPCC climate reference regions` and global `coastline`as Simple Feature (`sf`) Object. Terra can also be used to import shapefiles as vectors using the function `vect`. However, `sf` is more versatile, especially for shapefiles. Notice how the attributes of the `sf` objects resemble an Excel data sheet.
```{r, message = FALSE, warning = FALSE}
library(sf)
library(terra)
##~~~~ Use sf package for shapefile
# Import the shapefile of global IPCC climate reference regions (only for land)
IPCC_shp = read_sf("./SampleData-master/CMIP_land/CMIP_land.shp")
# View attribute table of the shapefile
print(IPCC_shp) # Notice the attributes look like a data frame
##~~~~ Use terra package for shapefile
IPCC_shp = vect("./SampleData-master/CMIP_land/CMIP_land.shp") # Reference regions
print(IPCC_shp)
dim(IPCC_shp) # Notice the dimensions of the shapefile
# Load global coastline shapefile
coastlines = vect("./SampleData-master/ne_10m_coastline/ne_10m_coastline.shp")
print(coastlines)
# Notice the difference in the geometry of coastlines and IPCC_shp
# Alternatively, download global coastlines from the web
# NOTE: May not work if the online server is down
# download.file("https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_coastline.zip?version=4.0.1",destfile = 'ne_110m_coastline.zip')
# # Unzip the downloaded file
# unzip(zipfile = "ne_110m_coastline.zip",exdir = 'ne-coastlines-110m')
# Interactive polygon map with mapview
library(mapview)
mapview(IPCC_shp) # Scroll over the region of interest and find the attributes for the region
```
Subsetting vector data is similar to selecting a row from a data frame.
```{r, message = FALSE, warning = FALSE}
ENA_poly=IPCC_shp[4,] # Subset shapefile for Eastern North-America (ENA)
ENA_poly # Polygon contents - notice it has 4 attributes
# Plot ENA polygon using terra
terra::plot(ENA_poly, main="Polygon for Eastern North-America")
```
Let us see how to combine rasters with overlapping vectors with `terra` . Remember to execute the lines with `add=TRUE` with the base plot.
```{r, message = FALSE, warning = FALSE}
terra::plot(sm, col=mypal2)
terra::plot(IPCC_shp, add=TRUE, col="transparent", border = "black")
terra::plot(coastlines, add=TRUE, col="transparent", border = "black")
```
Lets mark the location of Baton Rouge on this map. We will first make a base plot, over which coastlines and spatial location will be added. Again, remember to run these lines together.
```{r, message = FALSE, warning = FALSE}
#~~~ Add spatial point to shapefile/ raster
#~~ Make base map
terra::plot(IPCC_shp[c(3,4,6,7),], # IPCC land regions for Contiguous US.
col = "lightgray", # Background color
border = "black") # Border color
#~~ Add coastline
terra::plot(coastlines, col= "maroon", add=TRUE)
#~~ Add spatial point to the plot
Long=-91.0; Lat=30.62 # Lat- Long of Baton Rouge, LA
points(cbind(Long,Lat), # Lat-Long as Spatial Points
col="blue", pch=16, cex=1.2) # Shape, size and color of point
```
### Reprojection of rasters using `terra::project`
A coordinate reference system (CRS) is used to relate locations on Earth (which is a 3-D spheroid) to a 2-D projected map using coordinates (for example latitude and longitude). Projected CRSs are usually expressed in Easting and Northing (x and y) values corresponding to long-lat values in Geographic CRS.
A good description of coordinate reference systems and their importance can be found here:
- <https://docs.qgis.org/3.4/en/docs/gentle_gis_introduction/coordinate_reference_systems.html>
- <https://datacarpentry.org/organization-geospatial/03-crs/>
<br>
<center>{width="60%"}</center>
<br>
In R, the coordinate reference systems or `CRS` are commonly specified in `EPSG` (European Petroleum Survey Group) or PROJ4 format (See: <https://epsg.io/>). Few commonly used projection systems and their codes are summarized below:
------------------------------------------------------------------------
| | | |
|---------------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------------------|-------------|
| *Projection system* | *PROJ.4 code* | *EPSG code* |
| NAD83 (North American Datum 1983) | "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +type=crs" | EPSG:4269 |
| Mercator | "+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs" | ESRI:53004 |
| WGS 84 (World Geodetic System 1984) | "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs" | EPSG:3395 |
| WGS 84 / Pseudo-Mercator -- Spherical Mercator ( used in Google Maps, OpenStreetMap) | "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=`@null` +wktext +no_defs +type=crs" \| | EPSG:3857 |
| Robinson (better for plotting global data due to minimal distortion, except for higher latitudes) | "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs" | ESRI:54030 |
------------------------------------------------------------------------
The `SpatRaster` reprojection process is done with `project()` from the `terra` package.
```{r, message = FALSE, warning = FALSE,results='asis'}
# Importing SMAP soil moisture data
sm=rast("./SampleData-master/raster_files/SMAP_SM.tif")
#~~ Projection 1: NAD83 (EPSG: 4269)
sm_proj1 = terra::project(sm, "epsg:4269")
terra::plot(sm_proj1,
main = "NAD83", # Title of the plot
col = mypal2, # Colormap for the plot
axes = FALSE, # Disable axes
box = FALSE, # Disable box around the plots
asp = NA, # No fixed aspect ratio; asp=NA fills plot to window
legend=FALSE) # Disable legend
#~~ Projection 2: World Robinson projection (ESRI:54030)
sm_proj2 = terra::project(sm, "ESRI:54030")
terra::plot(sm_proj2,
main = "Robinson", # Title of the plot
col = mypal2, # Colormap for the plot
axes = FALSE, # Disable axes
box = FALSE, # Disable box around the plots
asp = NA, # No fixed aspect ratio; asp=NA fills plot to window
legend=FALSE) # Disable legend
```
Let us now plot a map of global surface soil moisture, reprojected to Robinson projection. Notice that to add the raster and vector to the plot, we use the functions `geom_spatraster` and `geom_spatvector` respectively. Another package `spData` provides several important datasets for global mapping, including, US states polygons (`us_states`), World country polygons (`world`), global elevation (`elev.tif`), among others (<https://jakubnowosad.com/spData/>). We will use world country polygons (`world`) in our map.
```{r, message = FALSE, warning = FALSE}
library(spData)
# Import global political boundaries from spData package
WorldSHP=terra::vect(spData::world)
# Generate plot
RobinsonPlot <- ggplot() +
geom_spatraster(data = sm)+ # Plot SpatRaster layer
geom_spatvector(data = WorldSHP,
fill = "transparent") + # Add world political map
ggtitle("Robinson Projection") + # Add title
scale_fill_gradientn(colors=mypal2, # Use user-defined colormap
name = "Soil Moisture", # Name of the colorbar
na.value = "transparent",# Set color for NA values
lim=c(0,0.8))+ # Z axis limit
theme_minimal()+ # Select theme. Try 'theme_void'
theme(plot.title = element_text(hjust =0.5), # Place title in the middle of the plot
text = element_text(size = 12))+ # Adjust plot text size for visibility
coord_sf(crs = "ESRI:54030", # Reproject to World Robinson
xlim = c(-152,152)*100000, # Convert x-y limits from decimal Deg. to meter
ylim = c(-55,90)*100000)
print(RobinsonPlot)
# Save high resolution plot to disk
# ggsave(
# "globalSM.png", # Name of the file to be created
# plot = RobinsonPlot, # Plot to be exported
# bg = "white", # Plot background
# height = 6, # Height
# width = 10, # Width
# units = c("in") # Units ("in", "cm", "mm" or "px")
# )
```
Run the commented script above, and view the exported plot saved on the disk.
### Reclassify and reproject \`SpatRasters\`
So far we have used continuous color scales for plotting rasters. Now we will reclassify the raster in discrete domains based on the cell values. This operation is called *RATifying* the raster i.e. adding the Raster Attribute Table (RAT) to the file. This allows use of discrete color scales for plotting.
```{r, message = FALSE, warning = FALSE}
library(spData)
# Import global political boundaries from spData package
WorldSHP=terra::vect(spData::world)
# Import raster file
sm=rast("./SampleData-master/raster_files/SMAP_SM.tif")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Defite category breaks for the raster
brk= seq(0,0.8, by=0.1)
# Add attribute table to the raster
library(dplyr)
sm_rat=dplyr::mutate(sm, # input raster
rat = cut(SMAP_SM, # Raster attribute
breaks = seq(0,0.8, by=0.1), # Breaks for each class
labels = levels(cut(values(sm$SMAP_SM), breaks=brk)) # Labels for each class
))
# levels(cut(values(sm$SMAP_SM), breaks=brk)) is same as typing
# c("(0,0.1]", "(0.1,0.2]", "(0.2,0.3]",
# "(0.3,0.4]", "(0.4,0.5]", "(0.5,0.6]",
# "(0.6,0.7]", "(0.7,0.8]")
plot(sm_rat$rat) # Simple plot with ratified values
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Generate global plot in Robinson projection for the RATifies raster
RobinsonPlot_cat <- ggplot() +
geom_spatraster(data = sm_rat, # Plot SpatRaster layer
aes(fill = rat))+ # aes = the attribute to use for plotting
geom_spatvector(data = WorldSHP,
fill = "transparent") + # Add world political map
ggtitle("Robinson Projection") + # Add title
scale_fill_manual(values = brewer.pal(length(brk)-1,"Spectral"), # Discrete palette
na.value="transparent", # Set transparent fill values for NA cells
name="SM", # Name of the legend bar
labels=levels(cut(values(sm$SMAP_SM), breaks=brk)))+
theme_minimal()+ # Select theme. Try 'theme_void'
theme(plot.title = element_text(hjust =0.5), # Place title in the middle of the plot
text = element_text(size = 12))+ # Adjust plot text size for visibility
coord_sf(crs = "ESRI:54030", # Reproject to World Robinson
xlim = c(-152,152)*100000, # Convert x-y limits from decimal Deg. to meter
ylim = c(-55,90)*100000)
print(RobinsonPlot_cat)
```
**Other resources:**
- More excellent examples on making maps in R can be found here: <https://bookdown.org/nicohahn/making_maps_with_r5/docs/introduction.html>.
- Quintessential resource for reference on charts and plots in R: <https://www.r-graph-gallery.com/index.html>.