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 <[email protected]>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 > > ______________________________________________ > [email protected] 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]] ______________________________________________ [email protected] 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.
