CoveMaps.R

The process of generating the bathymetric maps is performed by three R scripts (R is a computer language):

  1. ExtractRegions.R - This script extracts the data sets specific for a small defined region of the lake.
  2. CoveMaps.R - This script generates the actual map of a cove using the data generated by the script specified under (1)
  3. BatchRun.R - This script runs all of the small defined regions and generates all of the maps (30 as of this writing). This is the only process that is somewhat time consuming on the computer.

This note describes the second of the two R scripts.

Appendix A has the full script that produces the bathymetry map for one of the coves.

There are several things that one should be aware of:

  1. This is a volunteer effort. This did not cost anyone anything. This was all done on my personal time.
  2. This is still a work in progress because there are a number of technical issues that need to be resolved to make the most accurate maps possible. A lot more attention has to be paid to the boundary conditions and other interpolation algorithms for the irregularly spaced measures need to be tried. Furthermore, some older data should be folded into the analysis.
  3. The R language was being learned as the work was being done and hence one will come across ways of doing things with R that are more efficient and more elegant, and perhaps more accurate.

Waiting for a process to finish on the computer was ever an issue. A lot of Google searches were performed and some of the most useful, or groundbreaking solutions, are put as comments in the script.

The maps that are on the website are produced entirely with R. Three data sets were used:

  1. The contour of the lake
  2. The bathymetry measurements
  3. The boat slip locations

Each of these files comprise data for the whole lake. The ExtractRegion.R script extracts only those data that are pertinent to the cove or lake area under investigation. That cove or area is defined by a rectangular bounding box with the coordinates of the box and name of the area encoded in a small file.

An example of the final result is shown in the following figure:

Sample bathymetric image

PLV: 2/8/2013

This is the header part of the CoveMaps.R script.. it basically describes some of the highlights of the development process.

# Title:       Contours using Akima interpolation and Loess smoothing
# Description: Provide a nice looking plot of contours for a specified cove or lake section
# The bounding box is used to define a matrix of points within which the cove is located.
# Date:        10/9/2012
# Project:     Bathymetry of Deep Creek lake
# File:        CoveMapsV10x.R
# Author:      Pete Versteegen
# Revised:     11/15/12; this has been under development for a while; it's final now; some small fixings to work on
# Revised:     11/16/12; new inout data descriptor file
# Revised:     11/21/12; added a graphics descriptor file to specify graphics parameters unique to each cove
#              It's now pretty much final.
#======================================================================================================

R was learned on the fly, most it via Googling for answers to specific issues or problems that need to be solved. The following references were instrumental in helping figure out how to solve a particular problem.

# Reference Material
# http://stackoverflow.com/questions/11666197/create-colour-grid-over-map-in-r
# http://stackoverflow.com/questions/10047870/plotting-interpolated-data-on-map/10058844#10058844
# http://stackoverflow.com/questions/8563334/average-values-of-a-point-dataset-to-a-grid-dataset/8563683#8563683
# http://www.slideshare.net/tjagger/r-spatial-analysis-using-sp
# http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/geoR/html/polygrid.html
# http://stackoverflow.com/questions/10047870/plotting-interpolated-data-on-map
# par( "usr" ) returns a vector with the extreme coordinates of the user’s plotting region
# Ref: http://www.statmethods.net/advgraphs/axes.html
# Ref: http://sphaerula.com/legacy/R/placingTextInPlots.html
# Ref: http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html
# Ref: http://onertipaday.blogspot.com/2008/09/fitting-text-under-plot.html
# Reference for the following:  http://www.image.ucar.edu/GSP/Software/Fields/Help/image.plot.html
# http://stat.ethz.ch/R-manual/R-patched/library/graphics/html/legend.html
#======================================================================================================

This marks the beginning of the run. It's used later on to determine the running time of the problem.

begTime <- Sys.time()

Doing the various data manipulations requires functions defined in a number of passages. The following is a list of packages that were loaded. It's possible that some are redundant or not needed, but the effort has not been put in to clean up thai portion of the code.

#======================================================================================================
# Load the needed libraries - quit the output when loading, supposedly...

