Spatial interpolation

Libraries for current session

library(XML)
library(rvest)
library(tidyverse)
library(sf)
library(gplots)
library(lubridate)
library(ggspatial)
library(ggmap)
library(gstat)
library(knitr)
library(ggthemes)
library(tmap)
library(spatstat)
library(rgeos)

In this session we look how to interpolate the data spatially. Very often data is collected in discrete locations but the phenomena itself is continuous (air temperature, air pressure etc). This means that visualizing data only for observation sites may not be enough to get full picture about the spatial patterns.

Spatial interpolation of weather conditions in Estonia

Download latest weather observation data from Estonina Weather Service.
In previous sessions you have seen how to download data in various various formats (csv, Excel, shp, sdmx). One very common format is XML. From here you can read about the reasons: Why should I use XML?
XML format is not human-readable format. It is compromise between human-readability and machine-readability. And skills to use this format may be the first step from “ordinary people” to data scientist.

In R the import of XML is simple. Basically you don’t have to know anything about XML structure. weather data can be downloaded with 4 lines of code:

weatherDownTime <- Sys.time()
weather <- download.file("http://www.ilmateenistus.ee/ilma_andmed/xml/observations.php", "weather.xml")
xmlfile <- xmlParse("weather.xml")

# convert XML to data frame:
weather <- xmlToDataFrame(xmlfile)

# check the result:
glimpse(weather)

## Rows: 155
## Columns: 17
## $ name              <chr> "Kuressaare linn", "Tallinn-Harku", "Pakri", "Kunda"…
## $ wmocode           <chr> "", "26038", "26029", "26045", "26046", "26058", "26…
## $ longitude         <chr> "22.48944444411111", "24.602891666624284", "24.04008…
## $ latitude          <chr> "58.26416666666667", "59.398122222355134", "59.38950…
## $ phenomenon        <chr> "", "Overcast", "", "", "", "", "", "", "", "Overcas…
## $ visibility        <chr> "", "35.0", "31.0", "45.0", "20.0", "20.0", "20.0", …
## $ precipitations    <chr> "", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"…
## $ airpressure       <chr> "", "1007.3", "1006.5", "1008.8", "1009.7", "1010.1"…
## $ relativehumidity  <chr> "81", "81", "84", "81", "87", "84", "89", "81", "83"…
## $ airtemperature    <chr> "7.6", "7.7", "7.7", "7.1", "5.9", "6.3", "5.3", "5.…
## $ winddirection     <chr> "", "187", "169", "199", "197", "188", "209", "", "1…
## $ windspeed         <chr> "", "4.9", "4.5", "4.9", "4.7", "5.2", "5.6", "", "3…
## $ windspeedmax      <chr> "", "7.7", "8.6", "7.1", "6.9", "7.8", "7.9", "", "6…
## $ waterlevel        <chr> "", "", "", "", "", "", "", "", "64", "", "", "", ""…
## $ waterlevel_eh2000 <chr> "", "", "", "37", "", "", "", "", "", "", "", "", ""…
## $ watertemperature  <chr> "", "", "", "7.7", "", "", "", "", "", "", "", "", "…
## $ uvindex           <chr> "", "0.0", "", "", "", "", "", "", "", "", "", "", "…

Tidy up the table:

# convert columns from character to numeric:
weather <- weather %>% 
  mutate(longitude = as.numeric(as.character(longitude)), 
                              latitude = as.numeric(as.character(latitude)),
                              airtemperature = as.numeric(as.character(airtemperature)))

glimpse(weather)

## Rows: 155
## Columns: 17
## $ name              <chr> "Kuressaare linn", "Tallinn-Harku", "Pakri", "Kunda"…
## $ wmocode           <chr> "", "26038", "26029", "26045", "26046", "26058", "26…
## $ longitude         <dbl> 22.48944, 24.60289, 24.04008, 26.54140, 27.39827, 28…
## $ latitude          <dbl> 58.26417, 59.39812, 59.38950, 59.52141, 59.32902, 59…
## $ phenomenon        <chr> "", "Overcast", "", "", "", "", "", "", "", "Overcas…
## $ visibility        <chr> "", "35.0", "31.0", "45.0", "20.0", "20.0", "20.0", …
## $ precipitations    <chr> "", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"…
## $ airpressure       <chr> "", "1007.3", "1006.5", "1008.8", "1009.7", "1010.1"…
## $ relativehumidity  <chr> "81", "81", "84", "81", "87", "84", "89", "81", "83"…
## $ airtemperature    <dbl> 7.6, 7.7, 7.7, 7.1, 5.9, 6.3, 5.3, 5.9, 5.7, 5.7, 4.…
## $ winddirection     <chr> "", "187", "169", "199", "197", "188", "209", "", "1…
## $ windspeed         <chr> "", "4.9", "4.5", "4.9", "4.7", "5.2", "5.6", "", "3…
## $ windspeedmax      <chr> "", "7.7", "8.6", "7.1", "6.9", "7.8", "7.9", "", "6…
## $ waterlevel        <chr> "", "", "", "", "", "", "", "", "64", "", "", "", ""…
## $ waterlevel_eh2000 <chr> "", "", "", "37", "", "", "", "", "", "", "", "", ""…
## $ watertemperature  <chr> "", "", "", "7.7", "", "", "", "", "", "", "", "", "…
## $ uvindex           <chr> "", "0.0", "", "", "", "", "", "", "", "", "", "", "…

Map it:

ggplot()+
  geom_point(data = weather, aes(x = longitude, y=latitude, colour = airtemperature))+
  scale_color_gradientn(colours = topo.colors(20))+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime)+
  coord_fixed()

