R help - Jul 2012 - identifying local maxima

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?

Thanks.

Gary

	[[alternative HTML version deleted]]

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