suppressMessages(library(akima))
suppressMessages(library(maptools))
suppressMessages(library(PBSmapping))
suppressMessages(library(RColorBrewer))
suppressMessages(library(fields))
suppressMessages(library(sp))
suppressMessages(library(tgp))
suppressMessages(library(splancs))
suppressMessages(library(SDMTools))
suppressMessages(library(spatstat))
suppressMessages(library(raster))
suppressMessages(library(geoR))
suppressMessages(library(gtools))
suppressMessages(library(GISTools))    # North arrow and map scale
suppressMessages(library(lattice))
suppressMessages(library(rgeos))       # To manipulate polygons.
suppressMessages(library(calibrate))   # To label the points of interest

For each or the coves or area of the lake a small file specifies the name of the cove or are and the resolution of the desired map. The name and the resolution are part of the subsetted data produced by the ExtractRegions.R script. Also read in is the supplemental information pertaining to the named cove or area.

#======================================================================================================
# Read the problem specification

con      <- file("case_input.txt", "rt")
project  <- readLines(con,1)     # Root name of the file
resolut  <- readLines(con,2)     # Resolution in feet
close(con)
resolution <- as.numeric(resolut)
project

#======================================================================================================
# Read the supplemental information for the cove

fileName <- paste("./run_files/", project,"_runfile", "_", resolution,".txt", sep="")
fileConn <- file(fileName, "r")  # open an input file connection
project_data <- read.table(fileConn, header=FALSE, sep="\t", comment.char = "#")
project_data
close(fileConn)

projectp       <- trim(as.character(project_data[1,2]))
areaname       <- trim(as.character(project_data[2,2]))
nobathpoints   <- as.numeric(as.character(project_data[3,2]))
resolution     <- as.numeric(as.character(project_data[4,2]))
nshore         <- as.numeric(as.character(project_data[5,2]))
nslips         <- as.numeric(as.character(project_data[6,2]))
xminbox        <- as.numeric(as.character(project_data[7,2]))
xmaxbox        <- as.numeric(as.character(project_data[8,2]))
yminbox        <- as.numeric(as.character(project_data[9,2]))
ymaxbox        <- as.numeric(as.character(project_data[10,2]))
areap          <- as.numeric(as.character(project_data[11,2]))
milesp         <- as.numeric(as.character(project_data[12,2]))
nx             <- as.numeric(as.character(project_data[13,2]))
ny             <- as.numeric(as.character(project_data[14,2]))

The following sections accomplish the following:

  1. Set up the rectangular bounding box
  2. Read the outline of the cove or area
  3. Read in the locations of the boat slips
  4. Read in points of interest
  5. Read the bathymetric data
  6. Read the graphics parameters to locate the various elements on the final image
#======================================================================================================
# Define the bounding box as a polygon

xbox      <- c(xminbox, xmaxbox, xmaxbox, xminbox)
ybox      <- c(yminbox, yminbox, ymaxbox, ymaxbox)
box_cbind <- cbind(xbox, ybox)
deltax    <- (xmaxbox-xminbox)/nx
deltay    <- (ymaxbox-yminbox)/ny


#======================================================================================================
# Read the cove outline

fileName    <- paste("./run_files/", project, "_shoreline.txt", sep="")
fileConn    <- file(fileName, "r")  # open an input file connection
covename2   <- readLines(fileConn, n=1)
coveoutline <- read.table(fileConn, header=TRUE, sep="\t", comment.char = "#")
close(fileConn)

xcove      <- coveoutline$x
ycove      <- coveoutline$y
ncove      <- length(xcove)
cove_cbind <- cbind(xcove, ycove)


#======================================================================================================
# Read the boat slip locations

fileName  <- paste("./run_files/", project,"_boatslips.txt", sep="")
fileConn  <- file(fileName, "r")  # open an input file connection
covename3 <- readLines(fileConn, n=1)
docks     <- read.table(fileConn, header=TRUE, sep=" ", comment.char = "#")
close(fileConn)

xdock   <- docks$x
ydock   <- docks$y
ndocks  <- length(xdock)
dock.cb <- cbind(xdock, ydock)


#======================================================================================================
# Read the locations of points of interest

fileName  <- paste("./run_files/", project,"_poi.txt", sep="")
fileConn  <- file(fileName, "r")  # open an input file connection
covename4 <- readLines(fileConn, n=1)
pois     <- read.table(fileConn, header=TRUE, sep=" ", comment.char = "#")
close(fileConn)

poin   <- as.character(pois$name)
xpoi   <- as.numeric(pois$x)
ypoi   <- as.numeric(pois$y)
npoi   <- length(xpoi)
poi.cb <- cbind(xpoi, ypoi)
pois


