---
title: "The Potato and Urbanisation in Europe"
output: html_notebook
author: "PhD Unibo GIS class"
date: "April 2019"
---

## Some Technicalities in the Beginning

First, create an R project and put this file into the project folder. This has several advantages such as being able to reference to indirect file paths within the project directory (so if you share this with colleagues or copy it to another folder, you can just run the code as is and it works - no need to redefine all filepaths).

Loading packages and data.

```{r}
rm(list = ls()) # clear workspace
#install.packages(c("dplyr", "sf", "tmap", "ggplot2", "raster", "spatialEco"))
library(dplyr)
library(sf)
# define project CRS
# for plotting
library(tmap)
library(ggplot2)
library(raster)
library(spatialEco) # zonal statistics (masking raster::extract)
projcrs <- "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
```



Get in the Bairoch data (you know it already from the LMEC course when we computed effective distances on the Roman roadnetwork)

```{r}
Bairoch <- st_read("Bairoch_v2.gpkg")
tm_shape(Bairoch) + tm_dots()
head(Bairoch)
sum(Bairoch$X10, na.rm = T)
sum(Bairoch$X18)#, na.rm = T)
qplot(Bairoch$X10, binwidth = 5)
Bairoch %>% filter(X10 > 200) # the dplyr way
filter(Bairoch, X10 > 200) # the dplyr way
Bairoch[Bairoch$X10 > 200, ]
#test <- Bairoch
#st_geometry(test) <- NULL
#test[test$X10 > 200, ]

```


What if the data comes as a csv file? (just to demonstrate how that would work to make it an sf spatial dataframe immediately)

```{r}
Bairoch_v2_wide <- read_csv("../../../Bairoch_v2/Bairoch_v2_wide.csv")
Bairoch_v2_wide <- st_as_sf(Bairoch_v2_wide, coords = c("lon", "lat"), crs = projcrs)
```


Let's get rid of the out-of-Europe cities (creating a bounding box that defines the extent of our study area)

```{r}
bbox <- raster(ext = extent(-10.0, 19.0, 35.0, 57.0), res = 10, crs = projcrs)
bbox <- rasterToPolygons(bbox)
bbox <- st_as_sf(bbox)
Bairoch <- st_intersection(Bairoch, bbox)
```


Get in the potato suitability raster from FAO GAEZ

```{r}
potato <- raster("FAO-GAEZ_potato/data.asc")

potato <- crop(potato, raster(ext = extent(-10.0, 19.0, 35.0, 57.0), crs = projcrs))
plot(potato)
qtm(potato) + tm_shape(Bairoch) + tm_dots()

```

Before we dive into any more sophisticated analysis, let's just do the most basic thing and use one of the bread and butter GIS algorithms: the buffer
ATTENTION: CRS

```{r}
st_crs(Bairoch)
library(mapview)
mapview(st_buffer(Bairoch, .25))

buffersize <- .1 # in degrees? (1 something in the range of 100km, depending on the position of course)
buffersize_m <- 1e4 # in m for the raster::extract

# for the people that got an error, you need to project the raster
crs(potato)
#potato <- projectRaster(potato, crs = projcrs)
#potatobuffer <- zonal.stats(as(st_buffer(Bairoch, buffersize), "Spatial"), potato, mean, plot = FALSE)

# raster approach: drawing the buffer within extract() speeds up computation enormously
# buffer for raster::extract is in METERS!
potatobuffer_5 <- raster::extract(potato, Bairoch, fun = mean, buffer = 5000)
potatobuffer_10 <- raster::extract(potato, Bairoch, fun = mean, buffer = 10000)
potatobuffer_25 <- raster::extract(potato, Bairoch, fun = mean, buffer = 25000)
potatobuffer_30 <- raster::extract(potato, Bairoch, fun = mean, buffer = 30000)
potatobuffer_50 <- raster::extract(potato, Bairoch, fun = mean, buffer = 50000)
# now we attach this as a new column to the data frame (sf object)

#Bairoch$potato <- potatobuffer
Bairoch$potato5 <- potatobuffer_5
Bairoch$potato10 <- potatobuffer_10
Bairoch$potato25 <- potatobuffer_25
Bairoch$potato30 <- potatobuffer_30
Bairoch$potato50 <- potatobuffer_50
#qplot(Bairoch$potato10, Bairoch$potato25) + geom_abline(slope = 1, col = "red")
ggplot(data = Bairoch) + geom_point(aes(potato50, potato30)) + geom_abline(slope = 1, col = "red")

```

