## Packages:
  - **data.table**	data wrangling based on high performance [website]() 
  - **dplyr** 		data wrangling, using the *tidy data* approach [website]()
  - **rgdal**		bindings for GDAL and PROJ4 libraries
  - **rgeos**		
  - **sp**  		1st gen (2003) vector spatial data I/O and manipulation [website](https://edzer.github.io/sp/)
  - **sf**	 		2nd gen (2015) vector spatial data I/O and manipulation [website](https://r-spatial.github.io/sf/)
  - **raster**	 	raster spatial data I/O and manipulation [website]()
  - **mapedit**		interactive editing for spatial objects [blogpost-1](http://r-spatial.org/r/2017/01/30/mapedit_intro.html) [blogpost-2](http://r-spatial.org/r/2017/06/09/mapedit_0-2-0.html)
  - **ggplot2** 	generic data visualization, including some mapping functionalities, using the [*Grammar of Graphics*](http://vita.had.co.nz/papers/layered-grammar.pdf) approach [website](ggplot2.org/)
  - **tmap** 		static mapping, using the same *Grammar of Graphics* layers' approach as ggplot2 [GH](https://github.com/mtennekes/tmap)
  - **leaflet**		interactive mapping [website](http://rstudio.github.io/leaflet/)
  - **mapview**		interactive mapping [website](https://r-spatial.github.io/mapview/)

## Data I/O 
```
# load packages
pkgs <- c('data.table', 'dplyr', 'rgdal', 'rgeos', 'sp', 'sf', 'raster', 'mapedit', 'ggplot2', 'tmap', 'leaflet', 'mapview')
lapply(pkgs, require, char = TRUE)

# Read in a vector spatial points (in the form of a dataframe with two columns `x_lon` and `y_lat` mapping to coordinates) using sp, adding the [WGS84 CRS](https://epsg.io/4326)
sp1 <- read.csv('poi.csv')
coordinates(sp1) <- ~x_lon+y_lat
proj4string(dts) <- CRS("+init=epsg:4326")

# Read in a vector spatial points using sf
sf1 <- st_as_sf('poi.csv', coords = c('x_lon', 'y_lat'), crs = 4326)

# Read in a vector polygons shapefile using rgdal for sp object
sp2 <- readOGR('path/to/file', 'areas')

# Read in a vector polygons shapefile using sf
sp2 <- st_read('areas.shp')

# convert an sf object into an sp object (as sp has been so important for more than a decade, it'll remain a legacy for a while. For example, raster+vector analysis can't be done over *sf* objects)
sp2 <- as(sf2, Class = 'Spatial')

# convert an sp object into an sf object
sf2 <- st_as_sf(sp2)

# write an sp object
writeOGR(sp1, )

# write an sf Points object (notice that without the `layer_options` argument the output won't include the coordinates, but only the attributes)
st_write(sf1, layer_options = "GEOMETRY=AS_XY")

# write an sf Polygon object
st_write(sf2)

# Read in a raster single layer tiff file using raster
s3 <- raster('grid.tif')

# plot an object using R core
plot(sf1)

# plot multiple objects together using R core (remember to run the code for all layers at once rather than line-by-line)
plot(sf1)
plot(sf2, add = TRUE)
plot(s3, add = TRUE)
...

# Read in a raster multi band tiff file using raster
s4 <- brick('col_field.tif')

# plot multiple layers in a raster brick in the same chart using raster
plotRGB(s4)

# Join a spatial object with a dataframe (= add attributes to a spatial object)


# for more maps use tmap
tm_shape(s1) + 
    tm_polygons(alpha = 0.5) +
    tm_shape(s3) + 
    tm_rgb()
```

## The `sf` package
Vector data are ...
Their structure can take one of three form: point, line, polygon, or any collection of them. In the case of polygons there can even be a collection of collections.

  - `st_read('/path/to/file/fname.ext')` 
    The function can recognize lots of different file types, inferring the right one from the file extension.
    The output is a sort of dataframe with a metadata header, containing information about the spatial structure of the features, and an additional *geometry* column containing the actual spatial references. This means, in particular, that sf objects can be manipulated natively in the tidyverse. Notice though that when selecting columns, the *geometry* column is carried over implicitly.
We saw that *geometry* is a list column. A list-column behaves, to a certain extent, like any other R column. The main difference is that instead of an atomic value, each observation in that column is a piece of an R list and this list can be as complex as needed. The list column allows you to store far more information in a single variable, and `sf` takes advantage of this by storing all geographic information for each feature in the list.

  - `st_write(s, '/path/to/file/fname.ext')`

  - `st_geometry(s)` extract only the spatial feature from the `sf` object. This can be useful in order to plot (only) the outlines of the spatial object, as otherwise the plot will draw a choropleth for every attribute in the table.

  - `st_sf(s)` 

  - `st_crs(s)` get/set the coordinate reference system (CRS). When the geographic data you read in with sf already has a CRS defined, sf will recognize and retain it. When the CRS is not defined you will need to define it yourself using either an **EPSG** code, or a **proj4** string, a full set of parameters spelled out in a string. See: [IOGP's EPSG Geodetic Parameter Dataset](http://www.epsg.org/) and [Registry](http://www.epsg-registry.org/), [SpatialReference](http://spatialreference.org/), [EPSG.io: Coordinate Systems Worldwide](http://epsg.io).
     It is also worth noting that some spatial operations require the involved spatial object(s) to have a **projected** CRS, for the operation to work or the result to be meaningful. An **unprojected** CRS uses *latitude* and *longitude* coordinates, and
references the earth as a three-dimensional object. A **projected** CRS uses X and Y coordinates as a two-dimensional
representation of the earth. To determine if an object has a projected CRS, simply look at the value set for **+proj=** at the start of the projection string; if it is **longlat**, then the corresponding CRS is unprojected.

  - `st_transform(s, crs = )` transform (or *project*) the spatial object to the CRS specified with the `crs` argument, that can be specified as:
    - a EPSG code, 
    - a proj4 string, 
    - a CRS of another spatial object. In order to perform any spatial analysis with multiple layers, all involved layers should share the same CRS (for example, if the layers do not share a common CRS, they may not correctly align on a plot).
For more information on the required conversions and transformations formulas, see [this document](http://www.epsg.org/Portals/0/373-07-2.pdf).

  - Indexing, subsetting, filtering, and augmenting an sf object, through recoding, transformations or merging, could be easily done with `dplyr` as if the spatial object were instead a standard tibble. One important difference is that the resulting data frames will always include the geometry variable by default, unless it's explicitly dropped using one of the following commands:
```
st_set_geometry(s, NULL) 
set_geometry(s) <- NULL
```
that removes the geometry from the original spatial object, coercing it to a data.frame.

  - `st_simplify(s, preserveTopology = TRUE, dTolerance = t)` Often, for the purpose of statistical mapping, the much detailed geographical maps are not only ..., but also ineffective. The process of reducing the details of map This is a good way to reduce object size for computation or quick display on the web. Most times the user will barely see the difference, unless  
    - `preserveTopology = TRUE` allows the borders to stay aligned, 
    - `dTolerance = ` set the amount of simplification allowed in meters: the higher the value, the more simplified will be the result.

  - `st_cast(s, )`

  - To count the number of vertices for every polygon in a multipolygon spatial object:
	```
	ss <- st_cast(s$geometry, "MULTIPOINT")
	ss <- sapply(ss, length)
	sum(ss)
	```

For the following spatial analysis functionalities, any layer involved should use a projected CRS, and if multiple layers are involved in the same analysis, all of them should have the same CRS.

### Measures

  - `st_area(s)` returns the area of each feature in *s*. Notice that the the result is a *units* object, and requires some additional processing, like `unclass` it, in order to use it in subsequent calculations

  - `st_length(s)` returns the perimeter of each features in *s*

  - `st_distance(s1, s2, by_element = , )` calculates the distances between every combination of the features of both `s1` and `s2`. Provides a useful feature-to-feature distance matrix as output and can be used for most distance calculation needs. 

### Unary Operations

  - `st_buffer(s, dist = )` Calculates a **buffer layers** of `dist` radius, around every point in `s`. The units of dist is the projection unit.

  - `st_centroid(s, of_largest_polygon = FALSE)` Returns the geometric **centroid** for every polygon in `s`. When `of_largest_polygon ` is set to TRUE, return centroid of the largest (sub)polygon of a MULTIPOLYGON rather than of the whole MULTIPOLYGON

  - `st_bbox` Returns the coordinates for the **bounding box** of any vector spatial object, which is the smallest rectangle surrounding *every* feature included in the object. These are helpful to create polygons to clip layers to a common area or identify regions of influence.

  - `st_make_grid(s, n = )` Subdivides any spatial object `s` using a grid of `n * n` equal squares. When `n = 1`, returns the polygon corresponding to the *bounding box* of `s`. If `n = c(n1, n2)` the grid is made up of `n1 * n2` equal rectangles.

  - `st_combine(s)`

  - `st_convex_hull(s)` returns the *convex hull* of a `MULTI*` feature, which is the smallest convex set that contains all the features in the collection. Notice that `st_convex_hull` computes a tight box around each one of your features individually so if you want to create a convex hull around a group of features you'll need to use first `st_union()` to combine individual features into a single multi-feature.

  - `st_boundary(x)`

  - `st_simplify(x, preserveTopology = FALSE, dTolerance = 0)`

  - `st_triangulate(x, dTolerance = 0, bOnlyEdges = FALSE)`

  - `st_voronoi(x, envelope, dTolerance = 0, bOnlyEdges = FALSE)`

  - `st_polygonize(x)`

  - `st_line_merge(x)`

  - `st_point_on_surface(x)`

  - `st_node(x)`

  - `st_segmentize(x, dfMaxLength, ...)`



### Binary Operations
The output of most of these and related functions that test relationships between two sets of features, is a special kind of list (of class `sgbp`) whose first element `[[1]]` shows the features from the second object that *acts*, according to the function used, on the first object. Likewise, `[[2]]` would show the features from the first object that acts on the second object.

  - `st_union(s1, s2, ..., by_feature = FALSE)` dissolve multiple features into one. If the elements of the collection are distinct `POINT`, or non-overlapping `POLYGON`, the resulting feature is always a `MULTI*`, that can be used for subsequent tasks like computing the *convex hull*.

  - `st_intersects(s1, s2)` returns the indices of the features in `s2` that intersect with `s1`

  - `st_intersection(s1, s2)` clips `s2` based on the spatial shape of `s1`

  - `st_contains(s1, s2)` returns the indices of the features in `s2` that are completely inside any feature in `s1`

  - `st_disjoint(s1, s2 = s1)` 

  - `st_touches(s1, s2)` 

  - `st_crosses(s1, s2)` 

  - `st_within(s1, s2)` 

  - `st_contains_properly(s1, s2)` 

  - `st_overlaps(s1, s2)` 

  - `st_equals(s1, s2)` 

  - `st_covers(s1, s2)` 

  - `st_covered_by(s1, s2)` 

  - `st_equals_exact(s1, s2, par)` 

  - `st_is_within_distance(s1, s2, dist)` 

  - `st_difference(s1, s2)`

  - `st_sym_difference(s1, s2)`

  - `st_snap(s1, s2, tolerance)`

  - `st_join(s1, s2)` for **spatial joins**, where all involved objects are of spatial type, and are linked using their spatial features. This is used when there is no common id that can link the features in the usual way.

  - `st_` **spatial relationship**

### Mapping

#### ggplot2

  - `geom_sf` operates over an sf object like any other layer in ggplot2, mapping variables from the data into aesthetics on the map through the `aes` function. Leaving off the aesthetic mapping will map the geometry alone.

  - `stat_sf`


  - `coord_sf`


#### tmap
tmap operates similarly to ggplot2, in the sense that it starts with a function to set up the map, in this case `tm_shape(data)`, and then is followed by layers and other instructions separated with the plus sign. There are a few dedicated *geometries* according to the type of the spatial object:
  - `tm_dots` or `tm_bubbles` for POINTS
  - `tm_lines` for LINES
  - `tm_polygons` for POLYGONS
  - `tm_rgb` for RASTERS
To build a choropleth map just add the argument `col` with a char vector identifying the variable(s) to be mapped.




## The `raster` package
When working with raster data, one of the most important things to keep in mind is that the raw data can be of two format:
  - *single-band* rasters with a single layer of values, where each cell value represents some attribute at that location
  - *multi-band* rasters with more than one layer, where each cell group of values represents different aspects of the same attribute, or an attribute at different points in time (the typical example is a color aerial photo in which there would be one band each representing red, green or blue light).

  - `raster('/path/to/file/fname.ext')` to read single-band rasters

  - `brick('/path/to/file/fname.ext')` to read a multi-band raster

Instead of storing raster objects in data frames, the raster package stores spatial data and metadata in specially designed R classes that contain slots where the data and metadata are stored. The data and metadata can be accessed using a suite of dedicated functions

  - `nlayers(s)` the number of bands/layers 

  - `extent(s)` the bounding box

  - `crs(s)` the coordinate reference system 

  - `ncell(s)` the number of grid cells

  - `res(rst)` returns the resolution of the grid

  - `inMemory(s)` tests if the raster values are loaded in memory. Raster data can be very big, depending on the extent and resolution (grid size). In order to deal with this, both the raster() and brick() functions are designed to only read in the actual raster values as needed, partly or the full set, to save space in memory. The raster values will be read in by default when performing spatial analysis operations that require it.

  - `getValues(s)` to manually read in memory the values for a raster already defined.

  - `plotRGB(s)` plot all the raster layers in a brick as a single image. Using the core `plot` function with a multi-band raster to draw as many charts as the number of layers in the raster.

  - `crs(s)` get/set the coordinate reference system (CRS)

  - `projectRaster(s, crs = '', method = c('ngb'))` transform the CRS associated with the spatial object to another one indicated by the argument `crs`. The new CRS can be 
    - a proj4 string, 
    - a EPSG code, in the form "+init=epsg:xxxxx", 
    - a CRS of another spatial object. In this case, the parameter `asText = TRUE` is needed to force the `crs()` function to output the CRS as a string.

  - 'aggregate(rst, fact = 2, fun = mean)` reduces the complexity (resolution) of the grid by a factor equal to `fact`, using `fun` to calculate the value of each output cells. The higher the value of `fact`, the more pixelated the result.

  - `reclassify(rst, rcl = M)` allows to recode the grid cells using the rules found along the rows of the matrix `M`, which could be structured with 2 or 3 columns:
    - In the former case, the first element of each row represent the old value, while the second value represent the new one.
    - When in the 3-cols form, the first and second elements indicate the min and respectively the max of the range of values to be set at the value represented by the third element. 

### Raster + vector analysis

  - `mask(rst, mask = sp)` limit the scope of a raster using the shape of a vector layer as mask. Any raster cells outside of the vector object are assigned NA values. The result is a raster object whose extent is not changed. Remember that both layers must have the same CRS.

  - `crop(rst, sp)` the output raster is the result of cutting the extent of the input raster using the bounding box of the vector layer, so that the two object share the same extent. Notice that no raster cells are set to NA.

  - `extract(rst, sp, fun = )` extract the values from the raster object corresponding to the location in the vector object, applying `fun` if necessary for neighborhood aggregation 

  - `overlay(rst1, rst2, fun = )` 


## Use Case: Trees streets density and average canopy comparison in New York Neighborhoods

Data available are:
  - `neighborhoods` boundaries as vector Polygons
  - `trees` for the city as vector Points
  - `canopy` for the city as raster

```
# Add areas to the neighborhoods object (need to unclass!)
neighborhoods <- mutate(neighborhoods, area = unclass(st_area(neighborhoods)))
# Compute the counts of all trees by hood, possibly renaming the count field and removing the geometry
trees <- count(trees, hood) %>% st_set_geometry(NULL) %>% rename(n_trees = n)
# Join neighborhoods and counts
neighborhoods <- left_join(neighborhoods, trees, by = "hood")
# Replace NA values with 0
neighborhoods <- mutate(neighborhoods, trees = ifelse(is.na(trees), 0, trees))
# Compute the tree density
neighborhoods <- mutate(neighborhoods, tree_density = n_tree / area)
# Transform the neighborhoods CRS to match the canopy layer (This is actually not strictly necessary because extract() can transform CRS on the fly. But it will be needed for plotting and other operations later)
neighborhoods <- st_transform(neighborhoods, crs = crs(canopy, asText = TRUE))
# Convert neighborhoods object to a Spatial sp object
neighborhoods <- as(neighborhoods, 'Spatial')
# Compute the mean of canopy values by neighborhood, and add it to the neighborhoods object
neighborhoods <- mutate(neighborhoods, avg_canopy = extract(canopy, neighborhoods, fun = mean) )
# Compute the correlation between tree density and average canopy
cor(neighborhoods$tree_density, neighborhoods$avg_canopy)
# Create a choropleth (color-coded) map of the tree density using ggplot
ggplot(neighborhoods) + 
  geom_sf(aes(fill = tree_density)) + 
  scale_fill_gradient(low = "#edf8e9", high = "#005a32")
# Create a choropleth map of the average canopy using tmap
tm_shape(neighborhoods) + 
	tm_polygons('avg_canopy', 
		palette = 'Greens', 
        style = 'quantile', n = 7, 
        title = 'Average tree canopy (%)'
	)
# Create a choropleth map of both the tree densities and average canopy using tmap
tm_shape(neighborhoods) + 
	tm_polygons(c('tree_density', 'avg_canopy'), 
		palette = 'Greens', 
        style = 'quantile', n = 7, 
        title = c('Tree Density', 'Average tree canopy (%)')
	)
# Create a raster aerial photo of Manhattan with the neighborhoods overlaying for comparison
tm_shape(manhattan) + 
    tm_rgb() + 
    tm_shape(neighborhoods) + 
    tm_polygons(col = "black", lwd = 0.5, alpha = 0.5)
# Combine the last two maps into a single output
tmap_arrange(tm1, tm2, asp = NA)
```