#======================================================================================================
# Read the bathymetry

fileName   <- paste("./run_files/", project,"_bathymetry.txt", sep="")
fileConn   <- file(fileName, "r")  # open an input file connection
covename4  <- readLines(fileConn, n=1)
bathymetry <- read.table(fileConn, header=TRUE, sep=" ", comment.char = "#")
close(fileConn)


xbath   <- bathymetry$x
ybath   <- bathymetry$y
zbath   <- bathymetry$z
bath.cb <- cbind(xbath, ybath)


#======================================================================================================
# Read the graphics parameters

fileName <- paste("./graphing_parameters/", project,"_graphics_parameters.txt", sep="")
fileConn <- file(fileName, "r")  # open an input file connection
graph_parms <- read.table(fileConn, header=TRUE, sep="\t", comment.char = "#")
graph_parms
close(fileConn)

deltaxleft      <- as.numeric(as.character(graph_parms[1 ,2]))
deltaytop       <- as.numeric(as.character(graph_parms[2 ,2]))
linespacing     <- as.numeric(as.character(graph_parms[3 ,2]))
fontsize        <- as.numeric(as.character(graph_parms[4 ,2]))
legendlocation  <- trim(as.character(graph_parms[5 ,2]))
legendinset     <- as.numeric(as.character(graph_parms[6 ,2]))
minilakex0      <- as.numeric(as.character(graph_parms[7 ,2]))
minilakey0      <- as.numeric(as.character(graph_parms[8 ,2]))
scalefactor     <- as.numeric(as.character(graph_parms[9 ,2]))

scalearrow      <- as.numeric(as.character(graph_parms[10 ,2]))
xarrowpct       <- as.numeric(as.character(graph_parms[11 ,2]))
yarrowpct       <- as.numeric(as.character(graph_parms[12 ,2]))

scaleline       <- as.numeric(as.character(graph_parms[13 ,2]))
xscalepct       <- as.numeric(as.character(graph_parms[14 ,2]))
yscalepct       <- as.numeric(as.character(graph_parms[15 ,2]))


linespacing <- linespacing * (ymaxbox - yminbox)
deltaxleft  <- deltaxleft * (xmaxbox - xminbox)
deltaytop   <- deltaytop * (ymaxbox - yminbox)
deltaxleft

xmin <- xminbox
xmax <- xmaxbox
ymin <- yminbox
ymax <- ymaxbox

xarrow <- xmin + xarrowpct*(xmax-xmin)
yarrow <- ymin + yarrowpct*(ymax-ymin)
warrow <- scalearrow*(xmax-xmin)

xscaleline <- xmin + xscalepct*(xmax-xmin)
yscaleline <- ymin + yscalepct*(ymax-ymin)
ndivs <- 2
division <- scaleline/ndivs

The next part of the script defines various parameters need for computing the contour lines and setting up the data for plotting

#======================================================================================================
#======================================================================================================
# Compute the spatial intervals

x <- seq(xminbox, xmaxbox, as.integer((xmaxbox-xminbox)/nx +0.5))
y <- seq(yminbox, ymaxbox, as.integer((ymaxbox-yminbox)/ny +0.5))

grdxpoly <- c(xminbox, xmaxbox, xmaxbox, xminbox, xminbox)
grdypoly <- c(yminbox, yminbox, ymaxbox, ymaxbox, yminbox)
grdpoly  <- cbind(grdxpoly, grdypoly)

nxy <- (nx+1)*(ny+1)
grdptsbox <- gridpts(grdpoly, nxy)


#======================================================================================================
# Make a closed polygon of the cove geometry

covelength <- length(xcove)
xcovepoly  <- c(xcove, xcove[1])
ycovepoly  <- c(ycove, ycove[1])
covepoly   <- cbind(xcovepoly, ycovepoly)
p2 <- as(covepoly, "gpc.poly")


#======================================================================================================
# Make a closed polygon of the bounding box geometry

bboxlength <- length(xbox)
xbboxpoly <- c(xbox, xbox[1])
ybboxpoly <- c(ybox, ybox[1])
bboxpoly <- cbind(xbboxpoly, ybboxpoly)
p1 <- as(bboxpoly, "gpc.poly")

