identifying local maxima

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

identifying local maxima

wudadan
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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: identifying local maxima

wudadan
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
<[hidden email]>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
>
> ______________________________________________
> [hidden email] 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]]

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