[R] identifying local maxima

[Hans W Borchers]

Gary Dong <pdxgary163 <at> gmail.com> writes:
> 
> Dear R users,
> 
> I have created a Loess surface in R, in which x is relative longitude by
> miles, y is relative latitude by miles, and z is population density at the
> neighborhood level. The purpose is to identify some population centers in
> the region. I'm wondering if there is a way to determine the coordinates
> (x,y) of each center, so I can know exactly where they are.
> 
> Let me use the "elevation" data set in geoR to illustrate what have done:
> 
> require(geoR)
> 
> data(elevation)
> 
> elev.df <- data.frame(x = 50 * elevation$coords[,"x"], y = 50 *
> elevation$coords[,"y"], z = 10 * elevation$data)
> 
> elev.loess <- loess(z ~ x * y, data = elev.df, degree = 2, span = 0.1)
> 
> grid.x <- seq(10, 300, 1)
> grid.y <- seq(10, 300, 1)
> grid.mar <- list(x=grid.x, y=grid.y)
> elev.fit<-expand.grid(grid.mar)
> 
> z.pred <- predict(elev.loess, newdata = elev.fit)
> 
> persp(seq(10, 300, 5), seq(10, 300, 5), emp.pred, phi = 45, theta = 45,
> xlab = "X Coordinate (feet)", ylab = "Y Coordinate (feet)", main = "Surface
> elevation data")
> 
> With this, what's the right way to determine the coordinates of the local
> maixma on the surface?

If there are more local maxima, you probably need to start the optimization 
process from each point of your grid and see if you stumble into different 
maxima. Here is how to find the one maximum in your example.

    require(geoR); data(elevation)
    elev.df <- data.frame(x = 50 * elevation$coords[,"x"],
                          y = 50 *elevation$coords[,"y"], 
                          z = 10 * elevation$data)

    ##  Compute the surface on an irregular grid
    library(akima)
    foo <- function(xy) with( elev.df, 
            -interp(x, y, z, xy[1], xy[2], linear=FALSE, extrap=TRUE)$z )
    
    ##  Use Nelder-Mead to find a local maximum
    optim(c(150, 150), foo)
    # $par
    # [1] 195.8372  21.5866
    # $value
    # [1] -9705.536
    
    ##  See contour plot if this is the only maximum
    with(elev.df, 
        {A <- akima::interp(x, y, z, linear=FALSE)
         plot(x, y, pch=20, col="blue")
    	 contour(A$x, A$y, A$z, add=TRUE) })
    points(195.8372, 21.5866, col="red")

> Thanks.
> 
> Gary