https://github.com/cran/spatstat
Raw File
Tip revision: e575bb3736f6d70d2bd2b23ea2e4474cdbae00be authored by Adrian Baddeley on 11 November 2010, 13:09:43 UTC
version 1.21-1
Tip revision: e575bb3
ppp.R
#
#	ppp.S
#
#	A class 'ppp' to define point patterns
#	observed in arbitrary windows in two dimensions.
#
#	$Revision: 4.73 $	$Date: 2010/10/28 02:08:37 $
#
#	A point pattern contains the following entries:	
#
#		$window:	an object of class 'owin'
#				defining the observation window
#
#		$n:	the number of points (for efficiency)
#	
#		$x:	
#		$y:	vectors of length n giving the Cartesian
#			coordinates of the points.
#
#	It may also contain the entry:	
#
#		$marks:	a vector of length n
#			whose entries are interpreted as the
#			'marks' attached to the corresponding points.	
#	
#--------------------------------------------------------------------------
ppp <- function(x, y, ..., window, marks, check=TRUE ) {
  # Constructs an object of class 'ppp'
  #
  if(!missing(window))
    verifyclass(window, "owin")
  else
    window <- owin(...)

  if(missing(x) && missing(y)) 
    x <- y <- numeric(0)
  
  n <- length(x)
  if(length(y) != n)
    stop("coordinate vectors x and y are not of equal length")
  
  # validate x, y coordinates
  stopifnot(is.numeric(x))
  stopifnot(is.numeric(y))
  ok <- is.finite(x) & is.finite(y)
  if(any(!ok)) {
    nbg <- is.na(x) | is.na(y)
    if(any(nbg)) {
      howmany <- if(all(nbg)) "all" else paste(sum(nbg),  "out of", length(nbg))
      stop(paste(howmany, "coordinate values are NA or NaN"))
    }
    howmany <- if(!any(ok)) "all" else paste(sum(!ok),  "out of", length(ok))
    stop(paste(howmany, "coordinate values are infinite"))
  }

  names(x) <- NULL
  names(y) <- NULL
  
  # check (x,y) points lie inside window
  if(check && n > 0) {
    ok <- inside.owin(x, y, window)
    nout <- sum(!ok)
    if(nout > 0) {
      warning(paste(nout,
                    ngettext(nout, "point was", "points were"),
                    "rejected as lying outside the specified window"))
      rr <- ripras(x,y)
      bb <- bounding.box.xy(x,y)
      bb <- bounding.box(rr, bb, window)
      rejectwindow <-
        if(!is.null(rr)) rebound.owin(rr, bb) else bb
      rejects <- ppp(x[!ok], y[!ok], window=rejectwindow, check=FALSE)
      # discard illegal points
      x <- x[ok]
      y <- y[ok]
      n <- length(x)
    }
  } else nout <- 0
  # initialise ppp object
  pp <- list(window=window, n=n, x=x, y=y)
  # coerce marks to appropriate forma
  if(missing(marks))
    marks <- NULL
  if(is.hyperframe(marks)) 
    stop("Hyperframes of marks are not implemented for ppp objects; use ppx")
  if(is.matrix(marks)) 
    marks <- as.data.frame(marks)
  if(is.data.frame(marks)) {
    nc <- ncol(marks)
    if(nc == 0)
      marks <- NULL
    else if(nc == 1)
      marks <- marks[,,drop=TRUE]
  }
  # attach marks 
  if(is.null(marks)) {
    # no marks
    pp$markformat <- "none"
  } else if(is.data.frame(marks)) {
    # data frame of marks
    pp$markformat <- "dataframe"
    if(nout > 0) {
      marks <- marks[ok, ]
      marks(rejects) <- marks[!ok,]
    }
    if(nrow(marks) != n)
      stop("number of rows of marks != length of x and y")
    pp$marks <- marks
  } else {
    # should be a vector or factor
    # To recognise vector, strip attributes
    if(!is.factor(marks))
      attributes(marks) <- NULL
    if(!(is.vector(marks) || is.factor(marks)))
      stop("Format of marks not understood")
    # OK, it's a vector or factor
    pp$markformat <- "vector"
    if(nout > 0) {
      marks(rejects) <- marks[!ok]
      marks <- marks[ok]
    }
    if(length(marks) != n)
      stop("length of marks vector != length of x and y")
    names(marks) <- NULL
    pp$marks <- marks
  }
  class(pp) <- "ppp"
  if(check && any(duplicated(pp)))
    warning("data contain duplicated points")
  if(nout > 0) 
    attr(pp, "rejects") <- rejects
  pp
}

