How to create contour lines on a map chart?

How to create contour lines on a map chart?

How to create a more detailed map chart using the contourLines() function in TERR.

Attached (Filename: ContourLines.zip) is an example implementation of the contourLines() function. It includes Spotfire file "TERR contourLines sample.dxp", the sample source data file, TERR script file, as well as a screen shot of the final image.  

Here is the data function:

# ===============================================================
# Data function [TERR] Contour Lines
geospatialPredictTimeStamp =  "Thu Aug 28 17:45:08 2014"
# Inputs
#      x = numeric column (x coordinate)          Real([DataTable].[${xcolumn}])
#      y = numeric column (y coordinate)          Real([DataTable].[${ycolumn}])
#      z = numeric column (value to interpolate)  Real([DataTable].[${zcolumn}])
# Optional inputs
#      smooth.scale = numeric scalar, controlling level of detail
#      trimHull (0/1, controls whether to plot regions outside perimeter of points)
#   clipToData (1) - truncate surface to original data limits?
# Output
#   contourTable
#


TimeStamp = paste(date(),Sys.timezone())
datafun = '[TERR] Contour Lines'
tdir = 'C:/Demo' # place to store diagnostics if it exists (otherwise do nothing)
suppressWarnings(try(save(list=ls(), file=paste(tdir,'/ContourLines.in.RData',sep=''), RFormat=T )))
# remove(list=ls()); load( file = 'C:/Demo/ContourLines.in.RData' ); print(datafun); print(TimeStamp)


# x = runif(n=100, min=20, max=30); y = runif(n=100, min=20, max=30); z = runif(n=100, min=20, max=30)
# smooth.scale = 0.5; trimHull=1; clipToData=1
#-----------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------
#------- Define function for clip contour by convex hull ---------------------------
InOut.contour.points = function(x.cont0, y.cont0, xref.vec, yref.vec){
  # Peter Shaw
  # Compares the set of x,y points (xref.vec,yref.vec) with the
  # (non-closed) contour defined by (x.cont0,y.cont0), to evaluate whether
  # the points are inside, or outside the contour.  (The contour is first
  # closed upon itself to form a continuous loop).
 
  # Count the points.
  N.cont = length(x.cont0)
  N.ref  = length(xref.vec)
 
  # Close up the contours to create the working contour definitions.
  # First and last points will be the same.   
  x.cont = c( x.cont0, x.cont0[1] )
  y.cont = c( y.cont0, y.cont0[1] )
 
  u = 1:N.cont
 
  # Create pairs of vectors (vec1 and vec2) from reference points to
  #  movements along the contour (defining pie-shaped sections thus angles)
  vec1.x = rep(1,N.ref) %o% x.cont[u]   - xref.vec %o% rep(1,N.cont)
  vec2.x = rep(1,N.ref) %o% x.cont[u+1] - xref.vec %o% rep(1,N.cont)
  vec1.y = rep(1,N.ref) %o% y.cont[u]   - yref.vec %o% rep(1,N.cont)
  vec2.y = rep(1,N.ref) %o% y.cont[u+1] - yref.vec %o% rep(1,N.cont)
 
  # Form vector cross product and dot products:
  Vec1xVec2 = vec1.x * vec2.y - vec2.x * vec1.y
  Vec1dVec2 = vec1.x * vec2.x + vec1.y * vec2.y
 
  # 4-quadrant atan (angles may be > 90 degrees)
  Ang = atan2( Vec1xVec2, Vec1dVec2 )
 
  # Form angle sums. These should be 0 or +/-1 to roundoff.
  # Sign tells which direction around you went.  May be useful
  angle.sums = rowSums(Ang) / (2 * pi)
 
  # Round these to a reasonable roundoff and take abs
  return( round( abs( angle.sums ) * 1000 ) / 1000 )
}
#-----------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------

if(is.null(smooth.scale)) smooth.scale = 0.5
if(is.null(trimHull)) trimHull = 1
if(is.null(clipToData)) clipToData = 1


#  ------------ Make null objects -------------------
contourTable = data.frame(
  segment = character(0),
  level   = character(0),
  x       = numeric(0),
  y       = numeric(0),
  z       = numeric(0),
  order   = numeric(0),
  stringsAsFactors = F
)

library(RinR)