Of course we investigate the preliminary results a bit

```{r}
tm_shape(potato) + tm_raster() +
  tm_shape(Bairoch) + tm_dots(size = .1, col = "potato30", breaks = seq(0, 100, by = 10), palette = "OrRd")

# quick OLS
lm(X10 ~ potato, data = Bairoch)
summary(lm(X10 ~ potato1, data = Bairoch))
summary(lm(X18 ~ potato50, data = Bairoch))
```

Some more complicated descriptives, summing over all columns and check how much urban population (intensive margin) and how many cities (extensive margin) we have within our bounding box for each century

```{r}
# sums over the columns
as.data.frame(Bairoch) %>% dplyr::select(X10:X18.5) %>% summarise_all(funs(sum), na.rm = T)
as.data.frame(Bairoch) %>% dplyr::select(X10:X18.5) %>% summarise_all(funs(sum), na.rm = T)
# counts non-NA occurences over columns ()
as.data.frame(Bairoch) %>% dplyr::select(X10:X18) %>% sapply(function(x) sum(!is.na(x)))
as.data.frame(Bairoch) %>% dplyr::select(X10:X18) %>% summarise_all(funs(sum(!is.na(.))))
```

Let's compute urban growth across centuries (what should we do with NA cities that enter the dataset? discuss!)

```{r}
# the base R way
Bairoch$g1750 <- Bairoch$X17.5 / Bairoch$X17
Bairoch$g1800 <- Bairoch$X18 / Bairoch$X17.5
Bairoch$g1850 <- Bairoch$X18.5 / Bairoch$X18
# let's check the histograms
qplot(Bairoch$g1750)
qplot(Bairoch$g1800)
qplot(Bairoch$g1850)
# who are the crazy-growers?
Bairoch[order(Bairoch$g1850, decreasing = T), ]
# dplyr way
Bairoch %>% arrange(desc(g1850))
# which towns have the best potato suitability?
Bairoch %>% arrange(desc(potato))
# let's postpone the descriptives for when we added in country information
```

Now we can actually do a preliminary regression a la Nunn, Qian....

```{r}
summary(lm(g1750 ~ potato30, data = Bairoch))
summary(lm(g1800 ~ potato30, data = Bairoch))
summary(lm(g1850 ~ potato30, data = Bairoch))
summary(Bairoch$g1850)
summary(lm(g1850 ~ log(potato), data = Bairoch))
```

We still don't have any other time-invariant information than the potatoe suitability attached to each of our polygon features.
For further regressions we might want to exploit within-unit-variation only or cluster our standard errors. 
So let's begin with geographic information.


