Thanks, Hans. Much appreciated!

I'm also wondering, if a population center (local maxima) is defined as
a contiguous cluster of neighborhoods with population density that is
significantly higher than the neighborhoods surrounding it, will things
become easier? I mean, the identification of local maxima. Thanks.

Gary


On Thu, Jul 12, 2012 at 5:26 AM, Hans W Borchers
<hwborchers@googlemail.com>wrote:

> 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
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

	[[alternative HTML version deleted]]

