7 Geodata Visualization
Give a look to this interesting tutorial (vignette) written by Edzer Pebesma about Plotting Simple Features.
7.1 Basic plot function
The most important R function to plot data is called plot(). In most cases, it will manage particular classes by invoking the plot function available in the correspondent package.
7.1.1 Numeric Objects
A simple data frame can be plotted using base R:
x = rnorm(100)
y = rnorm(100)
plot(x, y,
main = "Scatter Plot of Two Normally Distributed Random Variables",
xlab = "X values",
ylab = "Y values")
7.1.2 Data frame
Create the data frame reading a table from file:
library(readr)
soil_profiles <- read_table("datasets/Profile2.txt")
soil_profiles
#> # A tibble: 113 × 29
#> SOIL_PROFILE X_Coord Y_Coord clay d_tel slp_tel
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 052P1 2485212. 4566397. 27.4 215 21.4
#> 2 052P2 2484993. 4567244. 30.1 306 8.84
#> 3 052P3 2484157. 4568470. 34.7 320 9.01
#> 4 052P4 2484006. 4566499. 26.1 207 9.72
#> 5 052P5 2485743. 4567682. 24.3 425 24.3
#> 6 052P6 2489992. 4563739. 25.3 125 33.3
#> 7 052P7 2484854. 4562926. 22.7 55 0
#> 8 052P8 2488026. 4564857. 14.5 156 5.66
#> 9 052P9 2479296. 4566807. 29.5 226 1.98
#> 10 052P10 2476558. 4560723. 23.6 59 4.76
#> # ℹ 103 more rows
#> # ℹ 23 more variables: aspct_tel <dbl>,
#> # aspct_tel_linear <dbl>, t20b_T200_F15_gl <dbl>,
#> # clc2k_tel <dbl>, il_acv_tel <dbl>, il_gchan_tel <dbl>,
#> # il_gpeak_tel <dbl>, il_gpit_tel <dbl>,
#> # il_gplai_tel <dbl>, il_gridg_tel <dbl>,
#> # il_meanc_tel <dbl>, il_north_tel <dbl>, …A data frame can be plotted using base R:
plot( soil_profiles$d_tel,
soil_profiles$clay,
main = "Scatter plot of clay vs elevation (Valle Telesina)",
xlab = "Elevation [m]",
ylab = "Clay [%]" )
7.1.3 Simple Feature
Simple feature objects (sf) can be plotted directly with plot(). Different layers (geometry, attributes) are shown automatically in base plots:
library(sf)
Campania_region <- st_read("datasets/Ca.geojson")
#> Reading layer `Ca' from data source
#> `/home/giuliano/lectures/datasets/Ca.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 9 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 13.76238 ymin: 39.99052 xmax: 15.80645 ymax: 41.50835
#> Geodetic CRS: WGS 84
plot(Campania_region$geometry, main = "Campania region boundary")
7.1.4 Raster
Raster layers from the terra package can also be visualized using plot():
library(terra)
r <- rast("datasets/ASTGTMV003_N02E042_dem.tif")
plot(r, main = "Elevation raster for Somalia (N02,E042)")
Figure 7.1: Elevation raster for Somalia (N02,E042)
We can also merge rasters to make a mosaic:
library(terra)
r1 <- rast("datasets/ASTGTMV003_N01E041_dem.tif")
r2 <- rast("datasets/ASTGTMV003_N01E042_dem.tif")
r3 <- rast("datasets/ASTGTMV003_N02E041_dem.tif")
r4 <- rast("datasets/ASTGTMV003_N02E042_dem.tif")
merged <- mosaic(r1, r2, r3, r4)
plot(merged, main = "Merged DEM tiles for Somalia")
Figure 7.2: Merged DEM tiles for Somalia
7.2 Advanced ggplot2 library
The ggplot2 library provides more advanced and layered visualization options. With the help of geom_sf() or other geospatial-aware geoms, you can build elegant and customizable maps.
7.2.1 Example: Plotting sf data
library(ggplot2)
library(sf)
# Read the GeoJSON file
Campania_region <- st_read("datasets/Ca.geojson")
#> Reading layer `Ca' from data source
#> `/home/giuliano/lectures/datasets/Ca.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 9 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 13.76238 ymin: 39.99052 xmax: 15.80645 ymax: 41.50835
#> Geodetic CRS: WGS 84
# Plot using ggplot2
ggplot(data = Campania_region) +
geom_sf(fill = "lightblue", color = "black") +
theme_minimal() +
labs(title = "Campania Region Boundary")
7.2.2 Example: Plotting a raster (as points)
For raster data, convert it to a data frame first:
library(terra)
library(ggplot2)
# Load the raster
r <- rast("datasets/ASTGTMV003_N02E042_dem.tif")
# Convert to data frame for ggplot2
r_df <- as.data.frame(r, xy = TRUE)
# Get the name of the raster layer (to use in aes)
names(r)
names(r_df)[3] <- "elevation"
# Plot using ggplot2
ggplot(r_df) +
geom_raster(aes(x = x, y = y, fill = elevation)) +
scale_fill_gradientn(colours = terrain.colors(10)) +
theme_minimal() +
labs(
title = "Elevation Raster for Somalia (N02, E042) using ggplot.",
x = "Longitude",
y = "Latitude",
fill = "Elevation (m)"
)
Figure 7.3: Elevation Raster for Somalia (N02, E042) using ggplot.
7.3 Advanced tmaps library
The tmap library allows the creation of either a static plot (mode=“plot”) or a dynamic plot (mode=“view”) with pan, zoom and data selection functionalities.
The mechanism is very simple since it work by adding layers one on top of the previous. Each layer is made by a combination of 2 functions:
-
tm_shape()is used to select the geodata (vector or raster) to show -
tm_*is used to set the specifics of the graphical representation of the selected geodata
Examples for the second function are grouped by the data model or geometry type:
- POINT:
tm_dots()ortm_bubbles() - POLYGON:
tm_polygons()ortm_borders - RASTER:
tm_raster()
Many arguments can be configured to customize the result of the map produced.
7.3.2 Create a map for POINT geometry
7.3.2.1 Import data
Create Air temperature measurements simple feature:
t_day <- read_csv("AirTemperature_day.csv")
t_day_sf <- st_as_sf(t_day,coords = c("lon","lat"), crs=4326)Which variables are available?
names(t_day_sf)
#> [1] "Name" "elev" "mean" "min" "max"
#> [6] "N" "geometry"7.3.2.3 Bubbles
Let’s create a map using the tmap library:
tmap::tmap_mode("view")
tm_shape(t_day_sf) +
tm_bubbles(size="max",col="elev")7.3.3 Create a map for (MULTI)POLYGON geometry
7.3.3.1 Borders
# Read the GeoJSON file
Campania_region <- st_read("datasets/Ca.geojson",quiet=TRUE)
tmap::tmap_mode("view")
tm_shape(Campania_region) +
tm_borders(col="red",lwd=3)7.4 Animating views
library(plotly)
library(htmlwidgets)
library(gapminder)
p <- plot_ly(x = rnorm(100))
data(gapminder, package = "gapminder")
gg <- ggplot(gapminder, aes(gdpPercap, lifeExp, color = continent)) +
geom_point(aes(size = pop, frame = year, ids = country)) +
scale_x_log10()
ggplotly(gg)7.4.1 Note
Note that you learnt how to create interactive maps using the "view" mode.
In the book GEOG3915 GeoComputation and Spatial Analysis practicals by Lex Comber more examples can be found.