Stations with no data disturb the pattern. Filter them out:

airTemp_Df <- weather %>% 
  dplyr::select(name, longitude, latitude, airtemperature) %>% 
  filter(!is.na(airtemperature))

airTemp_Df <- weather %>% 
  dplyr::select(longitude, latitude, airtemperature) %>% 
  filter(!is.na(airtemperature))

Plot:

ggplot()+
  geom_point(data = airTemp_Df, aes(x = longitude, y=latitude, fill = airtemperature), shape=21, size=5)+
  scale_fill_gradientn(colours = c("navyblue", "blue", "cyan", "green", "yellow", "orange", "red"))+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime)

Spatial context is still almost missing. Download the contour of Estonia:

download.file("http://aasa.ut.ee/GisMaps/data/est_wgs84.zip", destfile = "est_wgs84.zip") 
unzip("est_wgs84.zip")

est_cont <- st_read("eestimaa_wgs84.shp")

## Reading layer `eestimaa_wgs84' from data source 
##   `C:\ANTO\loengud\geopythonR\rspatial_Git\rspatial\eestimaa_wgs84.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 21.76431 ymin: 57.50931 xmax: 28.209 ymax: 59.82202
## Geodetic CRS:  WGS 84

est_cont_3301 <- st_transform(est_cont, 3301)

# generalization:
est_cont_3301_smpl <- st_simplify(est_cont_3301, preserveTopology = F, dTolerance = 200)

Plot:

ggplot()+
  geom_sf(data = est_cont, fill="grey", col="grey30", size = 0.25)+
  geom_point(data = airTemp_Df, aes(x = longitude, y=latitude, fill = airtemperature), shape=21, size=5)+
  scale_fill_gradientn(colours = c("navyblue", "blue", "cyan", "green", "yellow", "orange", "red"))+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime,
       fill = "Temperature in \nCelsius degrees")

Nearest neighbour interpolation (Voronoi diagram)

Creating the voronoi diagram is spatial operation and assumes that point data is converted to spatial object:

# create simple feature: 
airTemp_sf <- st_as_sf(airTemp_Df, coords= c("longitude", "latitude"), crs = 4326)

# Convert from lat-lon (WGS84) to Estonian official CRS (LEST-97 [epsg:3301])
airTemp_sf_3301 <- st_transform(airTemp_sf, 3301)

Create the Voronoi diagramm:

# Combines geometries into one: 
airTemp_sf_3301_gm <- st_union(airTemp_sf_3301)

# Create the Voronoi diagram: 
airTemp_sf_voronoi <- st_voronoi(airTemp_sf_3301_gm)
airTemp_sf_voronoi <- st_sf(airTemp_sf_voronoi)
glimpse(airTemp_sf_voronoi)

## Rows: 1
## Columns: 1
## $ airTemp_sf_voronoi <GEOMETRYCOLLECTION [m]> GEOMETRYCOLLECTION (POLYGON...

ggplot()+
  geom_sf(data = airTemp_sf_voronoi, fill="red", colour  = "black", size=0.5)

Now we have a voronoi diagram of Estonian meteorological observation sites. But we have to mask it with Estonian contour and it is just an empty geometrical object, the attributes (air temperature is missing). There are several options to solve this problem. May be the easiest way is spatial join of voronoi diagram with spatial points layer:

# st_intersection masks the voronoi with Estonian contour: 
airTemp_sf_voronoi <- st_intersection(st_cast(airTemp_sf_voronoi), st_union(st_buffer(est_cont_3301_smpl, dist = 500))) # clip to smaller box
airTemp_sf_voronoi <- st_sf(airTemp_sf_voronoi)


ggplot()+
  theme_void()+
  geom_sf(data = airTemp_sf_voronoi)

# Spatial join:
airTemp_sf_voronoi2 <- st_join(airTemp_sf_voronoi, airTemp_sf_3301)
glimpse(airTemp_sf_voronoi2)

## Rows: 96
## Columns: 2
## $ airtemperature <dbl> 8.4, 8.5, 8.2, 7.8, 7.9, 7.3, 7.6, 7.6, 7.9, 7.9, 8.2, …
## $ geometry       <GEOMETRY [m]> MULTIPOLYGON (((412459.8 65..., MULTIPOLYGON (…

# plot the result:
ggplot()+
  geom_sf(data = airTemp_sf_voronoi2, aes(fill = airtemperature), size=0.25, colour = "grey70")+
  scale_fill_gradientn(colours = c("navyblue", "blue", "cyan", "green", "yellow", "orange", "red"))+
  geom_sf_label(data = airTemp_sf_3301, aes(label = airtemperature), col="black",  size=3, alpha = 0.5)+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime)

if the observation network is dense, then the result looks also nice and logical. But very often the Voronoi interpolation is not smooth enough and we have to consider some other methods.

IDW

According to the Wikipedia is inverse distance weighting (IDW) a type of deterministic method for multivariate interpolation with a known scattered set of points. The assigned values to unknown points are calculated with a weighted average of the values available at the known points.
The name given to this type of methods was motivated by the weighted average applied, since it resorts to the inverse of the distance to each known point (“amount of proximity”) when assigning weights.

# Get the extent of interpolation area:
grd_extent <- st_bbox(est_cont_3301_smpl)
 
# create grid:
x.range <- as.numeric(c(grd_extent[1], grd_extent[3])) # min/max longitude of the interpolation area
y.range <- as.numeric(c(grd_extent[2], grd_extent[4])) # min/max latitude of the interpolation area
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 500), y = seq(from = y.range[1], to = y.range[2], by = 500)) # expand points to grid

Convert grid to spatial object:

coordinates(grd) <- ~x + y 
gridded(grd) <- TRUE

