[R] Using R map data to determine associated state for a coordinate?

Roger Bivand Roger.Bivand at nhh.no
Thu Sep 8 13:08:19 CEST 2005

```On Thu, 8 Sep 2005, Thomas Petzoldt wrote:

> Werner Wernersen wrote:
> > Hi!
> >
> > I have no idea if this is maybe an easy task utilizing
> > R since I read there is
> > geographical map data in some package:
> >
> > I have a huge number of geographical points with their
> > coordinates in Germany.
> > Now I want to determine for each point in which
> > "Bundesland" = state it is located.
> >
> > Can anybody tell me if this is doable relatively easy
> > in R and if so give me
> > some hints or links how to do it?
> >
> > Thanks a million,
> >    Werner
>
> Hello Werner,
>
> two building blocks, but don't know if the precision meets your needs.
>
> 1. Do you have a good map of Germany *with* Federal States?
>
> * If YES and if it's free:
> ==> I would be interested! Please post it's source.
>
> * If NO:
>      http://www.vdstech.com/map_data.htm
>
> 2. The following approach reads and converts a shapefile with functions
> from maptools and then follows the example of inside.owin() from the
> spatstat package.

The example code will work very well, but since yesterday, when we
released a new version of maptools depending on the sp package, it can
look like:

> library(maptools)
> plot(nc, lwd=2, border="grey")
> bbox(nc)
min       max
r1 -84.32385 -75.45698
r2  33.88199  36.58965
> x <- runif(1000, min=-84.32385, max=-75.45698)
> y <- runif(1000, min=33.88199, max=36.58965)
> xypts <- SpatialPoints(cbind(x, y))
> where_am_i <- overlay(xypts, nc)
> plot(xypts[is.na(where_am_i),], add=TRUE, pch=19, cex=0.2, col="grey80")
> summary(where_am_i)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
1.00   29.00   51.00   53.11   77.00  100.00  450.00

and in this case the points would be read into the SpatialPoints object
directly. Should overlay have trouble with the "huge" number of points,
you could take them in smaller batches, storing the intermediate results.
As Thomas said, you need the boundaries of the "Bundesland" first, and the
accuracy of your results will depend on the degree of detail of the
boundary polygons.

Roger Bivand

>
> Hope that helps
>
> Thomas Petzoldt
>
>
> ##########################################################
> library(maptools)
> library(spatstat)
>
> plot(ger)
>
> pger <- Map2poly(ger)
> sx<- pger[[13]]
> lines(sx, type="l", col="red") # Saxony ;-)
>
> ## Create an owin (observation window) object
> # direction of coordinates must be reversed, in some cases
> # if error message: remove rev()'s
> saxony <- owin(poly=list(x=rev(sx[,1]), y=rev(sx[,2])))
>
> # random points in rectangle
> x <- runif(1000, min= 6, max=15)
> y <- runif(1000, min=46, max=56)
>
> ok <- inside.owin(x, y, saxony)
>
> points(x[ok], y[ok])
> points(x[!ok], y[!ok], pch=".")
> ##########################################################
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
>

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of