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.