We use previously created sf-layer:

glimpse(airTemp_sf_3301)

## Rows: 96
## Columns: 2
## $ airtemperature <dbl> 7.6, 7.7, 7.7, 7.1, 5.9, 6.3, 5.3, 5.9, 5.7, 5.7, 4.8, …
## $ geometry       <POINT [m]> POINT (411346.6 6459155), POINT (534250.5 6584619…

In current example we look how perform the interpolation inside of the Estonian contour (not for full rectangular grid). For that we have to convert spatial layers from simple feature (sf) object to the Spatial Polygon Data Frame. Then we convert Spatial Points Data Frame (SPDF) to ppp object. We also define the study area inside of Estonian contour:

glimpse(airTemp_sf_3301)

## Rows: 96
## Columns: 2
## $ airtemperature <dbl> 7.6, 7.7, 7.7, 7.1, 5.9, 6.3, 5.3, 5.9, 5.7, 5.7, 4.8, …
## $ geometry       <POINT [m]> POINT (411346.6 6459155), POINT (534250.5 6584619…

airTemp_sp_3301 <- as_Spatial(airTemp_sf_3301)

est_cont_3301_sp <- as_Spatial(est_cont_3301)

# grid:
idw_grid <- spsample(est_cont_3301_sp, type = "regular", n = 4000) # gridding may be time consuming! (cell size is currently in meters)
st_crs(airTemp_sf_3301)

## Coordinate Reference System:
##   User input: EPSG:3301 
##   wkt:
## PROJCRS["Estonian Coordinate System of 1997",
##     BASEGEOGCRS["EST97",
##         DATUM["Estonia 1997",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4180]],
##     CONVERSION["Estonian National Grid",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",57.5175539305556,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",24,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",58,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6375000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["northing (X)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (Y)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Topographic mapping (large scale)."],
##         AREA["Estonia - onshore and offshore."],
##         BBOX[57.52,20.37,60,28.2]],
##     ID["EPSG",3301]]

# idw
rslt_pt <- gstat::idw(airtemperature~1, airTemp_sp_3301, newdata = idw_grid, idp = 2)

## [inverse distance weighted interpolation]

class(rslt_pt)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

rslt_df2 <- as.data.frame(rslt_pt)
glimpse(rslt_df2)

## Rows: 3,990
## Columns: 4
## $ x1        <dbl> 650877.5, 654174.5, 650877.5, 654174.5, 657471.4, 697035.0, …
## $ x2        <dbl> 6380051, 6380051, 6383348, 6383348, 6383348, 6383348, 638664…
## $ var1.pred <dbl> 5.324659, 5.306541, 5.306980, 5.290717, 5.247140, 4.560104, …
## $ var1.var  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

rslt_grd_df <- rslt_df2 %>% filter(!is.na(var1.pred))
glimpse(rslt_grd_df)

## Rows: 3,990
## Columns: 4
## $ x1        <dbl> 650877.5, 654174.5, 650877.5, 654174.5, 657471.4, 697035.0, …
## $ x2        <dbl> 6380051, 6380051, 6383348, 6383348, 6383348, 6383348, 638664…
## $ var1.pred <dbl> 5.324659, 5.306541, 5.306980, 5.290717, 5.247140, 4.560104, …
## $ var1.var  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

# if you want to convert "SpatialPointsDataFrame" to "raster":
#library()
rslt_df2_rst <- rslt_df2 %>% 
  dplyr::select(-var1.var) %>% 
  raster::rasterFromXYZ()

class(rslt_df2_rst)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

plot(rslt_df2_rst)

# you can convert raster to contours:
rslt_df2_cnt <- raster::rasterToContour(rslt_df2_rst)
# plot contours:
plot(rslt_df2_cnt)

# Another option to create colour palette (from library 'gplots'):
palett <- rich.colors(40, palette = "temperature", rgb = F)

ggplot()+
  theme_map()+
  theme(panel.grid = element_line(color = "grey70", size = 0.2), legend.position = "right")+
  geom_tile(data=rslt_grd_df, aes(x = x1, y = x2, fill = var1.pred))+
  geom_sf(data = est_cont_3301_smpl, fill = NA, size = 0.3, colour = "grey45")+
  scale_fill_gradientn(colours = palett)+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime, fill = "Celsius")

Or you are bored of traditional raster maps and want to use other regular grid? Hexagonal binning for example? For the hexagonal binning it’s better to use Estonian official CRS (LEST97, epsg:3301). This means we have to transform temperature data. In case of SpatialPointDataFrame we can’t use sf command st_transform().
We have to use spTransform:

# definition of CRS lest97 (epsg:3301):
lest97 <- CRS("+proj=lcc +lat_1=59.33333333333334 +lat_2=58 +lat_0=57.51755393055556 +lon_0=24 +x_0=500000 +y_0=6375000 +ellps=GRS80 +units=m +no_defs")

# definition of WGS84 (epsg:4326):
wgs84 <- CRS("+proj=longlat +ellps=WGS84")

# Make the hexagonal grid (points)
grid_hex <- spsample(est_cont_3301_sp, type = "hexagonal", cellsize = 5000) # gridding may be time consuming!
plot(grid_hex)

# convert it to polygons:
grid_hex <- HexPoints2SpatialPolygons(grid_hex, dx = 5000)

#define CRS:
proj4string(grid_hex) <- lest97
plot(grid_hex)

# Convert temperature observatins from sf to SPDF:
airTemp_sp_3301 <- as(airTemp_sf_3301, "Spatial")

proj4string(airTemp_sp_3301) <- lest97

# calculate IDW for hexagons:

P_idw_hex <- gstat::idw(airtemperature~1, airTemp_sp_3301, newdata = grid_hex, idp = 2)

## [inverse distance weighted interpolation]

# convert SpatialPolygonDataFrame to sf:
rslt_hex <- st_as_sf(P_idw_hex)

# plot it:
ggplot()+
  theme_map()+
  theme(panel.grid = element_line(color = "grey60", size = 0.15),
        legend.position = "right")+
  geom_sf(data = rslt_hex, aes(fill = var1.pred), col = "grey20", size = 0.1)+
  scale_fill_viridis_c(option = "cividis")+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime,
       fill = "Celsius")

Polynomial fit

In current example we use 2nd order polynomial. How the spatial interpolation work is explained in here: link

# Estonian contour from sf to spdf (SpatialPolygonsDataFrame):
est_cont_sp <- as_Spatial(est_cont)
class(est_cont_sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

# layer overview:
summary(est_cont_sp)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##        min      max
## x 21.76431 28.20900
## y 57.50931 59.82202
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
##       FID   
##  Min.   :0  
##  1st Qu.:0  
##  Median :0  
##  Mean   :0  
##  3rd Qu.:0  
##  Max.   :0

# define the grid:
grd <- as.data.frame(spsample(est_cont_sp, "regular", n = 5000)) 
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

#grd_bu <- grd # backup the grid (we use it later again)

Add projection information to the empty grid:

airTemp_Df_xy <- as(airTemp_sf, "Spatial")

# set CRS:
proj4string(airTemp_Df_xy) <- wgs84
proj4string(grd) <- wgs84

# extent of spatial interpolation:
airTemp_Df_xy@bbox <- as_Spatial(est_cont)@bbox

# 2nd polynomial model:
f.2 <- as.formula(airtemperature ~ X + Y + I(X*X)+I(Y*Y) + I(X*Y))

# Add X and Y to airTemp_Df_xy:
airTemp_Df_xy$X <- coordinates(airTemp_Df_xy)[,1]
airTemp_Df_xy$Y <- coordinates(airTemp_Df_xy)[,2]

# Run the regression model
lm.2 <- lm(f.2, data = airTemp_Df_xy)

# Use the regression model output to interpolate the surface
dat.2nd <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.2, newdata=grd))) 

# Clip the interpolated raster
r   <- raster::raster(dat.2nd)
r.m <- raster::mask(r, as_Spatial(est_cont))

class(r.m)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

Plot the result:

plot(r.m)

To migrate the plotting of the result to ggplot2 we have to convert it to data frame:

r.m_df <- raster::as.data.frame(r.m, xy = T)
class(r.m_df)

## [1] "data.frame"

r.m_df <- r.m_df %>% 
  filter(!is.na(var1.pred))

glimpse(r.m_df)

## Rows: 5,005
## Columns: 3
## $ x         <dbl> 25.50854, 25.69196, 25.72864, 25.50854, 25.72864, 25.76532, …
## $ y         <dbl> 59.64594, 59.64594, 59.64594, 59.60926, 59.60926, 59.60926, …
## $ var1.pred <dbl> 7.403732, 7.348715, 7.337483, 7.368508, 7.299827, 7.288113, …

# plot:
ggplot()+
  geom_tile(data = r.m_df, aes(x=x, y=y, fill=var1.pred))+
  scale_fill_gradientn(colours = c("navy", "darkgreen", "yellow", "orange", "red"))+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime)

In real life the colour palette should be more modest…

Sometimes we want to use contours instead of raster:

pr_cont <- raster::rasterToContour(r.m)
class(pr_cont)

## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"

plot(pr_cont)

Plot interpolated contours:

# convert to sf:
pr_cont <- st_as_sf(pr_cont)

# set CRS:
st_crs(pr_cont) <- 4326

# convert contour labels from character to numeric:
pr_cont <- pr_cont %>% 
  mutate(level = as.numeric(as.character(level)))

glimpse(pr_cont)

## Rows: 9
## Columns: 2
## $ level    <dbl> 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0
## $ geometry <MULTILINESTRING [°]> MULTILINESTRING ((27.24996 ..., MULTILINESTRING ((26.79366 ..…

# plot
ggplot()+
  theme_minimal()+
  geom_sf(data = est_cont, size=0.1)+
  geom_sf(data = pr_cont, aes(colour = level), size = 0.5)+
  geom_sf_label(data = pr_cont, aes(label = level), size = 2)+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime)+
  scale_colour_gradientn(colours = palett)

3rd order polynomial fit

# 3rd polynomial model:
f.3 <- as.formula(airtemperature ~ X + Y + I(X^X)+I(Y^Y) + I(X^Y))

# Run the regression model
lm.3 <- lm(f.3, data = airTemp_Df_xy)

# Use the regression model output to interpolate the surface
dat.3nd <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.3, newdata = grd))) 

# Clip the interpolated raster
r3   <- raster::raster(dat.3nd)
r3.m <- raster::mask(r3, as_Spatial(est_cont))

plot(r3.m)

r3.m_df <- raster::as.data.frame(r3.m, xy = T)
class(r3.m_df)

## [1] "data.frame"

r3.m_df <- r3.m_df %>% 
  filter(!is.na(var1.pred))

glimpse(r3.m_df)

## Rows: 5,005
## Columns: 3
## $ x         <dbl> 25.50854, 25.69196, 25.72864, 25.50854, 25.72864, 25.76532, …
## $ y         <dbl> 59.64594, 59.64594, 59.64594, 59.60926, 59.60926, 59.60926, …
## $ var1.pred <dbl> 7.566469, 7.462625, 7.442627, 7.449888, 7.324069, 7.303943, …

