GIS with R how to start?

class: center, middle, inverse, title-slide

# GIS with R <br/> how to start?
## 🗺
### Jakub Nowosad <br/> (<a href="https://nowosad.github.io" class="uri">https://nowosad.github.io</a>)
### GIS Learning Community, Cincinnati, 2017-11-01

---

### Prerequisites

- The latest versions of [R](https://cloud.r-project.org/) and [RStudio](https://www.rstudio.com/)
- A set of packages installed by running the below code in R:

```r
install.packages(c("raster", "sf", "rasterVis", "tmap", "leaflet",
                   "widgetframe", "tidyverse", "spData", "geofacet",
                   "tidycensus", "linemap", "gapminder", "tigris", 
                   "osmdata", "mapview", "opencage"))
```

- Datasets downloaded from [this webpage](https://github.com/Nowosad/gis_with_r_how_to_start/archive/master.zip).

---
### Prerequisites

Note: `tidycensus` and `opencage` require API keys.

#### To get a `tidycensus` API key

Go to http://api.census.gov/data/key_signup.html and complete the form there.
Replace text YOUR API KEY GOES HERE with your census API key. 
Execute the code below:

```r
library(tidycensus)
census_api_key("YOUR API KEY GOES HERE", install = TRUE)
```

#### To get a `opencage` API key

Go to https://geocoder.opencagedata.com/pricing, select the Free Trial, and complete the form there.
Replace text YOUR API KEY GOES HERE with your opencage API key. 
Execute the code below:

```r
cat("OPENCAGE_KEY=YOUR API CODE HERE", file=file.path(normalizePath("~/"), ".Renviron"), append=TRUE)
```

---
### Prerequisites

#### Check if the keys are installed

Restart R/RStudio (Session -> Restart R) and execute the code below:

```r
Sys.getenv("CENSUS_API_KEY")
Sys.getenv("OPENCAGE_KEY")
```

You should see your keys print to screen.

More info can be found [here](https://walkerke.github.io/tidycensus/articles/basic-usage.html) and [here](https://github.com/ropensci/opencage).

---
class: center, middle

# Part I - Examples

---
### Static visualization (I)

---
### Static visualization (II)

---
### Static visualization (III)

---
### Animation

---
### Interactive visualization

---
### Remote sensing calculations

---
### Terrain calculations

---
### Geocoding

What's the coordinates of New York, Los Angeles, Chicago, Houston, and Philadelphia?

---
class: center, middle
# Part II - An intro to R

---
### How to get R

- https://www.r-project.org/
- http://cran.rstudio.com/bin/linux/
- http://cran.rstudio.com/bin/windows/base/
- http://cran.rstudio.com/bin/macosx/

![](figs/R.png)

---
### Rstudio

- http://www.rstudio.com/ide/download/desktop
- This is an Integrated Development Environment (IDE) for R
- RStudio has many useful features, such as text editor, syntax highlighting, suggested code auto-completion, and many more

![](figs/rstudio.png)

---
### How to get help?

```r
# if you know a function name
?crop
```

- You can also look for help using the help window or **F1** key

## Online help
- [stackoverflow.com](http://stackoverflow.com/questions/tagged/r) - R related questions
- [gis.stackexchange.com](https://gis.stackexchange.com/questions/tagged/r) - GIS and R related questions
- [twitter](http://twitter.com/) - `#rstats`
- Web search engines [Rseek](http://www.rseek.org/), [Duckduckgo](http://duckduckgo.com/), [Google](http://google.com/), [Bing](http://bing.com/), etc.

---
### How to start with R?

- [Intro to R](https://github.com/Nowosad/Intro_to_R)
- [Beginner's guide to R: Introduction](https://www.computerworld.com/article/2497143/business-intelligence/business-intelligence-beginner-s-guide-to-r-introduction.html)
- [R for cats](https://rforcats.net/)
- [Introducing R to a non-programmer in one hour](http://alyssafrazee.com/introducing-R.html)
- [DataCamp: Introduction to R](https://www.datacamp.com/courses/free-introduction-to-r)
- [try R](http://tryr.codeschool.com/)
- [R news and tutorials contributed by R bloggers](https://www.r-bloggers.com/)
- [RStudio Cheat Sheets](https://www.rstudio.com/resources/cheatsheets/)
- [R for Data Science](http://r4ds.had.co.nz/)
- [Efficient R programming](https://csgillespie.github.io/efficientR/)
- [60+ R resources to improve your data skills](http://www.computerworld.com/article/2497464/business-intelligence/business-intelligence-60-r-resources-to-improve-your-data-skills.html)

---
### R packages

- A package is a group of functions
- `install.packages()` can be used to install packages from CRAN:

```r
install.packages("raster")
```

- To use a package, load it with function `library()`
- Unlike `install.packages()`, you need to load selected packages every time you run R!

```r
library(raster)
```

---
class: center, middle
# Part III - An intro to spatial R

---
### Spatial R

<table>
<caption>Differences in emphasis between the fields of Geographic Information Systems (GIS) and Geographic Data Science (GDS)</caption>
 <thead>
  <tr>
   <th style="text-align:left;"> Attribute </th>
   <th style="text-align:left;"> GIS </th>
   <th style="text-align:left;"> GDS </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Home disciplines </td>
   <td style="text-align:left;"> Geography </td>
   <td style="text-align:left;"> Geography, Computing, Statistics </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Software focus </td>
   <td style="text-align:left;"> Graphical User Interface </td>
   <td style="text-align:left;"> Code </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Reproduciblility </td>
   <td style="text-align:left;"> Minimal </td>
   <td style="text-align:left;"> Maximal </td>
  </tr>
</tbody>
</table>

- **sf** and **sp** are the most important R packages to handle vector data; **sf** is a successor of **sp** but its still evolving. Moreover, many other R packages depend on the functions and classes for the **sp** package
- **raster** is an extension of spatial data classes to work with rasters
- There are many ways to visualize spatial data in R, such as the **ggplot2**, **rasterVis**, **tmap**, **leaflet**, and **mapview** packages
- It is easy to connect R with a GIS software - GRASS GIS (**rgrass7**), SAGA (**RSAGA**), QGIS (**RQGIS** and **qgisremote**), and ArcGIS (**arcgisbinding**)

---
### The **sf** package

The **sf** package in an R implementation of [Simple Features](https://en.wikipedia.org/wiki/Simple_Features). This package incorporates:
- A new spatial data class system in R
- Functions for reading and writing data
- Tools for spatial operations on vectors

Most of the functions in this package starts with prefix `st_`.

```r
devtools::install_github("edzer/sfr") # development version
```

```r
install.packages("sf") # stable version
```

You need a recent version of the GDAL, GEOS, Proj.4, and UDUNITS libraries installed for this to work on Mac and Linux. More information on that at https://github.com/edzer/sfr.

```r
library(sf)
```

---
### Reading spatial data

```r
wrld = st_read("data/wrld.shp")
```

```
## Reading layer `wrld' from data source `/home/jakub/Science/misc/gis_with_r_how_to_start/data/wrld.shp' using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
```

```r
ham = st_read("data/hamilton_county.gpkg")
```

```
## Reading layer `hamilton_county.gpkg' from data source `/home/jakub/Science/misc/gis_with_r_how_to_start/data/hamilton_county.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 17 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -84.8203 ymin: 39.02153 xmax: -84.25651 ymax: 39.31206
## epsg (SRID):    4269
## proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
```

---
### Reading spatial data - text files

```r
ham_cities_sf = st_read("data/hamiltion_cities.csv", 
                        options = c("X_POSSIBLE_NAMES=X",
                                    "Y_POSSIBLE_NAMES=Y",
                                    "KEEP_GEOM_COLUMNS=NO"))
```

```
## options:        X_POSSIBLE_NAMES=X Y_POSSIBLE_NAMES=Y KEEP_GEOM_COLUMNS=NO 
## Reading layer `hamiltion_cities' from data source `/home/jakub/Science/misc/gis_with_r_how_to_start/data/hamiltion_cities.csv' using driver `CSV'
## Simple feature collection with 19 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -84.81995 ymin: 39.07862 xmax: -84.26384 ymax: 39.29034
## epsg (SRID):    NA
## proj4string:    NA
```

---
### Writing spatial data

```r
st_write(wrld, "data/new_wrld.shp")
```

```
## Writing layer `new_wrld' to data source `data/new_wrld.shp' using driver `ESRI Shapefile'
## features:       177
## fields:         10
## geometry type:  Multi Polygon
```

```r
st_write(wrld, "data/new_wrld.gpkg")
```

```
## Writing layer `new_wrld' to data source `data/new_wrld.gpkg' using driver `GPKG'
## features:       177
## fields:         10
## geometry type:  Multi Polygon
```

---
### **sf** structure

The **sf** object is a table (`data.frame`) with an additional geometry column - typically named `geom` or `geometry` and spatial metadata (`geometry type`, `dimension`, `bbox`, `epsg (SRID)`, `proj4string`).

```
## Simple feature collection with 5 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   iso_a2            name_long     continent                       geometry
## 1     AF          Afghanistan          Asia MULTIPOLYGON (((61.21081709...
## 2     AO               Angola        Africa MULTIPOLYGON (((16.32652835...
## 3     AL              Albania        Europe MULTIPOLYGON (((20.59024743...
## 4     AE United Arab Emirates          Asia MULTIPOLYGON (((51.57951867...
## 5     AR            Argentina South America MULTIPOLYGON (((-65.5 -55.2...
```

---
### Attributes - the **dplyr** package

- It is easy to use the **dplyr** package on `sf` objects:

```r
library(dplyr)
```

- `filter()`:

```r
wrld_fil = filter(wrld, pop < 297517)
```

```
## Simple feature collection with 3 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -73.297 ymin: -22.39998 xmax: 167.8449 ymax: 83.64513
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   iso_a2     name_long     continent region_un        subregion
## 1     GL     Greenland North America  Americas Northern America
## 2     NC New Caledonia       Oceania   Oceania        Melanesia
## 3     VU       Vanuatu       Oceania   Oceania        Melanesia
##                type   area_km2    pop  lifeExp gdpPercap
## 1           Country 2206644.44  56295       NA        NA
## 2        Dependency   23219.01 268000 77.57317        NA
## 3 Sovereign country    7490.04 258883 71.91832  2891.973
##                         geometry
## 1 MULTIPOLYGON (((-46.76379 8...
## 2 MULTIPOLYGON (((165.7799898...
## 3 MULTIPOLYGON (((167.8448767...
```

---
### Attributes - the **dplyr** package

- `mutate()`:

```r
wrld_mut = mutate(wrld, pop_density = pop/area_km2)
```

```
## Simple feature collection with 3 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 11.6401 ymin: -17.93064 xmax: 75.15803 ymax: 42.68825
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   iso_a2   name_long continent region_un       subregion              type
## 1     AF Afghanistan      Asia      Asia   Southern Asia Sovereign country
## 2     AO      Angola    Africa    Africa   Middle Africa Sovereign country
## 3     AL     Albania    Europe    Europe Southern Europe Sovereign country
##    area_km2      pop  lifeExp gdpPercap pop_density
## 1  652270.1 31627506 60.37446  1844.022    48.48836
## 2 1245463.7 24227524 52.26688  6955.960    19.45261
## 3   29694.8  2893654 77.83046 10698.525    97.44649
##                         geometry
## 1 MULTIPOLYGON (((61.21081709...
## 2 MULTIPOLYGON (((16.32652835...
## 3 MULTIPOLYGON (((20.59024743...
```

---
### Attributes - the **dplyr** package

- `summarize()`:

```r
wrld_sum1 = summarize(wrld, pop_sum = sum(pop, na.rm = TRUE), 
                      pop_mean = mean(pop, na.rm = TRUE), 
                      pop_median = median(pop, na.rm = TRUE))
```

```
## Simple feature collection with 1 feature and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##      pop_sum pop_mean pop_median                       geometry
## 1 7208968708 42910528   10521458 MULTIPOLYGON (((-159.208183...
```

---
### Attributes - the **dplyr** package

- `summarize()`:

```r
wrld_sum1 = wrld %>% 
        group_by(continent) %>% 
        summarize(pop_sum = sum(pop, na.rm = TRUE), 
                      pop_mean = mean(pop, na.rm = TRUE), 
                      pop_median = median(pop, na.rm = TRUE))
```

```
## Simple feature collection with 3 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 55.38525
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 3 x 5
##    continent    pop_sum pop_mean pop_median          geometry
##       <fctr>      <dbl>    <dbl>      <dbl>  <simple_feature>
## 1     Africa 1147005839 23895955   14129805 <MULTIPOLYGON...>
## 2 Antarctica          0      NaN         NA <MULTIPOLYGON...>
## 3       Asia 4306025131 95689447   18772481 <MULTIPOLYGON...>
```

---
### CRS and Reprojection

```r
st_crs(wrld)
```

```
## $epsg
## [1] 4326
## 
## $proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
## 
## attr(,"class")
## [1] "crs"
```

- The `st_transform()` can be used to transform coordinates

```r
wrld_3857 = st_transform(wrld, 3857)
st_crs(wrld_3857)
```

```
## $epsg
## [1] 3857
## 
## $proj4string
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
## 
## attr(,"class")
## [1] "crs"
```

---
### CRS and Reprojection

![](figs/coord_compare.png)

---
### Spatial operations

- Spatial subsetting
- Spatial joining/aggregation
- Topological relations
- Distances
- Spatial geometry modification

```r
continents = wrld %>% 
        st_set_precision(1000) %>% 
        st_union()
plot(continents)
```

---
### The **raster** package

The **raster** package consists of method and classes for raster processing. It allows to:

- Read and write raster data
- Perform raster algebra and raster manipulations
- Work on large datasets due to its ability to process data in chunks
- Visualize raster data
- Many more...

Basic **raster** operations are illustrated at https://rpubs.com/etiennebr/visualraster.

This package has three object classes:

- `RasterLayer` - for single-layer objects
- `RasterStack` - for multi-layer objects from separate files or a few layers from a single file
- `RasterBrick` - for multi-layer objects linked to a single file

```r
library(raster)
```
---
### Reading

```r
dem = raster("data/srtm.tif")
dem
```

```
## class       : RasterLayer 
## dimensions  : 720, 1200, 864000  (nrow, ncol, ncell)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : -85.00042, -84.00042, 38.79958, 39.39958  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jakub/Science/misc/gis_with_r_how_to_start/data/srtm.tif 
## names       : srtm 
## values      : 109, 323  (min, max)
```

---
### Writing

```r
writeRaster(dem, "data/new_dem.tif")
```

```r
writeRaster(dem, "data/new_dem2.tif", 
            datatype = "FLT4S", options=c("COMPRESS=DEFLATE"))
```

```r
writeFormats()
```

```
##       name        long_name                                
##  [1,] "raster"    "R-raster"                               
##  [2,] "SAGA"      "SAGA GIS"                               
##  [3,] "IDRISI"    "IDRISI"                                 
##  [4,] "IDRISIold" "IDRISI (img/doc)"                       
##  [5,] "BIL"       "Band by Line"                           
##  [6,] "BSQ"       "Band Sequential"                        
##  [7,] "BIP"       "Band by Pixel"                          
##  [8,] "ascii"     "Arc ASCII"                              
##  [9,] "CDF"       "NetCDF"                                 
## [10,] "big"       "big.matrix"                             
## [11,] "ADRG"      "ARC Digitized Raster Graphics"          
## [12,] "BMP"       "MS Windows Device Independent Bitmap"   
## [13,] "BT"        "VTP .bt (Binary Terrain) 1.3 Format"    
## [14,] "CTable2"   "CTable2 Datum Grid Shift"               
## [15,] "EHdr"      "ESRI .hdr Labelled"                     
## [16,] "ELAS"      "ELAS"                                   
## [17,] "ENVI"      "ENVI .hdr Labelled"                     
## [18,] "ERS"       "ERMapper .ers Labelled"                 
## [19,] "GPKG"      "GeoPackage"                             
## [20,] "GS7BG"     "Golden Software 7 Binary Grid (.grd)"   
## [21,] "GSBG"      "Golden Software Binary Grid (.grd)"     
## [22,] "GTiff"     "GeoTIFF"                                
## [23,] "GTX"       "NOAA Vertical Datum .GTX"               
## [24,] "HFA"       "Erdas Imagine Images (.img)"            
## [25,] "IDA"       "Image Data and Analysis"                
## [26,] "ILWIS"     "ILWIS Raster Map"                       
## [27,] "INGR"      "Intergraph Raster"                      
## [28,] "ISCE"      "ISCE raster"                            
## [29,] "ISIS2"     "USGS Astrogeology ISIS cube (Version 2)"
## [30,] "ISIS3"     "USGS Astrogeology ISIS cube (Version 3)"
## [31,] "KRO"       "KOLOR Raw"                              
## [32,] "LAN"       "Erdas .LAN/.GIS"                        
## [33,] "Leveller"  "Leveller heightfield"                   
## [34,] "MBTiles"   "MBTiles"                                
## [35,] "MRF"       "Meta Raster Format"                     
## [36,] "netCDF"    "Network Common Data Format"             
## [37,] "NITF"      "National Imagery Transmission Format"   
## [38,] "NTv2"      "NTv2 Datum Grid Shift"                  
## [39,] "NWT_GRD"   "Northwood Numeric Grid Format .grd/.tab"
## [40,] "PAux"      "PCI .aux Labelled"                      
## [41,] "PCIDSK"    "PCIDSK Database File"                   
## [42,] "PDF"       "Geospatial PDF"                         
## [43,] "PNM"       "Portable Pixmap Format (netpbm)"        
## [44,] "RMF"       "Raster Matrix Format"                   
## [45,] "ROI_PAC"   "ROI_PAC raster"                         
## [46,] "RST"       "Idrisi Raster A.1"                      
## [47,] "SAGA"      "SAGA GIS Binary Grid (.sdat)"           
## [48,] "SGI"       "SGI Image File Format 1.0"              
## [49,] "Terragen"  "Terragen heightfield"
```

---
### **raster** structure

```r
dem
```

```r
inMemory(dem)
```

```
## [1] FALSE
```

---
### Attributes

- the `getValues` function returns values from a raster object:

```r
values_dem = getValues(dem)
```

```
##  [1] 191 194 202 204 198 193 189 186 187 188 192 206 221 229 236 238 241
## [18] 236 226 217 209 205 201 198 203 211 218 220 225 232 237 236 238 242
## [35] 247 254 257 261 260 258 264 270 272 276 280 281 276 266 266 270
```

---
### Attributes

.pull-left[

```r
new_dem = dem + 50
```

<img src="gis_with_r_start_files/figure-html/unnamed-chunk-51-1.png" style="display: block; margin: auto;" />
]

.pull-right[

```r
new_dem2 = dem * new_dem
```

<img src="gis_with_r_start_files/figure-html/unnamed-chunk-53-1.png" style="display: block; margin: auto;" />
]

---
### CRS and Reprojection

```r
crs(dem)
```

```
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```

```r
dem3857 = projectRaster(dem, 
                         crs="+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0
                         +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0
                         +units=m +no_defs")
```

http://spatialreference.org/

---
### CRS and Reprojection

![](figs/coord_compare_raster.png)

---
### Extract

```r
ham_cities_sp = as(ham_cities_sf, "Spatial")
raster::extract(dem, ham_cities_sp)
```

```
##  [1] 276 160 157 161 188 226 260 202 184 261 171 243 258 260 165 178 173
## [18] 162 162
```

```r
ham_cities_sf$dem = raster::extract(dem, ham_cities_sp)
ham_cities_sf
```

```
## Simple feature collection with 19 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -84.81995 ymin: 39.07862 xmax: -84.26384 ymax: 39.29034
## epsg (SRID):    NA
## proj4string:    NA
##             name place                       geometry dem
## 1          Delhi  town POINT (-84.6052225 39.0950595) 276
## 2        Milford  town  POINT (-84.2958988 39.174625) 160
## 3      Covington  town  POINT (-84.508371 39.0836224) 157
## 4       Harrison  town POINT (-84.8199516 39.2619993) 161
## 5     Cincinnati  city POINT (-84.5124602 39.1014537) 188
## 6     Springdale  town  POINT (-84.4852213 39.287002) 226
## 7       Sycamore  town  POINT (-84.347427 39.2800883) 260
## 8        Norwood  town POINT (-84.4596641 39.1556149) 202
## 9       Loveland  town POINT (-84.2638406 39.2689562) 184
## 10   Fort Thomas  town POINT (-84.4483432 39.0786242) 261
## 11 Saint Bernard  town POINT (-84.4985541 39.1670033) 171
## 12    Montgomery  town POINT (-84.3539662 39.2282432) 243
## 13   Forest Park  town POINT (-84.5041108 39.2903353) 258
## 14      Blue Ash  town POINT (-84.3782817 39.2320073) 260
## 15        Ludlow  town POINT (-84.5474435 39.0925598) 165
## 16   Sharonville  town POINT (-84.4132779 39.2681145) 178
## 17       Reading  town POINT (-84.4421641 39.2236694) 173
## 18      Bellevue  town POINT (-84.4826482 39.1067078) 162
## 19       Newport  town POINT (-84.4919524 39.0889469) 162
```

---
### Crop

```r
ham84 = st_transform(ham, 4326)
ham_sp = as(ham84, "Spatial")
```

.pull-left[

```r
dem_crop = crop(dem, ham_sp)
```

<img src="gis_with_r_start_files/figure-html/unnamed-chunk-61-1.png" style="display: block; margin: auto;" />
]

.pull-right[

```r
dem_mask = mask(dem_crop, ham_sp)
```

<img src="gis_with_r_start_files/figure-html/unnamed-chunk-63-1.png" style="display: block; margin: auto;" />
]

---
### Reading spatial data - data packages

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Package name </th>
   <th style="text-align:left;"> Description </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> osmdata </td>
   <td style="text-align:left;"> Download and import of OpenStreetMap data </td>
  </tr>
  <tr>
   <td style="text-align:left;"> raster </td>
   <td style="text-align:left;"> The getData() function downloads and imports administrative country, SRTM/ASTER elevation, WorldClim data </td>
  </tr>
  <tr>
   <td style="text-align:left;"> rnaturalearth </td>
   <td style="text-align:left;"> Functions to download Natural Earth vector and raster data, including world country borders </td>
  </tr>
  <tr>
   <td style="text-align:left;"> tigirs </td>
   <td style="text-align:left;"> Census TIGER/Line shapefiles in R </td>
  </tr>
  <tr>
   <td style="text-align:left;"> USAboundaries </td>
   <td style="text-align:left;"> Contemporary state, county, and congressional district boundaries for the United States of America, as well as historical boundaries from 1629 to 2000 for states and counties </td>
  </tr>
  <tr>
   <td style="text-align:left;"> rnoaa </td>
   <td style="text-align:left;"> An interface to many NOAA data sources </td>
  </tr>
  <tr>
   <td style="text-align:left;"> rWBclimate </td>
   <td style="text-align:left;"> The World Bank climate data </td>
  </tr>
</tbody>
</table>

---
### Basic maps

- Basic maps of `sf` object can be quickly created using the `plot()` function:

```r
plot(wrld[0])
```

```r
plot(wrld["pop"])
```

![](figs/plot_compare.png)

---
### Basic maps

```r
plot(dem)
```

---
### **rasterVis**

- https://oscarperpinan.github.io/rastervis/
- http://www.colorbrewer.org

```r
library(rasterVis)
my_theme = rasterTheme(region=brewer.pal("RdYlGn", n = 9))
p = levelplot(dem_crop, margin = FALSE, par.settings = my_theme)
p = p + layer(sp.lines(ham_sp, lwd = 3, col = "darkgrey"))
p + layer(sp.points(ham_cities_sp, pch = 19, col = "black"))
```

---
### **tmap**

- https://cran.r-project.org/web/packages/tmap/vignettes/tmap-nutshell.html

```r
library(tmap)
tm_shape(wrld, projection="wintri") +
        tm_polygons("lifeExp", style="pretty", palette="RdYlGn",
                    auto.palette.mapping=FALSE, title=c("Life expectancy")) +
        tm_style_grey()
```

---
### **leaflet**

```r
library(leaflet)
leaflet(ham_sp) %>%
        addPolygons() %>% 
        addMarkers(data=ham_cities_sp, popup=~as.character(name))
```

<div id="htmlwidget-7e5f633d3c2f0c8bf9c9" style="width:100%;height:400px;" class="widgetframe html-widget"></div>
<script type="application/json" data-for="htmlwidget-7e5f633d3c2f0c8bf9c9">{"x":{"url":"gis_with_r_start_files/figure-html//widgets/widget_unnamed-chunk-71.html","options":{"xdomain":"*","allowfullscreen":false,"lazyload":false}},"evals":[],"jsHooks":[]}</script>

---
class: center, middle
# Part IV - Examples

---
### Static visualization (I)

---
### Static visualization (I)

```r
library(tidyverse)
library(sf)
library(spData)

world = st_read(system.file("shapes/world.gpkg", package = "spData"))
world_globe = sf::st_transform(world,
        "+proj=ortho +x_0=0 +y_0=0 +lon_0=-100 +lat_0=40")

p1 = ggplot(world_globe) +
        geom_sf(data = world_globe, aes(fill = lifeExp))  +
        scale_fill_viridis_c(name = "Life expectancy: ") +
        theme_void()
p1
ggsave(p1, file="figs/01_world_globe.png", 
       width = 70, height = 60, units = "mm", scale = 2)
```

---
### Static visualization (II)

---
### Static visualization (II)
https://hafen.github.io/geofacet/

```r
library(geofacet)
library(ggplot2)
data("state_unemp", package = "geofacet")
p2 = ggplot(state_unemp, aes(year, rate)) +
        geom_line() +
        facet_geo(~state, grid = "us_state_grid2", label = "name") +
        scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
        labs(title = "Seasonally Adjusted US Unemployment Rate 2000-2016",
             caption = "Data Source: bls.gov",
             x = "Year",
             y = "Unemployment Rate (%)") +
        theme(strip.text.x = element_text(size = 2)) + 
        theme_bw()
ggsave(p2, file="figs/02_state_unemp.png", 
       width = 70, height = 60, units = "mm", scale = 4)        
```

---
### Static visualization (III)

---
### Static visualization (III)
https://github.com/rCarto/linemap

```r
library(sf)
library(tidycensus)
library(linemap)
# A key can be acquired at http://api.census.gov/data/key_signup.html
# census_api_key("your key")
oh_pop = get_acs(geography = "tract", variables = "B01003_001E", 
                  state = "OH", geometry = TRUE) 
plot(oh_pop$geometry)
oh_pop
oh_pop2 = st_transform(oh_pop, 26917)

oh_pop2_grid = getgrid(x = oh_pop2, cellsize = 2500, var = "estimate")

png("figs/03_ohio_linemap.png", width = 650, height = 600)
opar = par(mar = c(0,0,0,0))
plot(st_geometry(oh_pop2), col="lightsteelblue3", border = NA, bg = "lightsteelblue1")
linemap(x = oh_pop2_grid, var = "estimate", k = 3, threshold = 1,
        col = "lightsteelblue3", border = "white", lwd = 0.8,
        add = TRUE)
par(opar)
dev.off()
```

---
### Animation

<img src="figs/04_anim.gif" width="550px" style="display: block; margin: auto;" />
---
### Animation

```r
library(gapminder)
library(tidyverse)
library(sf)
library(tmap)
library(spData)
europe_area = st_polygon(list(rbind(c(-13, 34), c(-13, 72),
                                    c(43, 72), c(43, 34), c(-13, 34)))) %>% 
        st_sfc(crs = st_crs(world)) 
europe = st_intersection(world, europe_area) 
europe_gdp = europe %>% 
        inner_join(gapminder, by = c("name_long" = "country"))
m1 = tm_shape(europe) + 
        tm_fill("grey") + 
     tm_shape(europe_gdp) + 
        tm_fill("gdpPercap.y", title = "GDP per capita: ") +
        tm_borders() + 
        tm_facets(by = "year", nrow = 1, ncol = 1, drop.units = TRUE) + 
        tm_legend(text.size = 1, title.size = 1.2, 
                  position = c("left", "TOP"), height = 0.3)

animation_tmap(m1, filename = "figs/04_anim.gif", width = 650, height = 600)
```

---
### Interactive visualization

---
### Interactive visualization

```r
library(tigris)
library(leaflet)
library(sf)
library(tidyverse)
library(osmdata)
library(widgetframe)
oh = counties(state="OH")
hamilton = oh %>% 
        st_as_sf() %>% 
        filter(NAME == "Hamilton") %>% 
        st_transform(4326)
pubs = opq(bbox = "Cincinnati OH") %>% 
  add_osm_feature(key = "amenity", value = "pub") %>% 
  osmdata_sf()
l1 = leaflet(hamilton) %>%
        addTiles(group = "OSM") %>%
        addProviderTiles(providers$Stamen.Watercolor, group = "Watercolor") %>% 
        addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>% 
        addLayersControl(    baseGroups = c("OSM", "Watercolor", "Toner Lite")) %>% 
        addPolygons() %>% 
        addMarkers(data = pubs$osm_points, popup=~as.character(name))
saveWidget(l1, file = "05_leaflet.html", selfcontained = FALSE)
```

---
### Remote sensing calculations

---
### Remote sensing calculations

```r
library(raster)
# download_landsat8 = function(destination_filename, band){
#         filename = paste0("http://landsat-pds.s3.amazonaws.com/L8/038/034/LC80380342015230LGN00/LC80380342015230LGN00_B", band, ".TIF")
#         download.file(filename, destination_filename, method='auto')
# }
# download_landsat8("data/landsat_b4.tif", band = 4)
# download_landsat8("data/landsat_b5.tif", band = 5)
b4 = raster("data/landsat_b4.tif")
b5 = raster("data/landsat_b5.tif")

ndvi_calc = function(b4, b5){
        ndvi = (b5 - b4)/(b5 + b4)
        ndvi
}

my_ndvi = overlay(b4, b5, fun = ndvi_calc)

png("figs/06_ndvi.png", width = 650, height = 500)
par(mfrow=c(1,2))
plot(my_ndvi, main = "NDVI")
hist(my_ndvi, main = "", xlab = "NDVI")
dev.off()
```

---
### Terrain calculations

---
### Terrain calculations

```r
library(raster)
# getData('SRTM', lon = 5, lat = 45, path = "data")
srtm = raster("data/srtm_38_04.tif")
my_area = extent(c(8, 10, 41.3, 43.1))
my_srtm = crop(srtm, my_area)
plot(my_srtm)

library(mapview)
mapView(my_srtm)

my_terrain = terrain(my_srtm, opt = c("slope", "aspect", "tpi", "tri", "roughness", "flowdir"))

png("figs/07_terrain.png", width = 650, height = 500)
plot(my_terrain)
dev.off()
```

---
### Geocoding

What's the coordinates of New York, Los Angeles, Chicago, Houston, and Philadelphia?

---
### Geocoding

```r
# http://happygitwithr.com/api-tokens.html
library(opencage)
library(tidyverse)
library(sf)
library(mapview)
library(htmlwidgets)

my_places = c("New York", "Los Angeles", "Chicago",
              "Houston", "Philadelphia")

output = my_places %>%
        map(opencage_forward, limit = 1) %>% 
        map_df(1)

output_sf = st_as_sf(output, 
                     coords = c("geometry.lng", "geometry.lat"), crs = 4326)
plot(output_sf$geometry, axes = TRUE)

m1 = mapview(output_sf)
mapshot(m1, "08_geocoding.html")
```

---
class: center, middle
# Part V - The end

---

![](figs/geocompr.png)

The online version of the book is at http://robinlovelace.net/geocompr/ and its source code at https://github.com/robinlovelace/geocompr.

We encourage contributions on any part of the book, including:

- Improvements to the text, e.g. clarifying unclear sentences, fixing typos (see guidance from [Yihui Xie](https://yihui.name/en/2013/06/fix-typo-in-documentation/))
- Changes to the code, e.g. to do things in a more efficient way
- Suggestions on content (see the projects [issue tracker](https://github.com/Robinlovelace/geocompr/issues) and the [work-in-progress](https://github.com/Robinlovelace/geocompr/tree/master/work-in-progress) folder for chapters in the pipeline)

Please see [our_style.md](https://github.com/Robinlovelace/geocompr/blob/master/our_style.md) for the books style.

---
### References

- Lovelace Robin, Nowosad Jakub, Meunchow Jannes. 2018. Geocomputation with R. CRC Press.
- Pebesma, Edzer. 2017. Sf: Simple Features for R. https://github.com/r-spatial/sf/.
- Pebesma, Edzer, and Roger Bivand. 2017. Sp: Classes and Methods for Spatial Data. https://CRAN.R-project.org/package=sp.
- Hijmans, Robert J.. 2016b. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.
- Hadley Wickham (2017). tidyverse: Easily Install and Load 'Tidyverse' Packages. R package version 1.1.1.
  https://CRAN.R-project.org/package=tidyverse
  
---
### References

- Ryan Hafen. 2017. geofacet: 'ggplot2' Faceting Utilities for Geographical Data. R package version 0.1.5. https://CRAN.R-project.org/package=geofacet
- Oscar Perpinan Lamigueiro and Robert Hijmans. 2016, rasterVis. R package version 0.41. 
- Tennekes, Martijn. 2017. Tmap: Thematic Maps. https://CRAN.R-project.org/package=tmap.
- Cheng, Joe, Bhaskar Karambelkar, and Yihui Xie. 2017. Leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library. https://CRAN.R-project.org/package=leaflet.
- Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2017. Mapview: Interactive Viewing of Spatial Objects in R. https://CRAN.R-project.org/package=mapview.
- Timothée Giraud. 2017. linemap: Line Maps. R package version 0.1.0.
- Bhaskar Karambelkar. 2017. widgetframe: 'Htmlwidgets' in Responsive 'iframes'. R package version 0.3.0.
  https://CRAN.R-project.org/package=widgetframe

---
### References

- Kyle Walker. 2017. tidycensus: Load US Census Boundary and Attribute Data as 'tidyverse' and 'sf'-Ready Data Frames. R package
  version 0.3. https://CRAN.R-project.org/package=tidycensus
- Roger Bivand, Jakub Nowosad and Robin Lovelace. 2017. spData: Datasets for Spatial Analysis. R package version 0.2.6.0.
  https://github.com/Nowosad/spData  
- Jennifer Bryan. 2015. gapminder: Data from Gapminder. R package version 0.2.0. https://CRAN.R-project.org/package=gapminder
- Kyle Walker. 2017. tigris: Load Census TIGER/Line Shapefiles into R. R package version 0.5.5. https://github.com/walkerke/tigris
- Mark Padgham, Bob Rudis, Robin Lovelace and Maëlle Salmon. 2017. osmdata: Import 'OpenStreetMap' Data as Simple Features or Spatial
  Objects. R package version 0.0.5. https://CRAN.R-project.org/package=osmdata
- Maëlle Salmon. 2017. opencage: Interface to the OpenCage API. R package version 0.1.2. https://CRAN.R-project.org/package=opencage