errorCondition = FALSE
if( length(x)==0 | length(y)==0 | length(z)==0 ) errorCondition=TRUE
if( !errorCondition ) { 
  #xydata.pred.in = data.frame( ID=ID, x=x, y=y, z.orig=z, stringsAsFactors=F ) # for prediction.
  #bad = is.na(x) | is.na(x)
  #xydata.pred.in = xydata.pred.in[ !bad, ]
 
  xyzdata = data.frame( x=x, y=y, z=z )
  bad = is.na(x) | is.na(x) | is.na(z)
  xyzdata = xyzdata[ !bad, ]

  #xyzdata.lo = try(loess(z ~ x * y, data=xyzdata, span=smooth.scale, na.action=na.exclude))
  xyzdata.lo = REvaluate(
    expr = { try(loess(z~x*y, data=data, span=span, na.action=na.exclude)) },
    data = list(data=xyzdata, span=smooth.scale)
  )
 
  if(class(xyzdata.lo)=='Error') errorCondition=TRUE #2nd way to get here
  if(!errorCondition){
    xvec = seq(from=min(x,na.rm=T), to=max(x,na.rm=T), length=200) #Lon.vec
    yvec = seq(from=min(y,na.rm=T), to=max(y,na.rm=T), length=200) #Lat.vec
    xyzdata.lo.pred = REvaluate(
      expr = { predict( xyzdata.lo, newdata=expand.grid(x=xvec, y=yvec)) },
      data = list(xyzdata.lo=xyzdata.lo, xvec=xvec, yvec=yvec)
    )
    x.mtx = xvec %o% rep(1,length(yvec))
    y.mtx = rep(1, length(xvec)) %o% yvec
   
    # calculate convex hull
    points.chull = REvaluate(
      expr = { chull(xyzdata$x, xyzdata$y) }, #calculate convex hull of data points
      data = list(xyzdata=xyzdata)
    )
   
    if( trimHull==1 ) { # Possibly trim points:
      mtx.inout.vec = InOut.contour.points(
        x.cont0 = xyzdata$x[points.chull],
        y.cont0 = xyzdata$y[points.chull],
        xref.vec = as.numeric(x.mtx),
        yref.vec = as.numeric(y.mtx)
      )
      xyzdata.lo.pred[ mtx.inout.vec==0 ] = NA
    }
   
    if( clipToData==1 ) { # Possibly limit surface to original data limits.
      xyzdata.lo.pred = pmin(xyzdata.lo.pred,max(z,na.rm=T))
      xyzdata.lo.pred = pmax(xyzdata.lo.pred,min(z,na.rm=T))
    }
  }
} # end of last if(!errorCondition)

if( !errorCondition ) { # break things up
  # -------------------------------- contour table work -------------------------------------
  # calculate convex hull
#  points.chull = REvaluate(
#    expr = {chull(x, y)}, #calculate convex hull of data points
#    data = list(x=x, y=y)
#  )

#  x.cont0 = xyzdata$x[points.chull],
#  y.cont0 = xyzdata$y[points.chull],
 
 
 
  # Set up data frame with the chull
  points.chull.closed = c(points.chull,points.chull[1])
  x.chull = xyzdata$x[points.chull.closed]
  y.chull = xyzdata$y[points.chull.closed]
  N = length(points.chull.closed)
  borderTable = data.frame(
    segment = rep('border', N),
    level = rep('border', N),
    x   = x.chull + 0.0001*runif(n=N, min = -1, max=1), #helps SF ordering
    y   = y.chull + 0.0001*runif(n=N, min = -1, max=1),
    z   = rep(NA, N),
    stringsAsFactors = F
  )
 
  # Minor Contours
  contourList = REvaluate(
    expr = { contourLines(x=x, y=y, z=z, levels=levels) },
    data = list(x=xvec, y=yvec, z=xyzdata.lo.pred, levels=pretty(xyzdata.lo.pred, 20))
  )
  if( length(contourList) > 0 ) {
    contourTable.pure = do.call(rbind,
                                lapply(X=1:length(contourList),
                                       FUN=function(iseg,contourList){
                                         this.seg=contourList[[iseg]]
                                         N=length(this.seg$x)
                                         table.part = data.frame(
                                           segment  = rep( paste(this.seg$level, iseg, sep='-'), N ),
                                           level    = rep("Minor", N),
                                           x      = this.seg$x,
                                           y      = this.seg$y,
                                           z    = rep(this.seg$level, N),
                                           stringsAsFactors = F
                                         )
                                       }, contourList = contourList
                                )
    )
    contourTable.minor = rbind(borderTable, contourTable.pure)
  } else {
    contourTable.minor = borderTable
  }
 
  # Major Contours
  contourList = REvaluate(
    expr = { contourLines(x=x, y=y, z=z, levels=levels) },
    data = list(x=xvec, y=yvec, z=xyzdata.lo.pred, levels=pretty(xyzdata.lo.pred, 5))
  )
  if(length(contourList)>0){
    contourTable.pure = do.call(rbind,
                                lapply(X=1:length(contourList),
                                       FUN=function(iseg,contourList){
                                         this.seg=contourList[[iseg]]
                                         N=length(this.seg$x)
                                         table.part = data.frame(
                                           segment  = rep( paste(this.seg$level,iseg, sep='-'), N),
                                           level    = rep("Major", N),
                                           x      = this.seg$x,
                                           y      = this.seg$y,
                                           z    = rep(this.seg$level, N),
                                           stringsAsFactors = F
                                         )
                                       }, contourList = contourList
                                )
    )
    contourTable.pure$x = contourTable.pure$x + 0.0001*runif(nrow(contourTable.pure),min=-1,max=1)
    contourTable.pure$y = contourTable.pure$y + 0.0001*runif(nrow(contourTable.pure),min=-1,max=1)
    contourTable = rbind(contourTable.minor, contourTable.pure)
  } else {
    contourTable = contourTable.minor
  }
 
 
  contourTable$order = 1:nrow(contourTable)
} # end of errorCondition

# One last check: if incoming z is constant, no contours will be present.
# Spotfire doesn't like all-missings in the contourTable, so replace the NAs in the border with mean of z if present.
if(all(is.na(contourTable$z))) {
  contourTable$z = mean(z,na.rm=T) # will be NaN if all of z is NA
}

# ===============================================================