#
#--------------------------------------------------------------------------
#

is.ppp <- function(x) { inherits(x, "ppp") }

#
#--------------------------------------------------------------------------
#

as.ppp <- function(X, ..., fatal=TRUE) {
  UseMethod("as.ppp")
}

as.ppp.ppp <- function(X, ..., fatal=TRUE) {
  check <- resolve.defaults(list(...), list(check=FALSE))$check
  return(ppp(X$x, X$y, window=X$window, marks=X$marks, check=check))
}

as.ppp.quad <- function(X, ..., fatal=TRUE) {
  return(union.quad(X))
}

as.ppp.data.frame <- function(X, W = NULL, ..., fatal=TRUE) {
  check <- resolve.defaults(list(...), list(check=TRUE))$check
  if(!(ncol(X) %in% 2:3)) {
    if(fatal)
      stop("X must be a two-column or three-column data frame")
    else
      return(NULL)
  }

  if(is.null(W)) {
    if(fatal)
      stop("x,y coords given but no window specified")
    else
      return(NULL)
  } 
    
  if(is.function(W))
    Z <- cobble.xy(X[,1], X[,2], W, fatal)
  else {
    win <- as.owin(W)
    Z <- ppp(X[,1], X[,2], window = win, check=check)
  }

  # add marks from third column
  if(ncol(X) == 3)
    Z <- Z %mark% X[,3]

  return(Z)
}
    
as.ppp.matrix <- function(X, W = NULL, ..., fatal=TRUE) {
  check <- resolve.defaults(list(...), list(check=TRUE))$check
  if(!verifyclass(X, "matrix", fatal=fatal)
     || !is.numeric(X)) {
    if(fatal)
      stop("X must be a two-column matrix of numbers")
    else
      return(NULL)
  }
  if(!(ncol(X) %in% 2:3)) {
    if(fatal)
      stop("X must be a two-column or three-column matrix of numbers")
    else
      return(NULL)
  }

  if(is.null(W)) {
    if(fatal)
      stop("x,y coords given but no window specified")
    else
      return(NULL)
  } 
    
  if(is.function(W))
    Z <- cobble.xy(X[,1], X[,2], W, fatal)
  else {
    win <- as.owin(W)
    Z <- ppp(X[,1], X[,2], window = win, check=check)
  }

  # add marks from third column
  if(ncol(X) == 3)
    Z <- Z %mark% X[,3]

  return(Z)
}
    
as.ppp.default <- function(X, W=NULL, ..., fatal=TRUE) {
	# tries to coerce data X to a point pattern
	# X may be:
	#	1. a structure with entries x, y, xl, xu, yl, yu
	#	2. a structure with entries x, y, area where
        #                    'area' has entries xl, xu, yl, yu
	#	3. a structure with entries x, y
        #       4. a vector of length 2, interpreted as a single point.
	# The second argument W is coerced to an object of class 'owin' by the 
	# function "as.owin" in window.S
        # If X also has an entry X$marks
        # then this will be interpreted as the marks vector for the pattern.
	#
  check <- resolve.defaults(list(...), list(check=TRUE))$check
  if(checkfields(X, c("x", "y", "xl", "xu", "yl", "yu"))) {
		xrange <- c(X$xl, X$xu)
		yrange <- c(X$yl, X$yu)
		if(is.null(X$marks))
			Z <- ppp(X$x, X$y, xrange, yrange, check=check)
		else
			Z <- ppp(X$x, X$y, xrange, yrange, 
				marks=X$marks, check=check)
		return(Z)
        } else if(checkfields(X, c("x", "y", "area"))
                  && checkfields(X$area, c("xl", "xu", "yl", "yu"))) {
                win <- as.owin(X$area)
                if (is.null(X$marks))
                  Z <- ppp(X$x, X$y, window=win, check=check)
                else
                  Z <- ppp(X$x, X$y, window=win, marks = X$marks, check=check)
                return(Z)
	} else if(checkfields(X, c("x", "y"))) {
                if(is.function(W))
                  return(cobble.xy(X$x, X$y, W, fatal))
		if(is.null(W)) {
                  if(fatal)
                    stop("x,y coords given but no window specified")
                  else
                    return(NULL)
                }
		win <- as.owin(W)
		if(is.null(X$marks))
                  Z <- ppp(X$x, X$y, window=win, check=check)
                else
                  Z <- ppp(X$x, X$y, window=win, marks=X$marks, check=check)
                return(Z)
        } else if(is.vector(X) && length(X) == 2) {
                win <- as.owin(W)
                Z <- ppp(X[1], X[2], window=win, check=check)
                return(Z)
	} else {
          if(fatal)
            stop("Can't interpret X as a point pattern")
          else
            return(NULL)
        }
}

