1 Introduction

Link to introductive video

Start libraries we need:

library(tidyverse)
library(sf)
library(tmap)
library(tmaptools)
library(lubridate)
library(ggthemes)
library(ggspatial)
library(raster)
library(rnaturalearth)
library(rgeos)

2 Physical map of Georgia

To create the map we need geolayers:

  • elevation
  • rivers
  • roads
  • populated places

You can find this kind of datasets from many different places. One option for example is OpenStreetMap data: http://download.geofabrik.de
Full zipped OSM dataset for Georgia is 44.9 MB. To speed up the process you can download selected data from here:

download.file("http://aasa.ut.ee/tbilisi/osm_georgia_selected.zip", destfile="osm_georgia_selected.zip")
unzip("osm_georgia_selected.zip")

The list of selected dataset contains 4 layers:

  • buildings
  • places
  • roads
  • waterways

Layers are quite big, which can make the processing a bit slow. Administrative borders and elevation are missing from downloaded archive.

2.0.1 Roads and rivers from OSM dataset

Lets look at roads and waterways:

# Import
roads <-  st_read("gis_osm_roads_free_1.shp")
## Reading layer `gis_osm_roads_free_1' from data source `C:\ANTO\loengud\tbilisi\2020\web\gis_osm_roads_free_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 156942 features and 10 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 39.5847 ymin: 40.94172 xmax: 46.724 ymax: 43.97154
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
#check the data:
glimpse(roads)
## Observations: 156,942
## Variables: 11
## $ osm_id   <fct> 5119612, 15950870, 17637354, 17637699, 17639488, 17673797,...
## $ code     <int> 5114, 5113, 5112, 5112, 5112, 5132, 5132, 5112, 5113, 5133...
## $ fclass   <fct> secondary, primary, trunk, trunk, trunk, trunk_link, trunk...
## $ name     <fct> NA, Ardahan-Posof yolu, NA, NA, NA, <U+10D7><U+10D1><U+10D...
## $ ref      <fct> NA, D 955, <U+10E1> 1, <U+10E1> 1, <U+10E1> 1, <U+10E1> 4,...
## $ oneway   <fct> B, F, F, B, B, F, F, B, F, F, B, B, B, T, B, B, B, B, F, B...
## $ maxspeed <int> 0, 20, 90, 0, 0, 0, 90, 90, 50, 0, 40, 50, 0, 0, 0, 0, 50,...
## $ layer    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ bridge   <fct> F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F...
## $ tunnel   <fct> F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F...
## $ geometry <LINESTRING [arc_degree]> LINESTRING (40.31842 43.230..., LINEST...

Roads layer contains 156942 road segments. Most probably we don’t need all of them. Let’s selected only major roads.

# Road types are stored in column 'fclass'. Get the list of unique records:
unique(roads$fclass)
##  [1] secondary      primary        trunk          trunk_link     primary_link  
##  [6] residential    tertiary       motorway_link  motorway       secondary_link
## [11] service        unclassified   track          unknown        living_street 
## [16] footway        track_grade3   path           track_grade2   track_grade4  
## [21] steps          pedestrian     track_grade1   tertiary_link  track_grade5  
## [26] bridleway      cycleway      
## 27 Levels: bridleway cycleway footway living_street motorway ... unknown

We have 27 different road types. Select step by step till you are satisfied:

roads_slct <- roads %>% 
  filter(fclass == "primary" | fclass == "secondary")

In queries we can use AND and OR operators: &and |.

In previous sessions we used geolayers of Georgia from The Humanitarian Data Exchange. But there are many different sources. For example you can download data directly from raster library:
Download country contour (level = 0):

gadm_0_georgia <- raster::getData('GADM', country = 'GEO', level = 0)

gadm_0_georgia_sf <- st_as_sf(gadm_0_georgia)

Plot country contour and roads layer:

ggplot()+
  geom_sf(data = gadm_0_georgia_sf)+
  geom_sf(data = roads_slct, col = "red")

Result is not very nice: the filtered road layer contains some obvious errors and unconnected roads.
Lets check the rivers layer:

rivers <-  st_read("gis_osm_waterways_free_1.shp")
## Reading layer `gis_osm_waterways_free_1' from data source `C:\ANTO\loengud\tbilisi\2020\web\gis_osm_waterways_free_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 17123 features and 5 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 39.52382 ymin: 40.93108 xmax: 46.8273 ymax: 45.74497
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
glimpse(rivers)
## Observations: 17,123
## Variables: 6
## $ osm_id   <fct> 17639034, 17810571, 17811047, 19854951, 19855058, 20449331...
## $ code     <int> 8103, 8101, 8101, 8102, 8101, 8101, 8101, 8101, 8101, 8101...
## $ fclass   <fct> canal, river, river, stream, river, river, river, river, r...
## $ width    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ name     <fct> <U+0425><U+043E><U+0431>, <U+10DB><U+10E2><U+10D9><U+10D5>...
## $ geometry <LINESTRING [arc_degree]> LINESTRING (41.76159 42.652..., LINEST...
# River types are stored in column 'fclass'. Get the list of unique records:
unique(rivers$fclass)
## [1] canal  river  stream drain 
## Levels: canal drain river stream
# filter:
rivers_slct <- rivers %>% filter(fclass == "river")

ggplot()+
  geom_sf(data = gadm_0_georgia_sf)+
  geom_sf(data = rivers_slct, col = "blue")

Result is again not very good. Clip rivers outside of Georgia. For that we can use st_intersection:

rivers_slct <- st_intersection(rivers_slct, gadm_0_georgia_sf)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
ggplot()+
  geom_sf(data = gadm_0_georgia_sf, col = "red")+
  geom_sf(data = rivers_slct, col = "blue", size= 0.25)

Result is still not 100% ok, but we don’t have time to deal with it. The quality of the dat ais just not good enough. Look for other options.

2.0.2 Roads layer from Natural Earth

Let’s take now roads layer. For that we should try some other option. Lot of spatial layers are available via Natural Earth

ne_roads  <- ne_download(scale = "large", type = "roads", category = "cultural")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Administrator\AppData\Local\Temp\Rtmpi83Afu", layer: "ne_10m_roads"
## with 56601 features
## It has 31 fields
## Integer64 fields read as strings:  scalerank question
class(ne_roads)
## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"

To plot Spatial***DataFrame objects we can use simply plot. It’s functionality is quite limited compared with ggplot2 but it is faster and able to plot much bigger objectS:

plot(ne_roads)

It is a global layer. How can we extract Georgian roads?

For that we can use function called gIntersection. It is a function for determining the intersection between the two given geometries. In other words - that portion of geometry A and geometry B that is shared between the two geometries.

ne_roads_GE <- gIntersection(ne_roads, gadm_0_georgia, byid = TRUE, drop_lower_td = TRUE)
class(ne_roads_GE)
## [1] "SpatialLines"
## attr(,"package")
## [1] "sp"
plot(gadm_0_georgia)
plot(ne_roads_GE, add= T, col = "red")

Result looks logical: all the road segments are connected with each other and nothing appears outside of Georgian borders.

2.0.3 Populated places

And now populated places. In OSM dataset the Georgian alphabet is used. Author of current tutorial needs to learn it… But till this time we have to find something else. Maybe the Natural Earth can help us again?

