Packages

In this post we will use the following libraries:

#install the libraries if necessary
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("sf")) install.packages("sf")
if(!require("raster")) install.packages("raster")
if(!require("rnaturalearth")) install.packages("rnaturalearth")

#packages
library(rnaturalearth)
library(sf)
library(raster)
library(tidyverse)
library(RColorBrewer)

The coast of Iceland as an example

Our example in this post will be Iceland, and, as it is an island territory it will facilitate the tutorial showing the process in a simple manner. The rnaturalearth package allows you to import the boundaries of countries (with different administrative levels) from around the world. The data comes from the platform naturalearthdata.com. I recommend exploring the package, more info here. The ne_countries( ) function imports the country boundaries. In this case we indicate with the argument scale the resolution (10, 50 or 110m), with country we indicate the specific country of interest and with returnclass we determine which class we want (sf or sp), in our case sf (simple feature).

world <- ne_countries(scale=50) #world map with 50m resolution

plot(world) #sp class by default

#import the limits of Iceland
iceland <- ne_countries(scale=10,country = "Iceland",returnclass = "sf")

#info of our spatial vector object
iceland
## Simple feature collection with 1 feature and 94 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -24.53991 ymin: 63.39671 xmax: -13.50292 ymax: 66.56415
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##          featurecla scalerank labelrank sovereignt sov_a3 adm0_dif level
## 188 Admin-0 country         0         3    Iceland    ISL        0     2
##                  type   admin adm0_a3 geou_dif geounit gu_a3 su_dif
## 188 Sovereign country Iceland     ISL        0 Iceland   ISL      0
##     subunit su_a3 brk_diff    name name_long brk_a3 brk_name brk_group
## 188 Iceland   ISL        0 Iceland   Iceland    ISL  Iceland      <NA>
##      abbrev postal           formal_en formal_fr name_ciawf note_adm0
## 188 Iceland     IS Republic of Iceland      <NA>    Iceland      <NA>
##     note_brk name_sort name_alt mapcolor7 mapcolor8 mapcolor9 mapcolor13
## 188     <NA>   Iceland     <NA>         1         4         4          9
##     pop_est pop_rank gdp_md_est pop_year lastcensus gdp_year
## 188  339747       10      16150     2017         NA     2016
##                        economy           income_grp wikipedia fips_10_
## 188 2. Developed region: nonG7 1. High income: OECD        NA       IC
##     iso_a2 iso_a3 iso_a3_eh iso_n3 un_a3 wb_a2 wb_a3   woe_id woe_id_eh
## 188     IS    ISL       ISL    352   352    IS   ISL 23424845  23424845
##                       woe_note adm0_a3_is adm0_a3_us adm0_a3_un adm0_a3_wb
## 188 Exact WOE match as country        ISL        ISL         NA         NA
##     continent region_un       subregion             region_wb name_len
## 188    Europe    Europe Northern Europe Europe & Central Asia        7
##     long_len abbrev_len tiny homepart min_zoom min_label max_label
## 188        7          7   NA        1        0         2         7
##          ne_id wikidataid name_ar name_bn name_de name_en  name_es name_fr
## 188 1159320917       Q189    <NA>    <NA>  Island Iceland Islandia Islande
##     name_el name_hi name_hu  name_id name_it name_ja name_ko name_nl
## 188    <NA>    <NA>  Izland Islandia Islanda    <NA>    <NA> IJsland
##      name_pl  name_pt name_ru name_sv name_tr name_vi name_zh
## 188 Islandia Islândia    <NA>  Island Izlanda Iceland    <NA>
##                           geometry
## 188 MULTIPOLYGON (((-14.56363 6...
#here Iceland
plot(iceland)

By default, the plot( ) function with the class sf creates as many facets of the map as there are variables in it. To limit this behavior we can use either a variable name plot(iceland["admin"]) or the limit argument plot(iceland, max.plot = 1). With the argument max.plot = 1 the function uses the first available variable of the map.

In addition, we see in the information of the object sf that the projection is WGS84 with decimal degrees (EPSG code: 4326). For the calculation of distances it is more convenient to use meters instead of degrees. Because of this, the first thing we do is to transform the map of Iceland to UTM Zone 27 (EPSG code: 3055). More information about EPSG and projections here. For that purpose, we use the st_transform( ) function. We simply indicate the map and the EPSG code.

#transform to UTM
iceland <- st_transform(iceland,3055)

Create a fishnet of points

We still need the points where we want to know the distance. In our case it will be a regular fishnet of points in Iceland with a resolution of 5km. We do this with the function st_make_grid( ), indicating the resolution in the unit of the coordinate system (meters in our case) with the argument cellsize, and what geometry we would like to create what (polygons, centers or corners).