cobble.xy <- function(x, y, f=ripras, fatal=TRUE) {
  if(!is.function(f))
    stop("f is not a function")
  w <- f(x,y)
  if(!is.owin(w)) {
    gripe <- "Supplied function f did not return an owin object"
    if(fatal)
      stop(gripe)
    else {
      warning(gripe)
      return(NULL)
    }
  }
  return(ppp(x, y, window=w))
}
  

# --------------------------------------------------------------

"[.ppp" <-
  function(x, i, j, drop, ...) {

        verifyclass(x, "ppp")
        
        if(missing(i) && missing(j))
          return(x)

        if(!missing(i)) {
          if(inherits(i, "owin")) {
            # i is a window
            window <- i
            ok <- inside.owin(x$x, x$y, window)
            x <- ppp(x$x[ok], x$y[ok], window=window, #SIC
                     marks=marksubset(x$marks, ok),
                     check=FALSE)
          } else if(inherits(i, "im")) {
            # i is an image
            if(i$type != "logical")
              stop(paste("Subset operator X[i] undefined",
                         "when i is a pixel image",
                         "unless it has logical values"), call.=FALSE)
            # convert logical image to window
            e <- sys.frame(sys.nframe())
            window <- solutionset(i, e)
            ok <- inside.owin(x$x, x$y, window)
            x <- ppp(x$x[ok], x$y[ok], window=window, #SIC
                     marks=marksubset(x$marks, ok),
                     check=FALSE)
          } else {
            # assume i is a subset index
            if(x$n == 0)
              return(x)
            subset <- i
            x <- ppp(x$x[subset], x$y[subset], window=x$window,
                     marks=marksubset(x$marks, subset),
                     check=FALSE)
          } 
        }

        if(!missing(j))
          x <- x[j]   # invokes code above

        return(x)
}


# ------------------------------------------------------------------
#
#
scanpp <- function(filename, window, header=TRUE, dir="", multitype=FALSE) {
  filename <- if(dir=="") filename else
              paste(dir, filename, sep=.Platform$file.sep)
  df <- read.table(filename, header=header)
  if(header) {
    # check whether there are columns named 'x' and 'y'
    colnames <- dimnames(df)[[2]]
    xycolumns <- match(c("x", "y"), colnames, 0)
    named <- all(xycolumns > 0)
  } else {
    named <- FALSE
  }
  if(named) {
    x <- df$x
    y <- df$y
  } else {
    # assume x, y given in columns 1, 2 respectively
    x <- df[,1]
    y <- df[,2]
    xycolumns <- c(1,2)
  }
  if(ncol(df) == 2) 
      X <- ppp(x, y, window=window)
  else {
    marks <- df[ , -xycolumns]
    if(multitype) 
      marks <- factor(marks)
    X <- ppp(x, y, window=window, marks = marks)
  }
  X
}

#-------------------------------------------------------------------

"markspace.integral" <-
  function(X) {
  verifyclass(X, "ppp")
  if(!is.marked(X, dfok=TRUE))
    return(1)
  if(is.multitype(X))
    return(length(levels(marks(X))))
  else
    stop("Don't know how to compute total mass of mark space")
}

#-------------------------------------------------------------------

print.ppp <- function(x, ...) {
  verifyclass(x, "ppp")
  ism <- is.marked(x, dfok=TRUE)
  cat(paste(if(ism) "marked" else NULL,
            "planar point pattern:",
            x$n,
            ngettext(x$n, "point", "points"),
            "\n"))
  if(ism) {
    mks <- marks(x, dfok=TRUE)
    if(is.data.frame(mks)) {
      # data frame of marks
      cat(paste("Mark variables: ", paste(names(mks), collapse=", "), "\n"))
    } else {
      # vector of marks
      if(is.factor(mks)) {
        cat("multitype, with ")
        cat(paste("levels =", paste(levels(mks), collapse="\t"),"\n"))
      } else {
        cat(paste("marks are",
                  if(is.numeric(mks)) "numeric,",
                  "of type", sQuote(typeof(mks)), "\n"))
      }
    }
  }
  print(x$window)
  if(!is.null(rejects <- attr(x, "rejects"))) {
    nrejects <- rejects$n
    cat(paste("\n***",
              nrejects,
              ngettext(nrejects, "illegal point", "illegal points"),
              "stored in",
              paste("attr(,", dQuote("rejects"), ")", sep=""),
              "***\n"))
  }
  return(invisible(NULL))
}

