Orthophoto / satellite image as a base map

0.1 Orthophoto as base map

Start libraries:

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

options(scipen = 999)

In this session we continue with raster data. In previous session you saw that raster layers (e.g. elevation data) can be used to add spatial context to your thematic maps. Next to the elevation models, satellite images and orthophotos from various sources are available. And sometimes the orthophoto is very good alternative to the traditional base map.
In next example we look how to use orthophotos from Land Board Estonia as background image for your thematic map. The example itself is quite silly: we try to create map of postal offices and parcel machines in town Narva.
Possible workflow:

download relevant orthophotos from the web page of Land Board Estonia;
re-scale, merge and crop the images;
download locations of postal offices and parcel machines;
create the map.

First we download to your working directory relevant RGB orthophotos (you need 2 tiles: 65834 and 65932) from here: link

download

This is real life example! It means images are very big (resolution is 20 cm) and download takes time. But don’t worry the processing of images is even slower:).
Downloaded images are in zip-files. You can unzip them in File Folder but we can do it from RStudio as well:

unzip("65834_OF_RGB_GeoTIFF_2020_06_14.zip")

unzip("65932_OF_RGB_GeoTIFF_2020_06_14.zip")

It works if downloaded images are located in Working Directory (getwd()) of current workspace. In next step we import raster images to current workspace:

tif1 <- stack('65834.tif')
tif2 <- stack('65932.tif')

# plot:
plotRGB(tif1) # what happens when you plot it with plot()? Why?

Because the images are very big it is reasonable to reduce the resolution with rescale() (from 20 cm resolution to 2 m (factor = 10)):

For the first image aggregate() took 195 seconds… Be patient!

tif1 <- aggregate(tif1, fact = 10) %>% 
  stack()

tif2 <- aggregate(tif2, fact = 10) %>% 
  stack()

After time consuming operations it is reasonable to save the results:

#save rescaled rasters:
writeRaster(tif1, filename='65834_2m.tif', format="GTiff", overwrite = TRUE)
writeRaster(tif2, filename='65932_2m.tif', format="GTiff", overwrite = TRUE)

Now we can merge re-scaled images into one:

tiff_merged <- merge(tif1, tif2)

Plot the merged tif:

tm_shape(tiff_merged)+
    tm_rgb()+
    tm_grid()

We want to crop the orthophoto with border of Narva:

# download adnistrative borders:
download.file("https://geoportaal.maaamet.ee/docs/haldus_asustus/asustusyksus_shp.zip", destfile = "admin_borders.zip")

#unzip:
unzip("admin_borders.zip")

# import:
narva_contour <- st_read("asustusyksus_20221101.shp") # attention the date is changing every month!

## Reading layer `asustusyksus_20221101' from data source 
##   `C:\ANTO\loengud\geopythonR\rspatial_Git\rspatial\asustusyksus_20221101.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4713 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 369032.1 ymin: 6377141 xmax: 739152.8 ymax: 6634019
## Projected CRS: Estonian Coordinate System of 1997

# keep only town Narva:
narva_contour <- narva_contour %>% 
  filter(ANIMI == "Narva linn")

plot(st_geometry(narva_contour))
plotRGB(tif1, add=T)
plotRGB(tif2, add= T)

Crop and mask the tif by vector:

tiff_merged <- raster::crop(tiff_merged, narva_contour)
tiff_merged <- raster::mask(tiff_merged, narva_contour)

#Plot with tmap:
tm_shape(tiff_merged)+
  tm_rgb()+
tm_shape(narva_contour)+
  tm_borders("red")

Sometimes we want to used grey scale images:

tiff_merged_grey <- (tiff_merged$layer.1 + tiff_merged$layer.2 + tiff_merged$layer.3) / 3

# plot with cividis palette:
tm_shape(tiff_merged_grey)+
  tm_raster(palette = "cividis", legend.show = F)

# plot in grey scale:
colors <- gray.colors(255, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL)

tm_shape(tiff_merged_grey)+
  tm_raster(palette = c("grey10", "grey90"), legend.show = F)

You can save the image:

writeRaster(tiff_merged_grey, "trt_orto_gry.tif")

We have the base map!
In next step download locations of postal offices and parcel machines owned by Omniva:

post_locations <- read_delim("https://www.omniva.ee/locations.csv", 
                             ";", escape_double = FALSE, trim_ws = TRUE)

# check the structure:
glimpse(post_locations)

## Rows: 1,239
## Columns: 25
## $ ZIP                        <dbl> 96331, 10696, 96243, 69494, 96280, 96142, 9…
## $ NAME                       <chr> "1. eelistus minu.omniva.ee-s", "1it postip…
## $ TYPE                       <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1…
## $ A0_NAME                    <chr> "EE", "EE", "EE", "EE", "EE", "EE", "EE", "…
## $ A1_NAME                    <chr> "1. eelistus minu.omniva.ee-s", "Harju maak…
## $ A2_NAME                    <chr> "1. eelistus minu.omniva.ee-s", "Tallinn", …
## $ A3_NAME                    <chr> NA, "Kristiine linnaosa", "Abja-Paluoja lin…
## $ A4_NAME                    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ A5_NAME                    <chr> NA, "Sõpruse pst", "Pärnu mnt", "Pärnu mnt"…
## $ A6_NAME                    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ A7_NAME                    <chr> NA, "33", "13", "21", "8", "23", "77", "35"…
## $ A8_NAME                    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ X_COORDINATE               <dbl> 24.79337, 24.71429, 25.35581, 25.35393, 25.…
## $ Y_COORDINATE               <dbl> 59.43032, 59.42011, 58.12580, 58.12615, 59.…
## $ SERVICE_HOURS              <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TEMP_SERVICE_HOURS         <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TEMP_SERVICE_HOURS_UNTIL   <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TEMP_SERVICE_HOURS_2       <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TEMP_SERVICE_HOURS_2_UNTIL <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comment_est                <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comment_eng                <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comment_rus                <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comment_lav                <chr> NA, NA, NA, NA, NA, NA, NA, "Pa labi no T/C…
## $ comment_lit                <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Pa…
## $ MODIFIED                   <dttm> 2022-11-07 12:49:20, 2021-12-20 10:51:39, …

