How to Read a Shapefile in R

Open and Plot Shapefiles in R

Overview

Teaching: 10 min
Exercises: 0 min

Questions

• Getting started with betoken, line and polygon vector data.

Objectives

• Know the difference between point, line, and polygon vector elements.

• Empathize the differences betwixt opening bespeak, line and polygon shapefiles in R.

• Understand the components of a spatial object in R.

In this tutorial, we will open and plot signal, line and polygon vector information stored in shapefile format in R.

About Vector Data

Vector data are composed of discrete geometric locations (ten, y values) known as vertices that define the "shape" of the spatial object. The organization of the vertices, determines the type of vector that we are working with: point, line or polygon.

There are 3 types of vector objects: points, lines or polygons. Each object type has a different construction. Image Source: National Ecological Observatory Network (NEON)

• Points: Each individual point is defined by a single x, y coordinate. There tin exist many points in a vector bespeak file. Examples of signal information include: sampling locations, the location of individual trees or the location of plots.
• Lines: Lines are composed of many (at least two) vertices, or points, that are continued. For instance, a road or a stream may be represented by a line. This line is composed of a series of segments, each "bend" in the route or stream represents a vertex that has defined 10, y location.
• Polygons: A polygon consists of three or more vertices that are connected and "closed". Thus the outlines of plot boundaries, lakes, oceans, and states or countries are often represented by polygons. Occasionally, a polygon can have a hole in the center of it (like a doughnut), this is something to be aware of but not an event we will bargain with in this tutorial.

Data Tip

Sometimes, boundary layers such as states and countries, are stored as lines rather than polygons. Withal, these boundaries, when represented every bit a line, will not create a closed object with a defined "area" that can be "filled".

Shapefiles: Points, Lines, and Polygons

Geospatial data in vector format are often stored in a shapefile format. Because the structure of points, lines, and polygons are different, each individual shapefile can only contain ane vector type (all points, all lines or all polygons). Yous will not find a mixture of point, line and polygon objects in a single shapefile.

Objects stored in a shapefile oftentimes take a prepare of associated attributes that describe the information. For example, a line shapefile that contains the locations of streams, might contain the associated stream name, stream "order" and other information about each stream line object.

• More than virtually shapefiles can found on Wikipedia.

Import Shapefiles

We volition apply the sf package to work with vector data in R. Notice that the rgdal package automatically loads when sf is loaded. We will also load the raster package and then nosotros tin explore raster and vector spatial metadata using similar commands.

            # load required libraries # for vector work library(sf)                      

            Linking to GEOS 3.6.1, GDAL 2.i.three, proj.4 4.9.3                      

            # for metadata/attributes- vectors or rasters library(raster)                      

            Loading required package: sp                      

The shapefiles that nosotros will import are:

• A polygon shapefile representing our field site purlieus,
• A line shapefile representing roads, and
• A point shapefile representing the location of the Fisher flux belfry located at the NEON Harvard Forest field site.

The first shapefile that nosotros will open contains the purlieus of our study area (or our Expanse Of Involvement or AOI, hence the name aoiBoundary). To import shapefiles we employ the sf function st_read.

st_read() requires the file path to the shapefile.

Allow's import our AOI.

            # Import a polygon shapefile aoi_boundary_HARV <- st_read(   "data/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")                      

            Reading layer `HarClip_UTMZ18' from data source `/Users/ebecker/Box Sync/Carpentry_repos/datacarpentry-lessons/geospatial/r-raster-vector-geospatial/_episodes_rmd/information/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp' using driver `ESRI Shapefile' Simple feature collection with 1 feature and 1 field geometry type:  POLYGON dimension:      XY bbox:           xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359 epsg (SRID):    32618 proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs                      

When we import the HarClip_UTMZ18 shapefile layer into R (as our aoi_boundary_HARV object), the st_read() function automatically stores information near the information. Nosotros are especially interested in the geospatial metadata, describing the format, CRS, extent, and other components of the vector information, and the attributes which describe properties associated with each individual vector object.

Data Tip

