Calculating the distance to the sea in R
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The distance to the sea is a fundamental variable in geography, especially relevant when it comes to modeling. For example, in interpolations of air temperature, the distance to the sea is usually used as a predictor variable, since there is a casual relationship between the two that explains the spatial variation. How can we estimate the (shortest) distance to the coast in R?
Packages
In this post we will use the following libraries:
Library | Description |
---|---|
tidyverse | Collection of packages (visualization, manipulation): ggplot2, dplyr, etc. |
sf | Simple Feature: import, export and manipulate vector data |
raster | Import, export and manipulate raster |
rnaturalearth | Set of vector maps ‘natural earth’ |
RColorBrewer | Color palettes |
#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.
- 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.
We can extract the projection from the information of the map of Iceland.
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 functionas( )
and the argument “Spatial”.
- 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)) .
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.