The following adds some points to the bathymetry data set to form a better grounding for interpolation. What is needed here is to add the points defined by two surveyed contour lines that can be found on the buy-down maps. These would solidify the anchor points for the bathymetry. Unfortunately these contours are not available in digital form. Digitizing them is a considerable effort.

#======================================================================================================
# Add the shoreline points with 0 depth

xbath <- as.numeric(c(xbath, xcove))
ybath <- as.numeric(c(ybath, ycove))
zbath <- as.numeric(c(zbath, rep(0,length(xcove))))

The following lines define the color and contour palette


#======================================================================================================
# Set the color palette and legend data

breakpoints  <- c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 50, 100)
color_label  <- c("<2 ft","<4 ft", "<6 ft", "<8 ft", "<10 ft", "<12 ft", "<14 ft", "<16 ft", "<18 ft", "<20 ft", "<50 ft", "<100 ft")
#legendcolors <- c("#FFFFFF", "#FF7100", "#FFAA00", "#FFE300", "#4cff00", "#41de00", 
legendcolors <- c("#FF0000", "#FF7100", "#FFAA00", "#FFE300", "#4cff00", "#41de00", 
                  "#32ab00", "#00FFE3", "#00E3FF", "#00AAFF", "#0071FF", "#3900FF")
nlevels      <- length(legendcolors)

The following lines of code do the real work. Akima is a function that interpolates irregularly spaced data onto a grid with regular intervals. It does this interpolation on a convex hull, meaning that the regularly spaced grid is bounded by a polygon that looks like a plastic band wrapped around the geometry of the cove or lake. The side effect is that interpolated results show up on dry land; not desirable.

#======================================================================================================
#======================================================================================================
#======================================================================================================
# Interpolate the values with Akima without contour lines

print("Doing Akima interpolation")
pdf(file=paste("./case_maps/Akima/", project, "_", resolution,"_Akima", ".pdf", sep=""), onefile=TRUE)

akima <- interp(xbath, ybath, zbath, duplicate="strip", xo=seq(xminbox, xmaxbox, length = nx),
          yo=seq(yminbox, ymaxbox, length = ny), linear = TRUE, extrap=FALSE)