```{r}
shape_EUR <- st_read("ne_10m_admin_0_map_subunits/ne_10m_admin_0_map_subunits.shp")
tm_shape(shape_EUR) + tm_polygons()

# since we are there, let's just dive a bit into mapplotting
names(shape_EUR)
tm_shape(shape_EUR) + tm_polygons("POP_EST") # the tmap way
# we can also use different projections
library(tmaptools)
shape_proj <- set_projection(shape_EUR, "robin") #gall
tm_shape(shape_proj) + tm_polygons("POP_EST") # the tmap way

# let's kick out non-Europe
## the nice thing about an sf object is that we can just apply standard matrix/data.frame subsetting syntax (we saw that in Python last time)

# forget about the following few lines for the moment
#world <- shape_EUR %>% group_by(CONTINENT) %>% summarize(pop = sum(POP_EST, na.rm = T))
#ggplot() + geom_sf(data = world, aes(fill = pop)) # the grammar of graphics (ggplot) way
# just another way of projecting
#ggplot(data = shape_EUR) +
#    geom_sf() +
#    coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")


# subset it the dplyr way
shape_EUR <- shape_EUR %>% filter(CONTINENT == "Europe")
object.size(shape_EUR)
names(shape_EUR)
# this would be the base-R way
#shape_EUR[shape_EUR$CONTINENT == "Europe", ]
ggplot(data = shape_EUR) +
    geom_sf() +
    coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

# as I know that some people are afraid to go away from GUI-type-of-GIS-software because one cannot "touch" things....
# there is also several ways in R to zoom in and out and click around, e.g. this one (others would be the leaflet package, mapview, or tmap_leaflet)
tm_shape(shape_EUR) + tm_polygons("POP_EST")
tmap_mode("view")
last_map()
tmap_mode("plot")

```

Our spatial sf object got a lot smaller already, but it still completely cluttered with redundant columns that come with the Natural Earth data.
Let's change that.

```{r}
# dplyr way
shape_EUR %>% dplyr::select(ADM0_A3, ADMIN, SU_A3, SUBUNIT, SUBREGION)
# base R way (can also subset by column names, not just by positon or T/F)
shape_EUR[, 1:3]
shape_EUR[, c("ADM0_A3", "ADMIN", "SU_A3", "SUBUNIT", "SUBREGION")]
# a more sophisticated base R way is to use %in% (if you don't know orderings or not all column names actually exist etc)
shape_EUR <- shape_EUR[, names(shape_EUR) %in% c("ADM0_A3", "ADMIN", "SU_A3", "SUBUNIT", "SUBREGION")]

```

Very good. Size reduced a lot. Left with 6 columns: 5 with content and the last, as usual, containts the geometry information.
As discussed earlier, we can make any sf object a completely standard data.frame by just dropping the geometries with st_geometry() <- NULL

Now we pass over the geographic information from the NaturalEarth file to our Bairoch cities.

```{r}
#intersected <- st_intersection(Bairoch, shape_EUR)
# check it (view!)
# does the job!
# let's make it default
Bairoch <- st_intersection(Bairoch, shape_EUR)
rm(intersected)
head(Bairoch)
```


Doing some visual exploration beforehand...

```{r}
ggplot(data = Bairoch, aes(potato30, g1850)) + geom_point() + # aes(col = SUBREGION)
  stat_smooth(method = "loess", formula = y ~ x, col = "blue") +
  stat_smooth(method = "lm", formula = y ~ x, col = "red") +
  facet_wrap(vars(SUBREGION))
```

Towards thinking about the map from the 1980ies that shows potato useage across Europe.
So let's also check where the action in terms of our y (urban growth) happens in space...

```{r}
tm_shape(Bairoch[!is.na(Bairoch$g1850), ]) + tm_dots(size = .1, col = "g1850")
```



Now we finally are ready to ramp up our R-metrics game a bit:

```{r}
library(estimatr)
lm_robust(g1850 ~ potato50, data = Bairoch)
summary(lm_robust(g1850 ~ potato50, data = Bairoch))

# let's check what the lm_robust can do for us... ?lm_robust

# from 10km to 50km turns from neg to pos
summary(lm_robust(g1850 ~ potato50, data = Bairoch, fixed_effects = ADMIN))


```



