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

[Roger Bivand]

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:
> ==> Downloadable map data are available on:
>      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)
Loading required package: foreign
Loading required package: sp
> nc <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1])
> 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))
> plot(xypts, add=TRUE, pch=19, cex=0.2)
> 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)
> 
> ger <- read.shape("germany.shp")
> 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
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no