Now we do the plotting of the results. Because Akima does its calculation over a convex hull. This has the consequence that interpolated results show up outside of the lake boundaries. These points are masked out by the plot(setdiff(p1, p2)… statement. Discovering this approach took me a long time. Thus what's done here is as follows:

  1. Plot the contours
  2. Plot the locations of the boat slips (two plots, one a square, one a dot in the middle of the square)
  3. Plot the mask to get rid of things between the cove outline and the convex hull outline
  4. Plot the outline of the cove
  5. Plot the points of interest if there are any
  6. Label the points of interest
image(akima, xlim = extendrange(xbox), ylim = extendrange(ybox), breaks=breakpoints,
      col=legendcolors, lwd=0.05, lty=3, xlab="Easting (ft)", ylab="Northing(ft)")

#contour(akima, add=TRUE, labcex=0.7, drawlabels=TRUE)              # contours
points(xdock, ydock, cex=0.7, pch=15, col="white")                  # Boat slips - square outline
points(xdock, ydock, cex=0.3, pch=20, col="black")                  # Boat slips - black dot
plot(setdiff(p1, p2), poly.args = list(col = "grey"), add = TRUE)
lines(xcovepoly, ycovepoly, type="l", lwd=2, col="blueviolet")      # Cove outline
if(npoi != 0){
    points(xpoi, ypoi, cex=1.0, pch=17, col="black")                # Points of Interest - black triangle
    text(xpoi, ypoi, poin, pos=4, offset=0.25, cex=0.65)            # Label the points of interest
}

The next part of the script places the various text elements on the graph. The locations of these elements are read from a file. A couple of iterations, changing manually the location properties in the file, are needed in order to produce a tidy graph. What is being placed are the following items:

  1. The who, with what and when
  2. The interpolation method used
  3. The are of the cove in acres
  4. The length of the shoreline of the cove in miles
  5. The spatial resolution of the map in feet
  6. The number of intervals in both the x and y directions
  7. The legend marker for the oat slips
  8. The reference elevation for the depths, feet above mean sea level.
  9. The arrow pointing north.
  10. The scale of the map
  11. The title of the map
  12. the legend of the colors.
# To add various elements is done on the basis of %x and %y of the location of the lower left-hand corner of the element
# The overall size is uniformly adjusted according to a scale factor that applies to each element

usr <- par( "usr" )
text(usr[1]+deltaxleft, usr[4]-deltaytop, paste("PLV: Image created with R:",Sys.Date()), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-1*linespacing, paste("Interpolations with Akima"), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-2*linespacing, paste("Area of the cove:",areap,"acres."), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-3*linespacing, paste("Shore length of the cove:",milesp,"miles."), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-4*linespacing, paste("Spatial resolution:",resolution,"ft."), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-5*linespacing, paste(nx,"by",ny,"grid."), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-6*linespacing, paste("White squares & black dots - boat slips"), pos=4, cex=fontsize)
text(usr[1]+deltaxleft, usr[4]-deltaytop-7*linespacing, paste("Depths referenced to 2,462 ft ASL"), pos=4, cex=fontsize)


north.arrow(xarrow,yarrow,warrow,lab="NORTH",cex.lab=0.5,tcol="black",col="orange")
map.scale(xscaleline,yscaleline,scaleline,"ft",ndivs,subdiv=division,tcol="black",scol="black",sfcol="black")

title(paste("A Bathymetric Map of ", areaname), font = 4)

legend(legendlocation, color_label, col=legendcolors, lwd=10, title="Contours", cex=fontsize, inset=legendinset)

The next part of the script draws a miniature version of the lake and draws a circle around showing where on the lake the cove is located. The central coordinates of the cove are read from a file that contains the data for all of the coves.

#======================================================================================================
# Plot the lake and the cove in small version

shore <- read.table(paste("./basic_inputs/DCL_shoreline_SPSC_V2.txt", sep="\t"), header=TRUE, sep="\t", comment.char = "#")
xlake <- shore$Easting
ylake <- shore$Northing

# Check to see if closed polygon.  If not make it that way.
if((xlake[1] != xlake[length(xlake)]) || (ylake[1] != ylake[length(ylake)])){
    xlake[length(xlake)+1] = xlake[1]
    ylake[length(ylake)+1] = ylake[1]
}

# Compute the plot scaling factors
xmin <- min(xlake)
xmax <- max(xlake)
ymin <- min(ylake)
ymax <- max(ylake)

scalex <- (xmaxbox-xminbox)/(xmax-xmin)
scaley <- (ymaxbox-yminbox)/(ymax-ymin)


# Read the file with all of the cove locations
fileConn      <- file("./basic_inputs/ListOfCoves.txt", "rt")
covelocations <- read.table(fileConn, header=TRUE, sep="\t", comment.char = "#")
close(fileConn)

M <- length(covelocations$ID)

# Search for the area name to be projected on the lake map -  stored as "k"
narea <- 0
for(j in 1:M){
    if(as.character(covelocations$ID[j]) == project) narea <- j
}
if(narea != 0){
    x00 <- xminbox + minilakex0*(xmaxbox - xminbox)     # x0-coordinate of small box
    y00 <- yminbox + minilakey0*(ymaxbox - yminbox)     # y0-coordinate of small box
    smallcovex <- x00 + scalefactor*scalex*(covelocations$Easting[narea] - xmin)
    smallcovey <- y00 + scalefactor*scaley*(covelocations$Northing[narea] - ymin)

    xbb <- c(xmin, xmax, xmax, xmin, xmin)
    ybb <- c(ymin, ymin, ymax, ymax, ymin)
    lines(x00 + scalefactor*scalex*(xlake-xmin), y00 + scalefactor*scaley*(ylake-ymin), type="l", col="red", lwd=1)     # Lake outline
    lines(x00 + scalefactor*scalex*(xbb-xmin),   y00 + scalefactor*scaley*(ybb-ymin), type="l", col="blue", lwd=1)      # Bounding box


    polygon(x00+scalefactor*scalex*(xbb-xmin), y00+scalefactor*scaley*(ybb-ymin), col="gold")
    polygon(x00+scalefactor*scalex*(xlake-xmin), y00+scalefactor*scaley*(ylake-ymin), col="light blue")
    points(smallcovex, smallcovey, cex=1.2, col="red")      # Circle of cove on the minnie whole lake map

}

The following few statements closes the plotting device and computes the elapsed done since starting this case run. It also list any warnings that were generated by the execution of this script.

dev.off()
#===============================================================================================
begTime
runTime <- Sys.time()-begTime
runTime
print(sprintf("Finished with -%s- in -CoveMaps.R-", project))
warnings()