```{r BigtownSubset}
# kick out all the cities below x
Bairoch40k <- Bairoch %>% mutate_at(vars(X10:X18.5), funs(replace(., which(. < 20), NA)))

Bairoch40k$g1750 <- Bairoch40k$X17.5 / Bairoch40k$X17
Bairoch40k$g1800 <- Bairoch40k$X18 / Bairoch40k$X17.5
Bairoch40k$g1850 <- Bairoch40k$X18.5 / Bairoch40k$X18
# let's check the histograms
qplot(Bairoch40k$g1750)
qplot(Bairoch40k$g1800)
qplot(Bairoch40k$g1850)

summary(lm(g1750 ~ potato1, data = Bairoch40k))
summary(lm(g1800 ~ potato1, data = Bairoch40k))
summary(lm(g1850 ~ potato1, data = Bairoch40k))

```


# Grid Level


Create the grid within a "bounding box"

```{r}

Grid <- raster(ext = extent(-10.0, 19.0, 35.0, 57.0), res = 1, crs = projcrs)
Grid <- rasterToPolygons(Grid) # make the raster a polygon so we can feed each cell with information
Grid <- st_as_sf(Grid) # and convert immediately to the sf class (sp will soon be outdated)
Grid$id_bbox <- 1:nrow(Grid)
```

Let's see what we got

```{r}
tm_shape(Grid) + tm_polygons() + 
  tm_shape(Bairoch) + tm_dots()
```


Now we cut out a "European shape" out of our rectangular Grid (think of cookiecutting)

```{r}
intersected <- st_intersection(shape_EUR, Grid)
# let's see what we got
tm_shape(intersected) + tm_polygons() + tm_shape(Bairoch) + tm_dots()

```


We have to get rid of all the countryborders within the gridcells. The split them apart.
But at the same time we want to retain the information of the EUR shapefile within our gridcells (country, region,...)

```{r}
union <- intersected %>% group_by(id_bbox) %>% summarise()
tm_shape(union) + tm_polygons() + tm_shape(Bairoch) + tm_dots()
# now let's give the new grid without water and countryborderintersections it's id
union$id <- 1:nrow(union)
```


By aggregating all polygons with the same id_bbox we lost the other information, but it is still in the `intersected` object.
Let's use it...

```{r}
intersected <- intersected %>% dplyr::select(-id_bbox) # kick out the id_bbox column because we have that also in the union layer and that duplicate would cause troubles when it comes to the naming
# the dplyr:: in front of the select just tells R to use the select function from the dplyr package (because there is another package that also has a function with the same name)
union <- st_join(union, intersected)
sum(duplicated(union$id_bbox)) # this created  alot of duplicates stacked on top of each other
# but it is a viable way to retain information
# as a rule we want that the country/region gets assigned to a gridcell when it has the largest area within
union$area <- st_area(union)
# now we kick out the duplicates
Grid_Bairoch <- union %>% group_by(id_bbox) %>% slice(which.max(area))
sum(duplicated(Grid_Bairoch$id_bbox)) # this created  alot of duplicates stacked on top of each other
# let's convert the area
Grid_Bairoch$area <- st_area(Grid_Bairoch)/1000000
# get an id
Grid_Bairoch$id <- 1:nrow(Grid_Bairoch)
tm_shape(Grid_Bairoch) + tm_polygons("area") + tm_shape(Bairoch) + tm_dots()
# and let's kick out our temporary files rightaway
# rm(intersected, union)
```

Now we can start to feed the Bairoch city info in our grid

```{r}
intersected <- st_intersection(Bairoch, Grid_Bairoch)
sum(duplicated(Grid$idbb))
# now we know the gridcell id for each Bairoch city
# let's use the trick from before again and group by the id and summ over all values
# this gives us the total urban pop in gridcell i
summed <- intersected %>% group_by(id) %>% summarise_at(vars(X10:X18.5), funs(sum(., na.rm = T)))
head(summed)
```

This newly obtain frame we now just merge onto the existing grid (a very standard "non-spatial" task)

