Skip to content

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

3 messages · Werner Wernersen, Thomas Petzoldt, Roger Bivand

Original
Werner Wernersen
# 
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
Thomas Petzoldt
#
Werner Wernersen wrote:
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.

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=".")
##########################################################
Roger Bivand
#
On Thu, 8 Sep 2005, Thomas Petzoldt wrote:

            

      
      
      
    

        
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:

            

      
      
      
    

        
Loading required package: foreign
Loading required package: sp

            

      
      
      
    

        
min       max
r1 -84.32385 -75.45698
r2  33.88199  36.58965

            

      
      
      
    

        
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