# plot:
ggplot()+
  geom_tile(data = r3.m_df, aes(x=x, y=y, fill = round(var1.pred, 1)))+
  scale_fill_gradientn(colours = palett)+
  labs(title = "Air temperature in Estonia (3rd order polynomial)", subtitle = weatherDownTime, fill = "temperature")

Still too general…

Kriging

Accordign to Wikipedia: “In statistics, originally in geostatistics, kriging or Gaussian process regression is a method of interpolation for which the interpolated values are modeled by a Gaussian process governed by prior covariances. Under suitable assumptions on the priors, kriging gives the best linear unbiased prediction of the intermediate values. Interpolating methods based on other criteria such as smoothness (e.g., smoothing spline) need not yield the most likely intermediate values. The method is widely used in the domain of spatial analysis and computer experiments.”

Kriging is a little more involved than IDW as it requires the construction of a semivariogram model to describe the spatial autocorrelation pattern for your particular variable.

More about kriging (and other methods): link

# define the model:
P_gs <- gstat::gstat(formula = airtemperature ~ 1, locations = airTemp_Df_xy)

P_v <- variogram(P_gs, width = 1)
plot(P_v)

P_semivariog <- variogram(airtemperature ~ 1, locations = airTemp_Df_xy, data = airTemp_Df_xy)
plot(P_semivariog)

P_semivariog %>% kable()

np	dist	gamma	dir.hor	dir.ver	id
20	4.65979	0.0420000	0	0	var1
49	14.82298	0.1469388	0	0	var1
86	24.19556	0.1777326	0	0	var1
105	33.91800	0.2413333	0	0	var1
153	43.43339	0.2683660	0	0	var1
175	52.89657	0.3293714	0	0	var1
176	62.15702	0.3415341	0	0	var1
175	72.22992	0.3913714	0	0	var1
218	81.55998	0.5301147	0	0	var1
222	91.31620	0.4996622	0	0	var1
238	100.82620	0.5541176	0	0	var1
229	110.62807	0.7344323	0	0	var1
219	119.95181	0.7872831	0	0	var1
235	129.74137	0.8640426	0	0	var1
247	139.43843	1.0627328	0	0	var1

From the empirical semivariogram plot and the information contained in the semivariog gstat object, we can estimate the sill, range and nugget to use in our model semivariogram.

“https://mgimond.github.io/Spatial/spatial-interpolation.html”

In this case, the range (the point on the distance axis where the semivariogram starts to level off) is around the value of the last lag - 139.4384272 - so we’ll use Range = 140.

The Sill (the point on the y axis where the semivariogram starts to level off) is around 1.0627328.

The nugget looks to be around 0.042.

Using this information we’ll generate a model semivariogram using the vgm() function in gstat:

# currently "Gau" model (to see others: ?vgm)
psill <- P_semivariog$gamma[nrow(P_semivariog)]
nugget <- P_semivariog$gamma[1]
range <- ceiling(P_semivariog$dist[nrow(P_semivariog)])

P_model.variog <- vgm(psill = psill, model = "Mat", nugget = nugget, range = range)

# fit model:
P_fit.variog <- fit.variogram(P_semivariog, P_model.variog)
plot(P_semivariog, P_fit.variog) # check from Help what are the other available models and choose the best fit!

Kriging:

P_krig <- gstat::krige(formula = airtemperature ~ 1, locations = airTemp_Df_xy, newdata = grd, model = P_model.variog)

## [using ordinary kriging]

#plot(P_krig)

P_r <- raster::raster(P_krig)
P_r.m <- raster::mask(P_r, est_cont)

P_r.m_df <- raster::as.data.frame(P_r.m, xy=T)
P_r.m_df <- P_r.m_df %>% dplyr::filter(!is.na(var1.pred))
glimpse(P_r.m_df)

## Rows: 5,005
## Columns: 3
## $ x         <dbl> 25.50854, 25.69196, 25.72864, 25.50854, 25.72864, 25.76532, …
## $ y         <dbl> 59.64594, 59.64594, 59.64594, 59.60926, 59.60926, 59.60926, …
## $ var1.pred <dbl> 7.445071, 7.397589, 7.372434, 7.343304, 7.314017, 7.259391, …

ggplot()+
  theme_minimal()+
  geom_tile(data = P_r.m_df, aes(x = x, y = y, fill = var1.pred))+
  geom_sf(data = est_cont, fill=NA, size=0.3, colour="grey45")+
  scale_fill_gradientn(colours = c("navy", "purple",  "yellow", "orange", "red1"))+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime)

And finally put raster, isotherms and labels together:

# convert to raster
P_r.m_cont <- raster::rasterToContour(P_r.m)

# convert to sf
P_r.m_cont <- st_as_sf(P_r.m_cont)

# set CRS:
st_crs(P_r.m_cont) <- 4326

## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that

# convert contour labels from character to numeric:
P_r.m_cont <- P_r.m_cont %>% 
  mutate(level = as.numeric(as.character(level)))

glimpse(P_r.m_cont)

## Rows: 9
## Columns: 2
## $ level    <dbl> 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0
## $ geometry <MULTILINESTRING [°]> MULTILINESTRING ((27.72157 ..., MULTILINESTRING ((26.70116 ..…

# plot
ggplot()+
  theme_minimal()+
  geom_tile(data = P_r.m_df, aes(x = x, y = y, fill = var1.pred))+
  scale_fill_gradientn(colours = palett)+
  geom_sf(data = P_r.m_cont, colour = "grey60", size = 0.5)+
  geom_sf(data = est_cont, size = 0.1, fill = NA)+
  geom_sf_label(data = P_r.m_cont, aes(label = level), size = 2, alpha = 0.75)+
  labs(title = "Air temperature in Estonia", subtitle = weatherDownTime,
       fill = "Celsius")+
  scale_colour_gradientn(colours = palett)

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Are you satisfied with this? Sometimes it is very hard to find suitable interpolation method. In that case you can try to apply some simple statistical model. Let’s test how it could be done with multiple linear regression. Look again the air temperature observations:

# air temperature observations:
glimpse(airTemp_Df)

## Rows: 96
## Columns: 3
## $ longitude      <dbl> 22.48944, 24.60289, 24.04008, 26.54140, 27.39827, 28.10…
## $ latitude       <dbl> 58.26417, 59.39812, 59.38950, 59.52141, 59.32902, 59.38…
## $ airtemperature <dbl> 7.7, 7.9, 7.9, 7.2, 5.8, 6.3, 5.3, 5.9, 5.7, 5.8, 4.9, …

# hexagons layer (convert it to sf):
class(grid_hex)

## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"

grid_hex_sf <- st_as_sf(grid_hex)

# change the CRS (station coordinates are in WGS84):
grid_hex_sf_4326 <- st_transform(grid_hex_sf, 4326)

Before we run the model, we should check, is there correlation between air temperature and coordinates:

ggplot(data = airTemp_Df)+
  geom_point(aes(x= longitude, y = airtemperature))+
  stat_smooth(aes(x= longitude, y = airtemperature), method = "lm")

ggplot(data = airTemp_Df)+
  geom_point(aes(x= latitude, y = airtemperature))+
  stat_smooth(aes(x= latitude, y = airtemperature), method = "lm")

It is! Probably not very strong, but still. Create the multiple linear regression model, where air temperature is function and coordinates are arguments:

airtemperature_lm <- lm(airtemperature~longitude + latitude, data= airTemp_Df)

Check the model. Is it significant (P-value < 0.05)? Are all model members significant?:

# model summary:
summary(airtemperature_lm)

## 
## Call:
## lm(formula = airtemperature ~ longitude + latitude, data = airTemp_Df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.47277 -0.28908  0.04826  0.28936  1.10745 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -27.56380    5.34822  -5.154 1.43e-06 ***
## longitude    -0.51453    0.03073 -16.743  < 2e-16 ***
## latitude      0.80879    0.09132   8.857 5.35e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4837 on 93 degrees of freedom
## Multiple R-squared:  0.7826, Adjusted R-squared:  0.7779 
## F-statistic: 167.4 on 2 and 93 DF,  p-value: < 2.2e-16

# model coefficients:
coef(airtemperature_lm)

## (Intercept)   longitude    latitude 
## -27.5638009  -0.5145256   0.8087903

Model is significant. We can now calculate values for all the hexagons:

grid_hex_crd_4326 <- st_centroid(grid_hex_sf_4326) %>% 
  st_coordinates() %>% 
  as_tibble()

glimpse(grid_hex_crd_4326)

## Rows: 2,013
## Columns: 2
## $ X <dbl> 26.54872, 26.63216, 26.34258, 26.42614, 26.50969, 26.59324, 26.67677…
## $ Y <dbl> 57.55581, 57.55408, 57.59875, 57.59716, 57.59551, 57.59380, 57.59204…

grid_hex_sf_4326 <- bind_cols(grid_hex_sf_4326, grid_hex_crd_4326)


# calculate predictions:
grid_hex_sf_4326 <- grid_hex_sf_4326 %>% 
  mutate(airtemp_mod= X * coef(airtemperature_lm)[2] +
           Y * coef(airtemperature_lm)[3] +
           coef(airtemperature_lm)[1])

Plot the results:

tmap_mode("plot")

tm_shape(grid_hex_sf_4326)+
  tm_polygons("airtemp_mod", style = "pretty", palette= "-Spectral")

Not very nice? As a geographers we now that air temperature decreases by 0.6 Celsius degrees every 100 meters. May be we should add the elevation component to our model as well? Let’s try!

# download again elevation data:
est_elev <- raster::getData('alt', country = 'EST', mask = F) 
plot(est_elev)

# airtemp sf :
glimpse(airTemp_sf)

## Rows: 96
## Columns: 2
## $ airtemperature <dbl> 7.7, 7.9, 7.9, 7.2, 5.8, 6.3, 5.3, 5.9, 5.7, 5.8, 4.9, …
## $ geometry       <POINT [°]> POINT (22.48944 58.26417), POINT (24.60289 59.398…

# extract the elevation for every meterological station:
airTemp_sf_elev <- raster::extract(est_elev, airTemp_sf)
airTemp_sf_elev

##  [1]   7  34  17   4  75  29 120  34  70  69  79  59  84  63  53  24   5  NA  10
## [20]  NA   1   7   3  86   3   6   2   2  NA  NA  NA   3  NA  NA  43  11  12  57
## [39]  72  62  29  41  49  31  46  52  46  55 102  38  32  17  45  39  35  39  48
## [58]  55  35  32  45  48  40  64  85  14  74  29  13   8  22  40  17  28  22  27
## [77]  30  65  70  73  62  55  11  14  13  55  32  28  33  31  33  30  23  77 139
## [96] 233

airTemp_sf$elev <- airTemp_sf_elev

# plot it:
tm_shape(est_cont)+
  tm_polygons()+
tm_shape(airTemp_sf)+
  tm_dots("elev", size = 0.2, palette = "cividis", colorNA = "red")

## For some coastal stations the elevation model doesn't fit. We may assume, that elevation for those station is quite close to zero;
# replace NA with 0:
airTemp_sf$elev[is.na(airTemp_sf$elev)] <- 0

# plot it!
tm_shape(est_cont)+
  tm_polygons()+
