Shapefiles

In many instances data is available only in the form of shapefiles. A shapefile is a simple format for storing geometric location and attribute information of geographic features. The shapefile is a popular geospatial vector format developed and regulated by ESRI as a (mostly) open specification. Geographic features in a shapefile can be represented by points, lines, or polygons (areas).

“R” contains various functions to read and write shapefiles. Here I show how to convert the new 2462 ft contour obtained from Century Engineering for the shoreline of Deep Creek Lake from a standard text format to a shapefile. This data set was generated as part of Century Engineering’s surevying work done for the lake buydown program.

Analysis

The following are the first few lines from a text file of the 2462 ft ASL contour values of Deep Creek Lake, MD, Eastings and Northings in ft.

Easting	Northing
649177.015	698016.962
649182.756	698030.201
649204.067	698060.465
649222.207	698077.754
...

A shapefile needs a CRS specification (Coordinate Reference System). The CRS specification was extracted from an existing shapefile, one provided by the County in 2012 .The data is subsequently plotted.

Here is the code snippet for reading a shapefile and plotting its contents:

shp1file <- "../data/dcl_outline_carpenter/DCLShoreline.shp"
shp1 <- readOGR(shp1file, "DCLShoreline")
plot(shp1, col="red", main="Deep Creek Lake\nCounty Data", lwd=1)

The plot of the data is shown in Figure 1.

The CRS string extracted from this file looks like this:

+proj=lcc +lat_1=38.3 +lat_2=39.45 +lat_0=37.66666666666666 +lon_0=-77

Next is to read the century file and convert it to a shapefile. This data is shown in Figure 2.

After having read in the text file with the Northings and Easting data, the new shapefile is written by the following code snippet:

# To use a shapefile function we must first transform regular text data into
# a spatial points dataframe
# Note that the easting and northing columns are in columns 1 and 2
plot.df1 <- SpatialPointsDataFrame(df1[,1:2],  # The columns to use
                    df1,                       # The R data object to convert
                    proj4string = prjstr)      # Assign a CRS (extracted above)

# Plot the data to check if OK (blue color)
plot(plot.df1, col="blue", cex=0.1)

And here we’re checking to see if the shapefile that we created can also be read and give the right results

# Double check to see of the CRS has been transferred properly
crs(plot.df1)

# First delete the old result
# REF: https://statisticsglobe.com/remove-directory-using-r
shp1file <- "../results/dcl_shoreline_2462"
unlink(shp1file, recursive = TRUE)   # Apply unlink()

# Export the shapefile
writeOGR(plot.df1, shp1file, "century2462-2017", driver="ESRI Shapefile")


# Verify that the shapefile is correctly specified be reading it in.
shp2file <- "../results/dcl_shoreline_2462"

shp2 <- readOGR(shp2file, "century2462-2017")

# An plot to verify (green color)
plot(shp2, col="green", main="Deep Creek Lake - 2017\nVerified Shapefile of Century Data", cex=0.1, lwd=1)

To check that all is done well, we read in the shapefile and plot it. The result is shown in Figure 3.

To conclude: “The conversion back and forth worked successfully.”