ne_populated_places <- ne_download(scale = "large", type = "populated_places", category = "cultural")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Administrator\AppData\Local\Temp\Rtmpi83Afu", layer: "ne_10m_populated_places"
## with 7343 features
## It has 119 fields
## Integer64 fields read as strings:  wof_id ne_id
class(ne_populated_places)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
glimpse(ne_populated_places@data)
## Observations: 7,343
## Variables: 119
## $ SCALERANK  <int> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...
## $ NATSCALE   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ LABELRANK  <int> 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5,...
## $ FEATURECLA <chr> "Admin-1 capital", "Admin-1 capital", "Admin-1 capital",...
## $ NAME       <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ NAMEPAR    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ NAMEALT    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ DIFFASCII  <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ NAMEASCII  <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ ADM0CAP    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ CAPIN      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ WORLDCITY  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MEGACITY   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ SOV0NAME   <chr> "Uruguay", "Uruguay", "Uruguay", "Uruguay", "Uruguay", "...
## $ SOV_A3     <chr> "URY", "URY", "URY", "URY", "URY", "TGO", "TGO", "TUN", ...
## $ ADM0NAME   <chr> "Uruguay", "Uruguay", "Uruguay", "Uruguay", "Uruguay", "...
## $ ADM0_A3    <chr> "URY", "URY", "URY", "URY", "URY", "TGO", "TGO", "TUN", ...
## $ ADM1NAME   <chr> "Colonia", "Flores", NA, "Canelones", "Florida", "Kara",...
## $ ISO_A2     <chr> "UY", "UY", "UY", "UY", "UY", "TG", "TG", "TN", "TN", "T...
## $ NOTE       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ LATITUDE   <dbl> -34.479999, -33.543999, -33.138999, -34.538004, -34.0990...
## $ LONGITUDE  <dbl> -57.8400025, -56.9009966, -58.3039975, -56.2840015, -56....
## $ CHANGED    <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,...
## $ NAMEDIFF   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ DIFFNOTE   <chr> "Added missing admin-1 capital. Population from GeoNames...
## $ POP_MAX    <int> 21714, 21093, 23279, 19698, 32234, 61845, 21054, 61705, ...
## $ POP_MIN    <int> 21714, 21093, 23279, 19698, 32234, 61845, 21054, 61705, ...
## $ POP_OTHER  <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ RANK_MAX   <int> 7, 7, 7, 6, 7, 8, 7, 8, 6, 8, 8, 8, 8, 7, 7, 7, 8, 6, 9,...
## $ RANK_MIN   <int> 7, 7, 7, 6, 7, 8, 7, 8, 6, 8, 8, 8, 8, 7, 7, 7, 7, 6, 9,...
## $ GEONAMEID  <dbl> 3443013, 3439749, 3442568, 3443413, 3442585, 2367568, 23...
## $ MEGANAME   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ LS_NAME    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ LS_MATCH   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ CHECKME    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_POP10  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_POP20  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_POP50  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_POP300 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_POP310 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_NATSCA <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_AREAKM <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_AREAKM <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_AREAMI <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_AREAMI <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_PERKM  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_PERKM  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_PERMI  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_PERMI  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_BBXMIN <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_BBXMIN <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_BBXMAX <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_BBXMAX <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_BBYMIN <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_BBYMIN <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MIN_BBYMAX <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MAX_BBYMAX <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MEAN_BBXC  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ MEAN_BBYC  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ COMPARE    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ GN_ASCII   <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ FEATURE_CL <chr> "P", "P", "P", "P", "P", "P", "P", NA, "P", "P", NA, "P"...
## $ FEATURE_CO <chr> "PPL", "PPL", "PPL", "PPL", "PPL", "PPL", "PPL", NA, "PP...
## $ ADMIN1_COD <dbl> 4, 6, 12, 2, 7, 5, 22, 0, 31, 34, 0, 6, 0, 33, 22, 0, 0,...
## $ GN_POP     <dbl> 21714, 21093, 23279, 19698, 32234, 61845, 21054, 0, 1987...
## $ ELEVATION  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ GTOPO30    <dbl> 28, 134, 43, 29, 74, 463, 374, 0, 49, 247, 0, 141, 0, 31...
## $ TIMEZONE   <chr> "America/Montevideo", "America/Montevideo", "America/Mon...
## $ GEONAMESNO <chr> "Geonames ascii name + lat.d + long.d matching.", "Geona...
## $ UN_FID     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ UN_ADM0    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ UN_LAT     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ UN_LONG    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1950    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1955    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1960    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1965    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1970    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1975    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1980    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1985    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1990    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP1995    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2000    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2005    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2010    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2015    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2020    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2025    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ POP2050    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ CITYALT    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ min_zoom   <dbl> 9.0, 9.0, 9.0, 9.0, 7.0, 9.0, 7.0, 6.1, 9.0, 9.0, 8.0, 9...
## $ wikidataid <chr> "Q56064", "Q862508", "Q849835", "Q841342", "Q842472", "Q...
## $ wof_id     <chr> "421199749", "890444639", "890451703", "890444649", "890...
## $ CAPALT     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ name_en    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_de    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_es    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_fr    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_pt    <chr> "Colónia do Sacramento", "Trinidad", "Fray Bentos", "Can...
## $ name_ru    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_zh    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ label      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_ar    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_bn    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_el    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_hi    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_hu    <chr> "Colonia del Sacramento", NA, NA, NA, NA, "Bassar", NA, ...
## $ name_id    <chr> NA, "Trinidad", "Fray Bentos", "Canelones", "Florida", N...
## $ name_it    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_ja    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_ko    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ name_nl    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_pl    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_sv    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_tr    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ name_vi    <chr> "Colonia del Sacramento", "Trinidad", "Fray Bentos", "Ca...
## $ wdid_score <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,...
## $ ne_id      <chr> "1159112629", "1159112647", "1159112663", "1159112679", ...