tm_shape(airTemp_sf)+
  tm_dots("elev", size = 0.2, palette = "cividis", colorNA = "red")

#Better!

# extract average elevation to every hexagon:
grid_hex_sf_4326 <- raster::extract(est_elev, grid_hex_sf_4326, fun = mean, na.rm = TRUE, sp = T) # takes some time...

class(grid_hex_sf_4326)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

# convert it back to sf:
grid_hex_sf_4326 <- st_as_sf(grid_hex_sf_4326)

# plot elevation map:
tm_shape(grid_hex_sf_4326)+
  tm_fill("EST_alt")

# prepare data for regression model:
airTemp_df2 <- airTemp_sf %>% 
  st_drop_geometry()

airTemp_crd2 <- airTemp_sf %>% 
  st_coordinates() %>% 
  as_tibble()

airTemp_df2 <- bind_cols(airTemp_df2, airTemp_crd2)

#create the model:
airtemperature_lm2 <- lm(airtemperature ~ X + Y + elev, data= airTemp_df2)

summary(airtemperature_lm2)

## 
## Call:
## lm(formula = airtemperature ~ X + Y + elev, data = airTemp_df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23840 -0.14400 -0.02464  0.22595  0.96101 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -20.843805   4.063215  -5.130 1.60e-06 ***
## X            -0.387708   0.027186 -14.261  < 2e-16 ***
## Y             0.646653   0.070623   9.156 1.36e-14 ***
## elev         -0.011108   0.001281  -8.672 1.42e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3607 on 92 degrees of freedom
## Multiple R-squared:  0.8804, Adjusted R-squared:  0.8765 
## F-statistic: 225.7 on 3 and 92 DF,  p-value: < 2.2e-16

coef(airtemperature_lm2)

##  (Intercept)            X            Y         elev 
## -20.84380454  -0.38770753   0.64665299  -0.01110797

# calculate prediction for every polygon:
grid_hex_sf_4326 <- grid_hex_sf_4326 %>% 
  mutate(airtemp_mod_elev = X * coef(airtemperature_lm2)[2] +
           Y * coef(airtemperature_lm2)[3] +
           EST_alt  * coef(airtemperature_lm2)[4]+
           coef(airtemperature_lm2)[1])


# final plot:
tm_grid(col = "grey")+
tm_shape(grid_hex_sf_4326)+
  tm_polygons("airtemp_mod_elev", style = "pretty", palette = c("dodgerblue", "yellow", "orange"), title = "C")+
  tm_layout(title = "Air temperature in Estonia")+
  tm_scale_bar(position = c("left", "bottom"))+
  tm_credits("A. Aasa")

We have to be careful! Is it still the map of air temperature of is it map of elevation?

Heat maps

Wikipedia: “A heat map (or heatmap) is a graphical representation of data where the individual values contained in a matrix are represented as colors. A density function visualization is a heat map for representing the density of dots in a map.”

Transparent points

Sometimes we even don’t need complex time consuming calculations and just using the transparency works well. In our next example we use gps-points from Monterey city, California. For that we have to download monterey map from Stamen maps.

# Download GPS points for Moneterey:
url <-  "http://aasa.ut.ee/Rspatial/data/gps_us_monterey.zip" # define the dataset address
download.file(url, "gps_us_monterey.zip") # download file
unzip("gps_us_monterey.zip") # unzip files
file.remove("gps_us_monterey.zip") # tidy up by removing the zip file

## [1] TRUE

# read shp:
gpsData_monterey <- st_read("gps_us_monterey.shp")

## Reading layer `gps_us_monterey' from data source 
##   `C:\ANTO\loengud\geopythonR\rspatial_Git\rspatial\gps_us_monterey.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 11373 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -121.9415 ymin: 36.57022 xmax: -121.8795 ymax: 36.62333
## Geodetic CRS:  WGS 84

# convert shp to data frame:
gpsData_monterey_df <- as_data_frame(gpsData_monterey)

gpsData_monterey_crd <- st_coordinates(gpsData_monterey) %>% 
  as_data_frame()

gpsData_monterey_df <- gpsData_monterey_df %>% 
  dplyr::select(-geometry)

gpsData_monterey_df <- bind_cols(gpsData_monterey_df, gpsData_monterey_crd)

Map (OSM) download:

# calculate & define bounding box: 
base_bbox <- st_bbox(gpsData_monterey) %>% as.numeric()

# it can be defined as an independednt object:
bbox <- c(left = base_bbox[1], 
          bottom = base_bbox[2], 
          right = base_bbox[3],
          top = base_bbox[4])

# download the map:
monterey_stamen_map <- get_stamenmap(bbox, maptype = "toner-lite", zoom = 14, force = T)

#plot it:
ggmap(monterey_stamen_map, extent = "device")

# plot with points
ggmap(monterey_stamen_map, extent = "device")+
  geom_sf(data = gpsData_monterey, size = 0.5, colour = "red", inherit.aes = F)

Or you can plot it with geom_point. ‘alpha’-parameter defines the transparency (0 - transparent; 1 - not transparent)

ggmap(monterey_stamen_map, extent = "device")+
  geom_point(data = gpsData_monterey_df, aes(x=X, y=Y), size = 1, colour = "red", alpha = 0.01)

map_extent <- st_as_sfc(st_bbox(gpsData_monterey))

monterey_hex_sf <- st_make_grid(map_extent, cellsize = .001, square = FALSE)

monterey_hex_sf <- st_sf(monterey_hex_sf)
class(monterey_hex_sf)

## [1] "sf"         "data.frame"

glimpse(monterey_hex_sf)