#create the fishnet
grid <- st_make_grid(iceland,cellsize = 5000,what = "centers")

#our fishnet with the extension of Iceland
plot(grid)

#only extract the points in the limits of Iceland
grid <- st_intersection(grid,iceland)   

#our fishnet now
plot(grid)

Calculating the distance

To estimate the distance we use the st_distance( ) function that returns a vector of distances for all our points in the fishnet. But first it is necessary to transform the map of Iceland from a polygon shape (MULTIPOLYGON) to a line (MULTILINESTRING). More details with ?st_cast.

#transform Iceland from polygon shape to line
iceland <- st_cast(iceland,"MULTILINESTRING")

#calculation of the distance between the coast and our points
dist <- st_distance(iceland,grid)

#distance with unit in meters
head(dist)
## Units: [m]
## [1]  790.7906 1151.4360 1270.7603 3128.9057 2428.5677 4197.7472

Plotting the calculated distance

Once obtained the distance for our points, we can combine them with the coordinates and plot them in ggplot2. For this, we create a data.frame. The object dist is a matrix of one column, so we have to convert it to a vector with the function as.vector( ). In addition, we divide by 1000 to convert the distance in meters to km. The st_coordinates( ) function extracts the coordinates of our points. For the final visualization we use a vector of colors with the RdGy palette (more here).

#create a data.frame with the distance and the coordinates of the points
df <- data.frame(dist=as.vector(dist)/1000,
                    st_coordinates(grid))

#structure
str(df)
## 'data.frame':    4104 obs. of  3 variables:
##  $ dist: num  0.791 1.151 1.271 3.129 2.429 ...
##  $ X   : num  608796 613796 583796 588796 593796 ...
##  $ Y   : num  7033371 7033371 7038371 7038371 7038371 ...
#colors 
col_dist <- brewer.pal(11,"RdGy")


ggplot(df,aes(X,Y,fill=dist))+ #variables
         geom_tile()+ #geometry
           scale_fill_gradientn(colours=rev(col_dist))+ #colors for plotting the distance
             labs(fill="Distance (km)")+ #legend name
             theme_void()+ #map theme
              theme(legend.position = "bottom") #legend position

Export the distance as a raster

To be able to export the estimated distance to the sea of Iceland, we need to use the rasterize( ) function of the library raster.

1. First, it is necessary to create an empty raster. In this raster we have to indicate the resolution, in our case it is of 5000m, the projection and the extension of the raster.

1. We can extract the projection from the information of the map of Iceland.

2. The extension can be extracted from our grid points with the function extent( ). However, this last function needs the class sp, so we pass the object grid in sf format, only for this time, to the class sp using the function as( ) and the argument “Spatial”.

2. In addition to the above, the data.frame df, that we created earlier, has to be converted into the sf class. Therefore, we apply the function st_as_sf( ) with the argument coords indicating the names of the coordinates. Additionally, we also define the coordinate system that we know.

#get the extension
ext <- extent(as(grid,"Spatial"))

#extent object
ext
## class       : Extent 
## xmin        : 338795.6 
## xmax        : 848795.6 
## ymin        : 7033371 
## ymax        : 7383371
#raster destination
r <- raster(resolution=5000,ext=ext,crs="+proj=utm +zone=27 +ellps=intl +towgs84=-73,47,-83,0,0,0,0 +units=m +no_defs")

#convert the points to a spatial object class sf
dist_sf <- st_as_sf(df,coords=c("X","Y"))%>%
                      st_set_crs(3055)

#create the distance raster
dist_raster <- rasterize(dist_sf,r,"dist",fun=mean)

#raster
dist_raster
## class       : RasterLayer 
## dimensions  : 70, 102, 7140  (nrow, ncol, ncell)
## resolution  : 5000, 5000  (x, y)
## extent      : 338795.6, 848795.6, 7033371, 7383371  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=27 +ellps=intl +towgs84=-73,47,-83,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : 0.006124901, 115.1712  (min, max)
#plot the raster
plot(dist_raster)

#export the raster
writeRaster(dist_raster,file="dist_islandia.tif",format="GTiff", overwrite=TRUE)

The rasterize( ) function is designed to create rasters from an irregular grid. In case we have a regular grid, like this one, we can use an easier alternative way. The rasterFromXYZ( ) function converts a data.frame with longitude, latitude and the variable Z into a raster object. It is important that the order should be longitude, latitude, variables.

r <- rasterFromXYZ(df[,c(2:3,1)],crs="+proj=utm +zone=27 +ellps=intl +towgs84=-73,47,-83,0,0,0,0 +units=m +no_defs")

plot(r)

With the calculation of distance we can create art, as seen in the header of this post, which includes a world map only with the distance to the sea of all continents. A different perspective to our world (here more (spanish)) .