Downloaded object belongs to the class SpatialPointsDataFrame. Convert it to sf and filter locations in Georgia (iso_a2 code = GE).

# Here we can use 'filter' which works with sf-layers:
ne_populated_places_sf <- st_as_sf(ne_populated_places)
ne_populated_places_sf <- ne_populated_places_sf %>% filter(ISO_A2 =="GE")

ggplot()+
  geom_sf(data = gadm_0_georgia_sf)+
  geom_sf(data = ne_populated_places_sf, col = "red")+
  geom_sf_text(data = ne_populated_places_sf, aes(label = NAME), nudge_y = 0.1, nudge_x = 0.25)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Number of places is not very big, but adds some additional spatial context to the map.

2.0.4 Now get the elevation data!

For elevation we can use already familiar getDatafunction from library called raster. From there is available SRTM (Shuttle Radar Topography Mission) elevation data link. Instead of country ISO3 code we have to fix geographical coordinates of area we are interested.

srtm <- raster::getData('SRTM', lon = 44, lat = 42)

plot(srtm)
plot(gadm_0_georgia, add= T)

Elavation raster is download by tiles (tile which covers coordinates N42 and E44). from the plot we see, that it is not enough to cover whole Georgia. We have to donload one more tile

#Download one more tile:
srtm2 <- raster::getData('SRTM', lon = 46, lat = 42)

Downloaded rasters are store as different objects. We can bind them together:

#Mosaic/merge srtm tiles
srtmmosaic <- mosaic(srtm, srtm2, fun = mean)

plot(gadm_0_georgia)
plot(srtmmosaic, add = T)
plot(gadm_0_georgia, add= T)

In case we want to keep elevation data only for borders of Georgia, we can mask raster outside of the borders:

srtmmosaic_masked <- mask(srtmmosaic, gadm_0_georgia)

plot(gadm_0_georgia)
plot(srtmmosaic_masked, add = T)

Looks like we have now all the layers we need. And we can put all the layers in logical order together.

library(maptools) # enables function 'pointLabel'

plot(gadm_0_georgia, col = "grey")
plot(srtmmosaic_masked, add = T, col = (terrain.colors(20)))
plot(rivers_slct, add=T,  col="lightblue", lwd = 0.25)
plot(ne_roads_GE, add = T, col = "red", lwd = 0.1)
plot(ne_populated_places_sf, col = "black", add = T)

# Add place names
pointLabel(st_coordinates(ne_populated_places_sf),labels = ne_populated_places_sf$NAME)

2.1 Some examples of spatial analysis.

There are quite often the cases where we want to combine several spatial layers and derive new information from this combination.
Let us imagine the situation when we want to calculate average elevation above sea level for every region in Georgia.

# We need the regions layer in format of SpatialPolygonsDataFrame. Let's download it again:
gadm_1_sp <- raster::getData('GADM', country = 'GEO', level = 1)
class(gadm_1_sp)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
# extract mean from raster:
alt_GEO <- getData('alt', country='GEO', mask=TRUE)
region_Ave_Elev <- raster::extract(alt_GEO, gadm_1_sp, fun = mean, na.rm=TRUE, sp = T)

#Check the data:
region_Ave_Elev@data %>% 
  head()
##   GID_0  NAME_0   GID_1       NAME_1 VARNAME_1 NL_NAME_1                 TYPE_1
## 1   GEO Georgia GEO.1_1     Abkhazia   Sokhumi      <NA> Avtonomiuri Respublika
## 5   GEO Georgia GEO.2_1       Ajaria    Batumi      <NA> Avtonomiuri Respublika
## 6   GEO Georgia GEO.3_1        Guria  Ozurgeti      <NA>                 Region
## 7   GEO Georgia GEO.4_1      Imereti   Kutaisi      <NA>                 Region
## 8   GEO Georgia GEO.5_1      Kakheti    Telavi      <NA>                 Region
## 9   GEO Georgia GEO.6_1 Kvemo Kartli   Rustavi      <NA>                 Region
##             ENGTYPE_1 CC_1 HASC_1 GEO_msk_alt
## 1 Autonomous Republic <NA>  GE.AB   1050.9602
## 5 Autonomous Republic <NA>  GE.AJ   1142.7060
## 6              Region <NA>  GE.GU    605.3960
## 7              Region <NA>  GE.IM    681.9269
## 8              Region <NA>  GE.KA    914.1449
## 9              Region <NA>  GE.KK   1101.4145

