[This article was first published on r-spatial, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.* Introduction * Vector loading + Shapefile (.shp) + GeoPackage (.gpkg) + GeoJSON (.geojson) * Vector saving * Synopsis * Supplement + Zipped shapefile (.shz) + Virtual File Systems + SQL preselection - Columns filtering - Rows filtering - Spatial filtering
[view rawRmd]
Summary:
If you have spatial vector data and are wondering how to load / save itin R, this tutorial is the answer to your questions. It presentspractical examples for the most popular formats using thesf package. We will use freevector layers from Natural Earth asa data source.
IntroductionFor convenience, all necessary files are located in the GitHubrepository:
We can download the mentioned data and interactive notebook (.Rmd)manually from the repository (“Code” button > “Download ZIP”) or usethe following script.
url = "https://github.com/kadyb/sf\_load\_save/archive/refs/heads/main.zip"download.file(url, "sf\_load\_save.zip")unzip("sf\_load\_save.zip")
In the first step, we need to download the sf package using theinstall.packages() function, and then use the library() function toload it into the session.
install.packages("sf")library("sf")
Vector loadingShapefile (.shp)Let’s start by loading the shapefile format, which actually consists ofseveral files (e.g., .shp, .shx, .dbf, .prj). More information can befound on Wikipedia, butcurrently it is not recommended to use this format due to its manylimitations.
Generally, we can use the read_sf() function to load data. It requiresproviding a path to the file. The file path can be defined in two waysin R and this is the most common source of problems (errors like:Error: Cannot open "file.shp"; The file doesn't seem to exist.).
The first way (easier) is to provide an absolute path, i.e. we mustprovide the exact location where the file is located. For instance:
path = "C:/Users/Krzysztof/Documents/file.shp"
However, this isnotthe recommended method, as it makes it impossible to locate files ondifferent operating systems. The second way is to specify a relativepath. In this case, we specify the location of the file relative tothe current working directory (or project). To find out where theworking directory is, we can use the getwd() function, and to changeit the setwd() function. For instance:
getwd()#> "C:/Users/Krzysztof/Documents"path = "file.shp"
Let’s load the shapefile using a relative path (all data can be found inthe data folder).
countries = read\_sf("data/countries/countries.shp")
We can then print the metadata about this vector layer by referring tothe countries object.
countries## Simple feature collection with 52 features and 168 fields## Geometry type: MULTIPOLYGON## Dimension: XY## Bounding box: xmin: -17.53604 ymin: -34.82195 xmax: 51.41704 ymax: 37.3452## Geodetic CRS: WGS 84## # A tibble: 52 × 169## featurecla scalerank LABELRANK SOVEREIGNT SOV\_A3 ADM0\_DIF LEVEL TYPE TLC ## <chr> <int> <int> <chr> <chr> <int> <int> <chr> <chr>## 1 Admin-0 cou… 0 2 Ethiopia ETH 0 2 Sove… 1 ## 2 Admin-0 cou… 0 3 South Sud… SDS 0 2 Sove… 1 ## 3 Admin-0 cou… 0 6 Somalia SOM 0 2 Sove… 1 ## 4 Admin-0 cou… 0 2 Kenya KEN 0 2 Sove… 1 ## 5 Admin-0 cou… 0 6 Malawi MWI 0 2 Sove… 1 ## 6 Admin-0 cou… 0 3 United Re… TZA 0 2 Sove… 1 ## 7 Admin-0 cou… 0 5 Somaliland SOL 0 2 Sove… 1 ## 8 Admin-0 cou… 0 3 Morocco MAR 0 2 Sove… 1 ## 9 Admin-0 cou… 0 7 Western S… SAH 0 2 Inde… 1 ## 10 Admin-0 cou… 0 4 Republic … COG 0 2 Sove… 1 ## # ℹ 42 more rows## # ℹ 160 more variables: ADMIN <chr>, ADM0\_A3 <chr>, GEOU\_DIF <int>,## # GEOUNIT <chr>, GU\_A3 <chr>, SU\_DIF <int>, SUBUNIT <chr>, SU\_A3 <chr>,## # BRK\_DIFF <int>, NAME <chr>, NAME\_LONG <chr>, BRK\_A3 <chr>, BRK\_NAME <chr>,## # BRK\_GROUP <chr>, ABBREV <chr>, POSTAL <chr>, FORMAL\_EN <chr>,## # FORMAL\_FR <chr>, NAME\_CIAWF <chr>, NOTE\_ADM0 <chr>, NOTE\_BRK <chr>,## # NAME\_SORT <chr>, NAME\_ALT <chr>, MAPCOLOR7 <int>, MAPCOLOR8 <int>, …
We can see that this layer consists of 52 features (rows) and 168 fields(columns). The next information is about geometry type, dimension,spatial extent (bounding box) and coordinate reference system (CRS). Inaddition, the first 10 rows were printed.
After loading the data, it is a good idea to present it on a map. Asimple plot() function can be used for this purpose. The countriesobject has many fields (attributes), but to start with we only needgeometry. It can be obtained by using the st_geometry() function.
plot(st\_geometry(countries))
GeoPackage (.gpkg)The next dataset is rivers (linear geometry) saved in GeoPackageformat. It is loaded in exactly the sameway as the shapefile before. Note that this format can consist ofmultiple layers of different types. In this case, we must define whichlayer exactly we want to load. To check what layers are in thegeopackage, use the st_layers() function, and then specify it usingthe layer argument in read_sf(). If the file only contains onelayer, we don’t need to do this.
st\_layers("data/rivers.gpkg")## Driver: GPKG ## Available layers:## layer\_name geometry\_type features fields crs\_name## 1 rivers Multi Line String 228 38 WGS 84rivers = read\_sf("data/rivers.gpkg", layer = "rivers")
We can also display metadata as in the previous example.
rivers## Simple feature collection with 228 features and 38 fields## Geometry type: MULTILINESTRING## Dimension: XY## Bounding box: xmin: -16.54233 ymin: -34.34378 xmax: 49.46094 ymax: 35.12311## Geodetic CRS: WGS 84## # A tibble: 228 × 39## dissolve scalerank featurecla name name\_alt rivernum note min\_zoom name\_en## <chr> <int> <chr> <chr> <chr> <int> <chr> <dbl> <chr> ## 1 975River 9 River <NA> <NA> 975 <NA> 7.1 <NA> ## 2 976River 9 River Rung… <NA> 976 <NA> 7.1 Rungwa ## 3 977River 9 River Ligo… <NA> 977 <NA> 7.1 Ligonha## 4 978River 9 River Dong… <NA> 978 <NA> 7.1 Dongwe ## 5 979River 9 River Cuito <NA> 979 <NA> 7.1 Cuito ## 6 980Lake … 9 Lake Cent… <NA> <NA> 980 <NA> 7.1 <NA> ## 7 980River 9 River <NA> <NA> 980 <NA> 7.1 <NA> ## 8 981River 9 River Bagoé <NA> 981 <NA> 7.1 Bagoé ## 9 982River 9 River Hade… <NA> 982 <NA> 7.1 Hadejia## 10 983River 9 River Sous <NA> 983 <NA> 7.1 Sous ## # ℹ 218 more rows## # ℹ 30 more variables: min\_label <dbl>, ne\_id <dbl>, label <chr>,## # wikidataid <chr>, name\_ar <chr>, name\_bn <chr>, name\_de <chr>,## # name\_es <chr>, name\_fr <chr>, name\_el <chr>, name\_hi <chr>, name\_hu <chr>,## # name\_id <chr>, name\_it <chr>, name\_ja <chr>, name\_ko <chr>, name\_nl <chr>,## # name\_pl <chr>, name\_pt <chr>, name\_ru <chr>, name\_sv <chr>, name\_tr <chr>,## # name\_vi <chr>, name\_zh <chr>, name\_fa <chr>, name\_he <chr>, …
And make a visualization, but this time we will plot rivers against thebackground of country borders. Adding more layers to the visualizationis done with the add = TRUE argument in plot() function. Note thatthe order in which objects are added is important – the objects addedlast are displayed at the top. The col argument is used to set thecolor of the object.
plot(st\_geometry(countries))plot(st\_geometry(rivers), add = TRUE, col = "blue")
GeoJSON (.geojson)The last GeoJSON file contains cities in the world. In this case, wealso use the read_sf() function to load this file.
cities = read\_sf("data/cities.geojson")cities## Simple feature collection with 1287 features and 31 fields## Geometry type: POINT## Dimension: XY## Bounding box: xmin: -17.47508 ymin: -34.52953 xmax: 51.12333 ymax: 37.29042## Geodetic CRS: WGS 84## # A tibble: 1,287 × 32## scalerank natscale labelrank featurecla name namepar namealt nameascii## <int> <int> <int> <chr> <chr> <chr> <chr> <chr> ## 1 10 1 8 Admin-1 capital Bassar <NA> <NA> Bassar ## 2 10 1 8 Admin-1 capital Sotou… <NA> <NA> Sotouboua## 3 10 1 7 Admin-1 capital Meden… <NA> <NA> Medenine ## 4 10 1 7 Admin-1 capital Kebili <NA> <NA> Kebili ## 5 10 1 7 Admin-1 capital Tatao… <NA> <NA> Tataouine## 6 10 1 7 Admin-1 capital L'Ari… <NA> <NA> L'Ariana ## 7 10 1 7 Admin-1 capital Jendo… <NA> <NA> Jendouba ## 8 10 1 7 Admin-1 capital Kasse… <NA> <NA> Kasserine## 9 10 1 7 Admin-1 capital Sdid … <NA> <NA> Sdid Bou…## 10 10 1 7 Admin-1 capital Silia… <NA> <NA> Siliana ## # ℹ 1,277 more rows## # ℹ 24 more variables: adm0cap <int>, capalt <int>, capin <chr>,## # worldcity <int>, megacity <int>, sov0name <chr>, sov\_a3 <chr>,## # adm0name <chr>, adm0\_a3 <chr>, adm1name <chr>, iso\_a2 <chr>, note <chr>,## # latitude <dbl>, longitude <dbl>, pop\_max <int>, pop\_min <int>,## # pop\_other <int>, rank\_max <int>, rank\_min <int>, meganame <chr>,## # ls\_name <chr>, min\_zoom <dbl>, ne\_id <int>, geometry <POINT [°]>
In this dataset, there is the featurecla column that indicates thetype of city. So let’s try to print them and then select only statecapitals.
We can print a column (attribute) in two ways, i.e. by specifying thecolumn name in:
cities["featurecla"]## Simple feature collection with 1287 features and 1 field## Geometry type: POINT## Dimension: XY## Bounding box: xmin: -17.47508 ymin: -34.52953 xmax: 51.12333 ymax: 37.29042## Geodetic CRS: WGS 84## # A tibble: 1,287 × 2## featurecla geometry## <chr> <POINT [°]>## 1 Admin-1 capital (0.7890036 9.261)## 2 Admin-1 capital (0.9849965 8.557002)## 3 Admin-1 capital (10.4167 33.4)## 4 Admin-1 capital (8.971003 33.69)## 5 Admin-1 capital (10.4667 33)## 6 Admin-1 capital (10.2 36.86667)## 7 Admin-1 capital (8.749999 36.5)## 8 Admin-1 capital (8.716698 35.2167)## 9 Admin-1 capital (9.500004 35.0167)## 10 Admin-1 capital (9.383302 36.0833)## # ℹ 1,277 more rows# the `head()` function prints only the first 6 elementshead(cities[["featurecla"]])## [1] "Admin-1 capital" "Admin-1 capital" "Admin-1 capital" "Admin-1 capital"## [5] "Admin-1 capital" "Admin-1 capital"# or alternatively# head(cities$featurecla)
This layer contains 1287 different cities. To find out what types ofcities these are, we can use the table() function, which willsummarize them.
table(cities[["featurecla"]])## ## Admin-0 capital Admin-0 capital alt Admin-1 capital ## 54 6 609 ## Admin-1 region capital Populated place ## 19 599
We are interested in Admin-0 capital and Admin-0 capital alt typesbecause some countries have two capitals. We make selection as followsusing the | (OR) operator:
sel = cities$featurecla == "Admin-0 capital" | cities$featurecla == "Admin-0 capital alt"head(sel)## [1] FALSE FALSE FALSE FALSE FALSE FALSE
As a result of this operation, we got a logical vector with TRUE andFALSE values (if the city is / is not the capital). Now let’s create anew object named capitals, which will contain only capitals.
```
``
In the last step, we prepare the final visualization. We can add a title(mainargument), axes (axesargument) and change the backgroundcolor (bgcargument) of the figure. We can also change the pointsymbol (pchargument), set its size (cexargument) and fill color(bg` argument).
plot(st\_geometry(countries), main = "Africa", axes = TRUE, bgc = "deepskyblue", col = "burlywood")plot(st\_geometry(rivers), add = TRUE, col = "blue")plot(st\_geometry(capitals), add = TRUE, pch = 24, bg = "red", cex = 0.8)
Vector savingSaving vector data is as easy as loading. There is a dedicatedwrite_sf() function for this purpose and it requires two arguments:
For example, let’s save our capital object as a GeoPackage (.gpkg),but as an exercise you can save it in other formats as well (you justneed to change the extension).
write\_sf(capitals, "data/capitals.gpkg")
SynopsisThe sf package allows loading vector data with the read_sf()function and saving it with the write_sf() function in R. A list ofall supported vector formats can be found on the GDALwebsite.
For more information, see:
SupplementIn the previous part of the tutorial, we looked at simple examples ofloading vector data, while in this section we will check out moreadvanced ways.
Zipped shapefile (.shz)As we noted earlier, a shapefile consists of several files, which can becumbersome. Some solution is to use zipped shapefiles, which is de factoan archive. To create such a file, the extension .shz (or .shp.zip) andthe ESRI Shapefile driver are required. Loading is done in a standardway by specifying the path to the “.shz” file.
write\_sf(capitals, "data/capitals.shz", driver = "ESRI Shapefile")
Hooray, only one file on the disk!
Virtual File SystemsGDAL provides some facilities for loading files using some abstractionby Virtual FileSystems. In practice,this means that we can refer directly to the files without firstunpacking or downloading them in R. For example, we can directly openthe shapefile that is in the archive on the website. To do this, we mustuse two prefixes:
/vsicurl/ to download the file/vsizip/ to unpack the archive```
``` SQL preselectionWe can use SQL queries topre-filter features, so only selected objects / attributes will beloaded. This allows us to limit the size of the object in memory andspeed up the operation time. Moreover, we can also make spatialselection, i.e. limit the loading of data only to a selected area.
Columns filteringThe query argument in the read_sf() function is used to pass SQLqueries. Let’s go back to the countries dataset and load only thecolumn with the names of countries (NAME_LONG).
sql = "SELECT NAME\_LONG FROM countries"f = "data/countries/countries.shp"read\_sf(f, query = sql)## Simple feature collection with 52 features and 1 field## Geometry type: MULTIPOLYGON## Dimension: XY## Bounding box: xmin: -17.53604 ymin: -34.82195 xmax: 51.41704 ymax: 37.3452## Geodetic CRS: WGS 84## # A tibble: 52 × 2## NAME\_LONG geometry## <chr> <MULTIPOLYGON [°]>## 1 Ethiopia (((34.0707 9.454592, 34.06689 9.531176, 34.09821 9.679…## 2 South Sudan (((35.92084 4.619332, 35.85654 4.619603, 35.78122 4.61…## 3 Somalia (((46.46696 6.538292, 46.48805 6.558645, 46.50841 6.57…## 4 Kenya (((35.70585 4.619447, 35.70594 4.619962, 35.71152 4.66…## 5 Malawi (((34.96461 -11.57356, 34.65125 -11.57004, 34.61673 -1…## 6 Tanzania (((32.92086 -9.4079, 32.90546 -9.398185, 32.83074 -9.3…## 7 Somaliland (((48.93911 11.24913, 48.93911 11.13674, 48.93911 11.0…## 8 Morocco (((-8.817035 27.66146, -8.818449 27.6594, -8.81292 27.…## 9 Western Sahara (((-8.817035 27.66146, -8.816537 27.66147, -8.752562 2…## 10 Republic of the Congo (((18.62639 3.476869, 18.63455 3.449222, 18.64241 3.32…## # ℹ 42 more rows
Rows filteringWe can also select rows using a condition, e.g. population (POP_EST)greater than 25 million.
sql = "SELECT * FROM countries WHERE POP\_EST > 25000000"f = "data/countries/countries.shp"read\_sf(f, query = sql) # 17 countries## Simple feature collection with 17 features and 168 fields## Geometry type: MULTIPOLYGON## Dimension: XY## Bounding box: xmin: -17.01374 ymin: -34.82195 xmax: 50.50392 ymax: 37.09394## Geodetic CRS: WGS 84## # A tibble: 17 × 169## featurecla scalerank LABELRANK SOVEREIGNT SOV\_A3 ADM0\_DIF LEVEL TYPE TLC ## <chr> <int> <int> <chr> <chr> <int> <int> <chr> <chr>## 1 Admin-0 cou… 0 2 Ethiopia ETH 0 2 Sove… 1 ## 2 Admin-0 cou… 0 2 Kenya KEN 0 2 Sove… 1 ## 3 Admin-0 cou… 0 3 United Re… TZA 0 2 Sove… 1 ## 4 Admin-0 cou… 0 3 Morocco MAR 0 2 Sove… 1 ## 5 Admin-0 cou… 0 2 Democrati… COD 0 2 Sove… 1 ## 6 Admin-0 cou… 0 2 South Afr… ZAF 0 2 Sove… 1 ## 7 Admin-0 cou… 0 3 Sudan SDN 0 2 Sove… 1 ## 8 Admin-0 cou… 0 3 Ivory Coa… CIV 0 2 Sove… 1 ## 9 Admin-0 cou… 0 2 Nigeria NGA 0 2 Sove… 1 ## 10 Admin-0 cou… 0 3 Angola AGO 0 2 Sove… 1 ## 11 Admin-0 cou… 0 3 Algeria DZA 0 2 Sove… 1 ## 12 Admin-0 cou… 0 3 Mozambique MOZ 0 2 Sove… 1 ## 13 Admin-0 cou… 0 3 Uganda UGA 0 2 Sove… 1 ## 14 Admin-0 cou… 0 3 Cameroon CMR 0 2 Sove… 1 ## 15 Admin-0 cou… 0 3 Ghana GHA 0 2 Sove… 1 ## 16 Admin-0 cou… 0 2 Egypt EGY 0 2 Sove… 1 ## 17 Admin-0 cou… 0 3 Madagascar MDG 0 2 Sove… 1 ## # ℹ 160 more variables: ADMIN <chr>, ADM0\_A3 <chr>, GEOU\_DIF <int>,## # GEOUNIT <chr>, GU\_A3 <chr>, SU\_DIF <int>, SUBUNIT <chr>, SU\_A3 <chr>,## # BRK\_DIFF <int>, NAME <chr>, NAME\_LONG <chr>, BRK\_A3 <chr>, BRK\_NAME <chr>,## # BRK\_GROUP <chr>, ABBREV <chr>, POSTAL <chr>, FORMAL\_EN <chr>,## # FORMAL\_FR <chr>, NAME\_CIAWF <chr>, NOTE\_ADM0 <chr>, NOTE\_BRK <chr>,## # NAME\_SORT <chr>, NAME\_ALT <chr>, MAPCOLOR7 <int>, MAPCOLOR8 <int>,## # MAPCOLOR9 <int>, MAPCOLOR13 <int>, POP\_EST <dbl>, POP\_RANK <int>, …
Spatial filteringFinally, to perform spatial filtering, we must first define the spatialextent / bounding box (st_bbox() function) and specify its coordinatereference system (CRS). Then the bounding box needs to be converted intoa polygon using the st_as_sfc() function and finally converted to aWell-KnowTextrepresentation using st_as_text() function. Therefore, prepared textis passed to the wkt_filter argument. Follow the example below ofloading rivers only in southern Africa:
bbox = st\_bbox(c(xmin = 10, xmax = 40, ymax = -35, ymin = -20), crs = st\_crs(4326))bbox = st\_as\_text(st\_as\_sfc(bbox))bbox## [1] "POLYGON ((10 -20, 40 -20, 40 -35, 10 -35, 10 -20))"f = "data/rivers.gpkg"rivers\_south = read\_sf(f, wkt\_filter = bbox)plot(st\_geometry(rivers\_south), axes = TRUE)
To leave a comment for the author, please follow the link and comment on their blog: r-spatial.
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.Continue reading: How to load and save vector data in R