## Rows: 3,937
## Columns: 1
## $ monterey_hex_sf <POLYGON [°]> POLYGON ((-121.942 36.57051..., POLYGON ((-121…

monterey_hex_sf_rows <- nrow(monterey_hex_sf)

# create id-column:
monterey_hex_sf <- monterey_hex_sf %>% 
  mutate(hex_id = seq(1, monterey_hex_sf_rows, 1))

gpsData_monterey_hex <- st_join(gpsData_monterey, monterey_hex_sf)

gpsData_monterey_hex_df <- gpsData_monterey_hex %>% 
  st_drop_geometry() %>% 
  group_by(hex_id) %>% 
  summarise(n = n()) %>% 
  ungroup()


monterey_hex_sf <- left_join(monterey_hex_sf, gpsData_monterey_hex_df)

## Joining, by = "hex_id"

monterey_hex_sf <- monterey_hex_sf %>% 
  filter(!is.na(n))


ggmap(monterey_stamen_map, extent = "device")+
  geom_sf(data = monterey_hex_sf, aes(fill = log10(n)), inherit.aes = F, alpha = 0.75)+
  scale_fill_gradientn(colours = c("pink", "red", "red4"))

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Density

We can produce density map just as graph:

ggmap(monterey_stamen_map, extent = "device", darken = 0.5)+
  theme_void() +
  stat_density2d(data = gpsData_monterey_df, aes(x = X, y = Y, fill = ..level.., alpha = ..level..), geom = "polygon", bins = 30) +
  scale_fill_gradientn(colours = c("yellow", "orange", "firebrick2")) +
  scale_alpha(range = c(0.2, 1))

Kernel Density Estimates (KDE)

In some cases we don’t need only the visualization but also the raster surface for further analysis.

In this case we can use the functionality from library adehabitatHR:

library(adehabitatHR)

kde_dens <- kernelUD(as_Spatial(gpsData_monterey), h = "href", grid = 1000)
class(kde_dens)

## [1] "estUD"
## attr(,"package")
## [1] "adehabitatHR"

What does the argument h or grid?
Plot the result:

plot(kde_dens)

So plotted raster is hardly readable, because the spatial context is missing. We could delete or filter out the pixels where value is 0. Output is stored in specific estUD format. For the further analysis the raster class is more relevant. Convert it:

kde_dens_rst <- raster::raster(kde_dens)

Get the minimum and maximum value of raster:

raster::minValue(kde_dens_rst)

## [1] 0

raster::maxValue(kde_dens_rst)

## [1] 14343.75

Now we can plot the raster with defined range of values:

plot(kde_dens, zlim = c(.1, raster::maxValue(kde_dens_rst)))

Or we could extract the raster area were 99% of values are. The aim of the command getverticeshr is to calculate the home-range area of animals, but can be used in current case as well:

range99 <- getverticeshr(kde_dens, percent = 99)
class(range99)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

range99 <- range99 %>% 
  st_as_sf()

Sometimes the created sf-object is invalid. But it can be “fixed”. In current case we can use:

tmap_options(check.and.fix = T)

Read from help, what it does?
Or command st_make_valid()?

Now we can plot the full raster and area, where 99% of values are stored:

tm_shape(kde_dens)+
  tm_raster(palette = "cividis", style = "fisher", title = "Density")+
tm_shape(range99)+
  tm_borders(alpha=.7, col = "magenta", lwd = 2)+
  tm_fill(alpha = .5, col = "magenta")+
tm_layout(legend.outside = T)

## Warning: The shape range99 is invalid. See sf::st_is_valid

We crop and mask the raster with range99 polygon:

kde_dens_rst_msk <- raster::crop(kde_dens_rst, range99)

tm_shape(kde_dens_rst_msk)+
  tm_raster(palette = "cividis", style = "cont", title = "Density")+
  tm_layout(legend.outside = T)

and mask:

kde_dens_rst_msk <- raster::mask(kde_dens_rst_msk, range99)

tm_shape(kde_dens_rst_msk)+
  tm_raster(palette = "-cividis", style = "cont", title = "Density")+
  tm_layout(legend.position = c("left", "bottom"), title = "density of GPS points", title.position = c("right", "top"))

Still the spatial context is missing. Add basemap:

library(basemaps)

# The extent of map we need:
kde_dens_rst_msk_bbox <- st_bbox(kde_dens_rst_msk)

#download the desired map:
kdens_bmap <- basemap_stars(ext = kde_dens_rst_msk_bbox,
                map_service = "esri",
                map_type = "natgeo_world_map",
                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 'natgeo_world_map' from map service 'esri'...

tm_shape(kdens_bmap)+
  tm_rgb(alpha = .7)+
  tm_shape(kde_dens_rst_msk)+
  tm_raster(palette = c("pink", "red1"), style = "fisher", title = "Density", alpha = .5)+
  tm_layout(legend.position = c("left", "bottom"), title = "density of GPS points", title.position = c("right", "top"))

There are other map types available:

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"

Sometimes plotting the map in 3D can help to improve the impression:

raster::persp(kde_dens_rst, box = F, maxpixels = 1000)

Obviously it’s currently not the case…

Or the other option to delete raster pixels where value is > 0.0001:

kde_dens3 <- kde_dens_rst
kde_dens3[kde_dens3 < 0.0001] <- NA


kde_dens3 <- raster::mask(kde_dens_rst, kde_dens3)

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(kde_dens3)+
  tm_raster(style = "log10", alpha = 0.6,
            palette = c("snow", "yellow", "orange", "red4", "black"))

---

Author: Anto Aasa
Supervisors: Anto Aasa & Lika Zhvania
LTOM.02.041
Last update: 2022-11-07 17:10:10

.