Coordinate columns are stored as character. Convert those columns to numeric (the comand is not in here…)!

post_locations <- post_locations %>% 
  mutate(X_COORDINATE = as.numeric(X_COORDINATE),
         Y_COORDINATE = as.numeric(Y_COORDINATE)) %>% 
  filter(!is.na(X_COORDINATE))

Convert it to sf:

post_locations_sf <- st_as_sf(post_locations, coords = c("X_COORDINATE", "Y_COORDINATE"), crs = 4326)

and plot:

tm_shape(tiff_merged_grey)+
  tm_raster(palette = c("grey10", "grey90"), legend.show = F)+
tm_shape(narva_contour)+
  tm_borders("navy")+
tm_shape(post_locations_sf)+
  tm_symbols(col = "cyan", size= 0.3,
          border.col = "blue") +
tm_layout(title = "post offices & parcel machines\nin Narva", title.color = "cyan",
          bg.color = "black")+
  tm_credits("A. Aasa")+
  tm_scale_bar(bg.color = "grey")

This was short walk through the raster data in R. You can read more about it from here: link

0.2 Base map for tmap

In previous examples we saw, that making of base maps by your self is not very hard, but it can be really time consuming. And the main thing we usually don’t have is time! In next example we look how to use maps from third party services as a base map for tmap. You are already familiar with possibility to switch tmap to view mode (tmap_mode("view")) and all the thematic layers are mapped on top of OSM,

By default the Esri World Gray Canvas is used:

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(post_locations_sf)+
  tm_dots("red")

You can set it to OSM:

tm_basemap(server = "OpenStreetMap")+
tm_shape(post_locations_sf)+
  tm_dots("red")

or to Esri World topo Map:

tm_basemap(server = "Esri.WorldTopoMap")+
tm_shape(post_locations_sf)+
  tm_dots("red")

Looks ok and it’s very fast!
But it’s interactive… And very often we need the map in static form. The png will do the job!
We can have it as well!

The help comes from library called basemaps:

library(basemaps)

# The extent of map we need:
bounding_box <- st_bbox(post_locations_sf)

#download the desired map:
bmap <- basemap_stars(ext = bounding_box,
                map_service = "osm",
                map_type = "streets",
                map_res = 2)

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'streets' from map service 'osm'...

# How big is it?
object.size(bmap)

## 14704360 bytes

# plot it:
tm_shape(bmap)+
  tm_rgb()

tmap_mode("plot")

## tmap mode set to plotting

tm_shape(bmap)+
  tm_rgb()+
tm_shape(post_locations_sf)+
  tm_dots("magenta")

You can test other options as well:

get_maptypes()

## $osm
## [1] "streets"     "streets_de"  "topographic"
## 
## $osm_stamen
## [1] "toner"      "toner_bg"   "toner_lite" "terrain"    "terrain_bg"
## [6] "watercolor"
## 
## $osm_thunderforest
##  [1] "cycle"          "transport"      "landscape"      "outdoors"      
##  [5] "transport_dark" "spinal"         "pioneer"        "mobile_atlas"  
##  [9] "neighbourhood"  "atlas"         
## 
## $carto
##  [1] "light"                "light_no_labels"      "light_only_labels"   
##  [4] "dark"                 "dark_no_labels"       "dark_only_labels"    
##  [7] "voyager"              "voyager_no_labels"    "voyager_only_labels" 
## [10] "voyager_labels_under"
## 
## $mapbox
## [1] "streets"   "outdoors"  "light"     "dark"      "satellite" "hybrid"   
## [7] "terrain"  
## 
## $esri
##  [1] "natgeo_world_map"                     
##  [2] "usa_topo_maps"                        
##  [3] "world_imagery"                        
##  [4] "world_physical_map"                   
##  [5] "world_shaded_relief"                  
##  [6] "world_street_map"                     
##  [7] "world_terrain_base"                   
##  [8] "world_topo_map"                       
##  [9] "world_dark_gray_base"                 
## [10] "world_dark_gray_reference"            
## [11] "world_light_gray_base"                
## [12] "world_light_gray_reference"           
## [13] "world_hillshade_dark"                 
## [14] "world_hillshade"                      
## [15] "world_ocean_base"                     
## [16] "world_ocean_reference"                
## [17] "antarctic_imagery"                    
## [18] "arctic_imagery"                       
## [19] "arctic_ocean_base"                    
## [20] "arctic_ocean_reference"               
## [21] "world_boundaries_and_places_alternate"
## [22] "world_boundaries_and_places"          
## [23] "world_reference_overlay"              
## [24] "world_transportation"                 
## [25] "delorme_world_base_map"               
## [26] "world_navigation_charts"

Some of them are working nicely, some are not available for current region and some of them are bad.

In next home work you have to demonstrate skills you learned during the current session!

Author: Anto Aasa
Supervisors: Anto Aasa & Lika Zhvania
LTOM.02.041
Last update: 2022-11-16 11:09:54

logo

Orthophoto / satellite image as a base map

Anto Aasa

November 18, 2022

0.1 Orthophoto as base map

0.2 Base map for tmap