```{r}
# first we drop the geometry from summed
st_geometry(summed) <- NULL
# now we merge onto the Grid using the unique identifier (id) from the Grid (that we defined before for each Bairoch city)
Grid_Bairoch <- merge(Grid_Bairoch, summed, by = "id", all.x = T)
# from before we have a lot of 0's and no NA's....
# ... we have to change that because we need non-existing stuff to be dropped in the regression, and not treated as 0
Grid_Bairoch <- Grid_Bairoch %>% mutate_at(vars(X10:X18.5), funs(replace(., which(. == 0), NA)))
# let's see it
tm_shape(Grid_Bairoch) + tm_polygons("X18")
tm_shape(Grid_Bairoch) + tm_polygons("X18", breaks = c(10, 100, 200, 300, 400, 600, 1200)) + tm_shape(Bairoch) + tm_dots()
```


What if we wanted to know the biggest city (or whatever other information, like river, main crop, biggest factory/firm/employer,...) and want to pass this information (as a column of strings) onto the grid?

```{r}
# this is a crude trick that I cooked up with some R streetsmarts, there is for sure a more elegant way doing this with 1 line
citynames <- intersected[, names(intersected) %in% c("id", "citybb", "X18.5")]
# still has a geometry, let's kick it
st_geometry(citynames) <- NULL
# now our rule is that we want the name of the biggest city in 1850 in each cell
citynames <- as.data.frame(citynames) %>% arrange(desc(X18.5)) %>% group_by(id) 
citynames <- citynames[!duplicated(citynames$id), ]
citynames <- citynames[, !names(citynames) %in% "X18.5"]
Grid_Bairoch <- merge(Grid_Bairoch, citynames, by = "id", all.x = T)

# let's see if we did a good job (interactively)
tplot <- tm_shape(Grid_Bairoch) + tm_polygons("X18", popup.vars = c("citybb", "id", "X18"))
tmap_leaflet(tplot)

```

We could also remove very tiny cells if we wanted to...

```{r}

```

Add the geoposition explicitly to the layer. Probably we need it at some point...? (latitude might make sense as a control)
Note that the FAO GAEZ potato suitability captures that already a bit in a way.

```{r}
Grid_Bairoch$lon <- as.numeric(st_coordinates(st_centroid(Grid_Bairoch))[, 1]) # lon
Grid_Bairoch$lat <- as.numeric(st_coordinates(st_centroid(Grid_Bairoch))[, 2]) # lat
```


Rasterstats on the Grid

```{r}

potatomean <- raster::extract(potato, Grid_Bairoch, fun = mean)
Grid_Bairoch$potato_mean <- as.vector(potatomean)
qplot(Grid_Bairoch$potato_mean)
tm_shape(Grid_Bairoch) + tm_polygons("potato_mean")
# this produces NAs, and it's not obvious why

potatomean <- zonal.stats(as(Grid_Bairoch, "Spatial"), potato, mean, plot = FALSE)
Grid_Bairoch$potato_mean1 <- potatomean
tm_shape(Grid_Bairoch) + tm_polygons("potato_mean1")
# let's check
qplot(Grid_Bairoch$potato_mean, Grid_Bairoch$potato_mean1)
# everything alright, straight line - perfect match

```

## Regressions on the Grid

Now we start with grid regressions. Same story from above.
First we need to calculate our y's

```{r}

Grid_Bairoch$g1750 <- Grid_Bairoch$X17.5 / Grid_Bairoch$X17
Grid_Bairoch$g1800 <- Grid_Bairoch$X18 / Grid_Bairoch$X17.5
Grid_Bairoch$g1850 <- Grid_Bairoch$X18.5 / Grid_Bairoch$X18
# let's check the histograms
qplot(Grid_Bairoch$g1750)
qplot(Grid_Bairoch$g1800)
qplot(Grid_Bairoch$g1850)

summary(lm(g1750 ~ potato_mean1, data = Grid_Bairoch))
summary(lm(g1800 ~ potato_mean1, data = Grid_Bairoch))
summary(lm(g1850 ~ potato_mean1, data = Grid_Bairoch))

```