The Shapefile Metadata & Attributes in R tutorial provides more data on both metadata and attributes and using attributes to subset and plot data.

Key metadata for all shapefiles include:

1. Object Type: the class of the imported object.
2. Coordinate Reference System (CRS): the projection of the data.
3. Extent: the spatial extent (geographic area that the shapefile covers) of the shapefile. Note that the spatial extent for a shapefile represents the extent for ALL spatial objects in the shapefile.

We can view shapefile metadata using the st_geometry_type, st_crs and st_bbox methods:

            # view just the geometry type for the shapefile st_geometry_type(aoi_boundary_HARV)                      

            [ane] POLYGON 18 Levels: GEOMETRY Indicate LINESTRING POLYGON ... TRIANGLE                      

            # view just the crs for the shapefile st_crs(aoi_boundary_HARV)                      

            Coordinate Reference System:   EPSG: 32618    proj4string: "+proj=utm +zone=eighteen +datum=WGS84 +units=thousand +no_defs"                      

            # view just the extent for the shapefile st_bbox(aoi_boundary_HARV)                      

                          xmin      ymin      xmax      ymax   732128.0 4713208.7  732251.1 4713359.2                      

            # view all metadata at aforementioned time aoi_boundary_HARV                      

            Simple feature drove with ane feature and 1 field geometry type:  POLYGON dimension:      XY bbox:           xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359 epsg (SRID):    32618 proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=chiliad +no_defs   id                       geometry ane  1 POLYGON ((732128 4713359, 7...                      

Our aoi_boundary_HARV is an sf polygon object, in the CRS UTM zone 18N. The CRS is critical to interpreting the object extent values as information technology specifies units.

The spatial extent of a shapefile or R spatial object represents the geographic "edge" or location that is the furthest north, south east and west. Thus is represents the overall geographic coverage of the spatial object. Epitome Source: National Ecological Observatory Network (NEON)

Spatial Data Attributes

Each object in a shapefile has one or more than attributes associated with information technology. Shapefile attributes are similar to fields or columns in a spreadsheet. Each row in the spreadsheet has a set of columns associated with it that draw the row element. In the example of a shapefile, each row represents a spatial object - for example, a road, represented as a line in a line shapefile, will have one "row" of attributes associated with it. These attributes can include different types of information that describe objects stored within a shapefile. Thus, our route, may take a name, length, number of lanes, speed limit, type of road and other attributes stored with it.

Each spatial feature in an R spatial object has the same set of associated attributes that describe or characterize the feature. Attribute data are stored in a split up *.dbf file. Attribute information can be compared to a spreadsheet. Each row in a spreadsheet represents one feature in the spatial object. Epitome Source: National Ecological Observatory Network (NEON)

Nosotros can view the attributes of an sf object past printing it to the screen. The geometry of the object can exist dropped by turning the object into a data.frame.

            Simple feature drove with 1 characteristic and 1 field geometry blazon:  POLYGON dimension:      XY bbox:           xmin: 732128 ymin: 4713209 xmax: 732251.ane ymax: 4713359 epsg (SRID):    32618 proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=1000 +no_defs   id                       geometry one  one POLYGON ((732128 4713359, 7...                      

            # alternate mode to view attributes data.frame(aoi_boundary_HARV)                      

                          id                       geometry 1  1 POLYGON ((732128 4713359, seven...                      

In this example, our polygon object but has one attribute: id.

We tin can view a metadata & aspect summary of each shapefile by inbound the name of the R object in the console. Note that the metadata output includes the geometry type, the number of features, the extent, and the coordinate reference organisation (crs) of the R object. The output of summary() shows a preview of the R object attributes.

            # view a summary of metadata & attributes associated with the spatial object summary(aoi_boundary_HARV)                      

                          id             geometry  Min.   :1   POLYGON      :i    1st Qu.:ane   epsg:32618   :0    Median :ane   +proj=utm ...:0    Mean   :1                      third Qu.:i                      Max.   :1                      

Plot a Shapefile

Next, let's visualize the data in our sf object using plot().

            # create a plot of the shapefile # 'lwd' sets the line width # 'col' sets internal color # 'border' sets line colour plot(aoi_boundary_HARV, col = "cyan1", border = "black", lwd = 3,      primary = "AOI Purlieus Plot")                      

Challenge: Import Line and Point Shapefiles

Using the steps in a higher place, import the HARV_roads and HARVtower_UTM18N layers into R. Phone call the Harv_roads object lines_HARV and the HARVtower_UTM18N point_HARV.

Respond the following questions:

1. What type of R spatial object is created when you import each layer?

2. What is the CRS and extentfor each object?

3. Do the files contain, points, lines or polygons?

4. How many spatial objects are in each file?

Answers

                # import line shapefile lines_HARV <- st_read("information/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")                              

                Reading layer `HARV_roads' from information source `/Users/ebecker/Box Sync/Carpentry_repos/datacarpentry-lessons/geospatial/r-raster-vector-geospatial/_episodes_rmd/information/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp' using commuter `ESRI Shapefile' Uncomplicated characteristic collection with 13 features and 15 fields geometry type:  MULTILINESTRING dimension:      XY bbox:           xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260 epsg (SRID):    32618 proj4string:    +proj=utm +zone=xviii +datum=WGS84 +units=m +no_defs                              

                # import point shapefile point_HARV <- st_read("information/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")                              

                Reading layer `HARVtower_UTM18N' from information source `/Users/ebecker/Box Sync/Carpentry_repos/datacarpentry-lessons/geospatial/r-raster-vector-geospatial/_episodes_rmd/information/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp' using driver `ESRI Shapefile' Simple feature drove with 1 feature and 14 fields geometry type:  POINT dimension:      XY bbox:           xmin: 732183.2 ymin: 4713265 xmax: 732183.2 ymax: 4713265 epsg (SRID):    32618 proj4string:    +proj=utm +zone=eighteen +datum=WGS84 +units=g +no_defs                              

                Coordinate Reference System:   EPSG: 32618    proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"                              

                                  xmin      ymin      xmax      ymax   730741.2 4711942.0  733295.5 4714260.0                              

                Coordinate Reference Organization:   EPSG: 32618    proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"                              

                                  xmin      ymin      xmax      ymax   732183.2 4713265.0  732183.2 4713265.0                              

                # iii #lines_HARV contains only lines and point_HARV contains only 1 point  # four -> numerous ways to find this; lines_HARV=thirteen, nrow(lines_HARV)  #easiest, merely not previously taught                              

                lines_HARV  #expect at 'features'                              

                Elementary feature drove with 13 features and 15 fields geometry type:  MULTILINESTRING dimension:      XY bbox:           xmin: 730741.ii ymin: 4711942 xmax: 733295.5 ymax: 4714260 epsg (SRID):    32618 proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=g +no_defs First 10 features:    OBJECTID_1 OBJECTID       Type             NOTES MISCNOTES RULEID one          14       48 woods road Locust Opening Rd      <NA>      v ii          40       91   footpath              <NA>      <NA>      six three          41      106   footpath              <NA>      <NA>      6 4         211      279 rock wall              <NA>      <NA>      i five         212      280 stone wall              <NA>      <NA>      1 6         213      281 rock wall              <NA>      <NA>      1 7         214      282 stone wall              <NA>      <NA>      1 8         215      283 rock wall              <NA>      <NA>      1 9         216      284 stone wall              <NA>      <NA>      ane ten        553      674  boardwalk              <NA>      <NA>      two             MAPLABEL SHAPE_LENG             LABEL BIKEHORSE RESVEHICLE 1  Locust Opening Rd 1297.35706 Locust Opening Rd         Y         R1 2               <NA>  146.29984              <NA>         Y         R1 3               <NA>  676.71804              <NA>         Y         R2 4               <NA>  231.78957              <NA>      <NA>       <NA> 5               <NA>   45.50864              <NA>      <NA>       <NA> 6               <NA>  198.39043              <NA>      <NA>       <NA> 7               <NA>  143.19240              <NA>      <NA>       <NA> 8               <NA>   90.33118              <NA>      <NA>       <NA> ix               <NA>   35.88146              <NA>      <NA>       <NA> 10              <NA>   67.43464              <NA>         North         R3    RECMAP Shape_Le_1                            ResVehic_1 i       Y 1297.10617    R1 - All Research Vehicles Immune 2       Y  146.29983    R1 - All Inquiry Vehicles Allowed 3       Y  676.71807 R2 - 4WD/High Clearance Vehicles Only 4    <NA>  231.78962                                  <NA> 5    <NA>   45.50859                                  <NA> 6    <NA>  198.39041                                  <NA> 7    <NA>  143.19241                                  <NA> 8    <NA>   90.33114                                  <NA> 9    <NA>   35.88152                                  <NA> 10      N   67.43466              R3 - No Vehicles Allowed                     BicyclesHo                       geometry ane  Bicycles and Horses Allowed MULTILINESTRING ((730819.two ... ii  Bicycles and Horses Allowed MULTILINESTRING ((732040.2 ... iii  Bicycles and Horses Allowed MULTILINESTRING ((732057 47... 4                         <NA> MULTILINESTRING ((731903.6 ... v                         <NA> MULTILINESTRING ((732039.1 ... 6                         <NA> MULTILINESTRING ((732056.2 ... 7                         <NA> MULTILINESTRING ((731964 47... 8                         <NA> MULTILINESTRING ((732105.2 ... 9                         <NA> MULTILINESTRING ((732222.nine ... 10      DO Not Evidence ON REC MAP MULTILINESTRING ((732153.viii ...                              

Plot Multiple Shapefiles

The plot() function can exist used for bones plotting of spatial objects. We utilise the add = TRUE argument to overlay shapefiles on top of each other, equally we would when creating a map in a typical GIS application like QGIS.

We can use main = "" to give our plot a title. If nosotros desire the championship to span two lines, we use \n where the line should intermission.

            # Plot multiple shapefiles plot(aoi_boundary_HARV, col = "lightgreen",      principal = "NEON Harvard Forest\nField Site") plot(lines_HARV, add = TRUE)                      

            Alert in plot.sf(lines_HARV, add = TRUE): ignoring all but the showtime attribute                      

            # use the pch element to adjust the symbology of the points plot(point_HARV, add  = Truthful, pch = 19, col = "purple")                      

            Warning in plot.sf(point_HARV, add = TRUE, pch = xix, col = "purple"): ignoring all only the offset attribute                      

Data Tip

The pch argument specifies the point shape. A list of valid signal shapes can be found by viewing this graphic

Challenge: Plot Raster & Vector Data Together

You lot can plot vector data layered on elevation of raster data using the add together = TRUE plot attribute. Create a plot that uses the NEON AOP Canopy Height Model NEON_RemoteSensing/HARV/CHM/HARV_chmCrop.tif every bit a base of operations layer. On tiptop of the CHM, please add:

• The study site AOI.
• Roads.
• The tower location.

Be sure to give your plot a meaningful championship.

For assistance consider using the Shapefile Metadata & Attributes in R, the Plot Raster Data in R tutorials.

Answers

                # import CHM chm_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")  plot(chm_HARV,      main = "Map of Study Area\north w/ Canopy Superlative Model\nNEON Harvard Forest Field Site")  plot(lines_HARV,      add together = TRUE,      col = "blackness")                              

                Warning in plot.sf(lines_HARV, add = Truthful, col = "blackness"): ignoring all but the first attribute                              

                plot(aoi_boundary_HARV, border="grey20",      add = TRUE,      lwd = 4) plot(point_HARV, pch=viii,      add = True)                              

                Warning in plot.sf(point_HARV, pch = eight, add = True): ignoring all just the showtime aspect                              

Boosted Resources: Plot Parameter Options

For more on parameter options in the base R plot() function, check out these resources:

• Parameter methods in R.
• Color names in R

Key Points

viningtherear.blogspot.com

Source: https://erinbecker.github.io/r-raster-vector-geospatial/06-vector-open-shapefile-in-r/index.html

Share this post