Plot the results

# Convert SpatialpolygonDataFrame to sf:
region_Ave_Elev_sf <- st_as_sf(region_Ave_Elev)

ggplot()+
  geom_sf(data = region_Ave_Elev_sf, aes(fill = GEO_msk_alt), col = "grey", size= 0.25)+
  geom_sf_label(data = region_Ave_Elev_sf, aes(label= round(GEO_msk_alt, 0)), alpha = 0.5)+
  scale_fill_viridis_c()+
  labs(title = "Average elevation above sea level (m)", fill = "meters \n(a.s.l.)")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

You can save the previous map:

# Save previous plot:
ggsave("georgia_Mean_elev.png", dpi = 300, height = 4, width = 6, units = "in")

# Another posibility to save the map as raster image or pdf is:

png("filename.png", res= 300, units = "in", height = 4, width = 6)  # this will create an empty file and leaves it open for writing

ggplot()+ .........    #(the command, which creates the map)

 dev.off()     # This is very important! If you will not close the writing session all the next outputs will be written into "filename.png".

2.2 What is the warmest place in Georgia?

# download climate means (can be slow! (75MB)):
climate <- raster::getData('worldclim', var ='bio', res = 0.5, lon = 42, lat = 44)
climate_2 <- climate$bio1_17
climate_masked <- mask(climate_2, gadm_0_georgia)

# Temperature map:
climate_masked <- climate_masked / 10 # to get real values
plot(gadm_0_georgia)
plot(climate_masked, add= T, col = topo.colors(20))
plot(gadm_0_georgia, ad = T)

Good introduction to raster::getData: link

# we have some populated places:
ggplot()+
  geom_sf_text(data = ne_populated_places_sf, aes(label = NAME))
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

#sf to SpatialpointDataFrame:
ne_populated_places_sp <- as(ne_populated_places_sf, "Spatial")


# warmest place:
place_temp <- raster::extract(climate_masked, ne_populated_places_sp, fun = mean, na.rm=TRUE, sp = T)

library(knitr)
place_temp@data %>% 
  dplyr::select(NAME, bio1_17) %>% 
  arrange(-bio1_17) %>% 
  kable()
NAME bio1_17
Poti 14.4
Batumi 14.1
Sukhumi 14.1
Kutaisi 13.7
Rustavi 13.6
Tbilisi 13.1
Tskhinvali 8.6

Warmest place is Poti (14.4 Celsius degrees)!

2.3 Extra

2.3.1 On request: example of “calendar heatmap”

dat_tbilisi_airtemp <- read.csv2("http://aasa.ut.ee/tbilisi/tbilisi_airTemp_2017.csv")

dat_tbilisi_airtemp_day <- dat_tbilisi_airtemp %>% filter(!is.na(air_temp)) %>% mutate(date2 = as.Date(date)) %>% group_by(date2) %>% summarise(ave_airtemp = mean(air_temp)) %>% ungroup()

ggplot()+
  geom_line(data = dat_tbilisi_airtemp_day, aes(x= date2, y = ave_airtemp))

2.3.1.1 Create the calendar plot:

For that you need library “ggTimeseries”. To install this library you have to install first “devtools”. Also the installation process is a bit different:

install.packages("devtools")
library(devtools)
install_github("Ather-Energy/ggTimeSeries")

Plot:

library(ggTimeSeries)

ggplot_calendar_heatmap(dat_tbilisi_airtemp_day,
   'date2',
   'ave_airtemp')+
  scale_fill_viridis_c()+
  labs(title = "Mean daily air temperature in Tbilisi (2017)", fill= "C")

3 Homework number 4

Use previously downloaded Worldclim dataset. First example was based on annual mean air temperature. You task now is to create 2 thematic maps:

  1. Annual precipitation in Georgia,
  2. Average annual precipitation in Georgian regions.

Bind all the home works into one doc/pdf/html. Add script-files separately (in case of R-markdowm script is already in html file). Please add your name to each file name! E-mail results to: . Final deadline for all submissions is May 15th, 2020.

Possible results:

Author: Anto Aasa
in ISET, Tbilisi
Last update: 2020-04-22 14:33:04

.