summary.ppp <- function(object, ..., checkdup=TRUE) {
  verifyclass(object, "ppp")
  result <- list()
  result$is.marked <- is.marked(object, dfok=TRUE)
  result$n <- object$n
  result$window <- summary(object$window)
  result$intensity <- result$n/result$window$area
  if(checkdup)
    result$nduplicated <- sum(duplicated(object))
  if(result$is.marked) {
    mks <- marks(object, dfok=TRUE)
    if(result$multiple.marks <- is.data.frame(mks)) {
      result$marknames <- names(mks)
      result$is.numeric <- FALSE
      result$marktype <- "dataframe"
      result$is.multitype <- FALSE
    } else {
      result$is.numeric <- is.numeric(mks)
      result$marknames <- "marks"
      result$marktype <- typeof(mks)
      result$is.multitype <- is.multitype(object)
    }
    if(result$is.multitype) {
      tm <- as.vector(table(mks))
      tfp <- data.frame(frequency=tm,
                        proportion=tm/sum(tm),
                        intensity=tm/result$window$area,
                        row.names=levels(mks))
      result$marks <- tfp
    } else 
      result$marks <- summary(mks)
  }
  class(result) <- "summary.ppp"
  if(!is.null(rejects <- attr(object, "rejects"))) 
    result$rejects <- rejects$n
  return(result)
}

print.summary.ppp <- function(x, ..., dp=3) {
  verifyclass(x, "summary.ppp")
  cat(paste(if(x$is.marked) "Marked planar " else "Planar ",
            "point pattern: ",
            x$n,
            " points\n",
            sep=""))
  oneline <- resolve.defaults(list(...), list(oneline=FALSE))$oneline
  if(oneline) return(invisible(NULL))
  unitinfo <- summary(x$window$units)
  cat(paste("Average intensity",
            signif(x$intensity,dp),
            "points per square",
            unitinfo$singular,
            unitinfo$explain,
            "\n"))
  ndup <- x$nduplicated
  if((!is.null(ndup)) && (ndup > 0))
    cat("\n*Pattern contains duplicated points*\n")
  if(x$is.marked) {
    if(x$multiple.marks) {
      cat(paste("Mark variables: ", paste(x$marknames, collapse=", "), "\n"))
      cat("Summary:\n")
      print(x$marks)
    } else if(x$is.multitype) {
      cat("Multitype:\n")
      print(signif(x$marks,dp))
    } else {
      cat(paste("marks are ",
                if(x$is.numeric) "numeric, ",
                "of type", sQuote(x$marktype), "\n"))
      cat("Summary:\n")
      print(x$marks)
    }
  }
  cat("\n")
  print(x$window)
  if(!is.null(nrejects <- x$rejects)) 
    cat(paste("\n***",
              nrejects,
              ngettext(nrejects, "illegal point", "illegal points"),
              "stored in",
              paste("attr(,", dQuote("rejects"), ")", sep=""),
              "***\n"))
  return(invisible(x))
}

# ---------------------------------------------------------------

identify.ppp <- function(x, ...) {
  verifyclass(x, "ppp")
  if(!is.marked(x) || "labels" %in% names(list(...)))
    identify(x$x, x$y, ...)
  else {
    marx <- marks(x, dfok=FALSE)
    marques <- if(is.numeric(marx)) paste(signif(marx, 3)) else paste(marx)
    id <- identify(x$x, x$y, labels=marques, ...)
    mk <- marx[id]
    if(is.factor(marx)) mk <- levels(marx)[mk]
    cbind(id=id, marks=mk)
  }
}

rebound <- function(x, rect) {
  UseMethod("rebound")
}

rebound.ppp <- function(x, rect) {
  verifyclass(x, "ppp")
  x$window <- rebound.owin(x$window, rect)
  return(x)
}

as.data.frame.ppp <- function(x, row.names=NULL, ...) {
  df <- data.frame(x=x$x, y=x$y, row.names=row.names)
  marx <- marks(x, dfok=TRUE)
  if(is.null(marx))
    return(df)
  if(is.data.frame(marx))
    df <- cbind(df, marx)
  else
    df <- data.frame(df, marks=marx)
  return(df)
}

is.empty.ppp <- function(x) { return(x$n == 0) }

npoints <- function(x) {
  UseMethod("npoints")
}

npoints.ppp <- function(x) { x$n }
back to top