5

I imported a Shapefile holding a quite detailed world map from: http://www.diva-gis.org

library(sp)
library(rgdal)
world = readOGR('./countries_shp/', 'countries')

Using R with the sp-Packages, plotting is really easy:

spplot(world, 'UNREG2')
plot(world)

Now I want to

  • draw a customizable grid
  • add points derived from coordinates

For example importing coordinates of Vienna (48°12′32″N 16°22′21″E) can be done using char2dms of the sp-package:

vienna.N = char2dms("48d12'32\"N")
vienna.E = char2dms("16d22'21\"E")

But how can a draw a point at this exact location on my map?


Also adding a grid is straight forward:

plot(gridlines(world), add=TRUE)

assuming, that the last plot was created with "plot(world)". How can this be achieved with the spplot output?

underdark
  • 84,148
  • 21
  • 231
  • 413
FloE
  • 241
  • 2
  • 10
  • 1
    Just posting a link to some useful R tricks from a few months ago here: http://gis.stackexchange.com/questions/3310/what-is-the-most-useful-spatial-r-trick You can plot points with dismo and ggplot2 (I don't know how to do it in sp). – djq Apr 05 '12 at 14:43
  • 1
    http://r-spatial.sourceforge.net/gallery/#fig06.R gives an example. You can't mix base graphics with spplot (grid) if that is what your asking. – Andy W Apr 05 '12 at 14:59
  • Thanks for the fast answers! I do not want to mix base and sp (lattice) graphics, but draw a customizable grid over/on an spplot. It is very easy using base-graphics (plot) so there just has to be a similar way for the sp package? – FloE Apr 05 '12 at 15:24
  • This sounds like a lattice question, pure and simple: you might get a fast authoritative answer over on SO. – whuber Apr 05 '12 at 15:45
  • I will try to get an answer there, but I think that the the transformation/projection of the grid is actually a GIS related topic... – FloE Apr 05 '12 at 16:14
  • Please try to use meaningful thread titles. – underdark Apr 05 '12 at 16:42

2 Answers2

5

I figured out how to print a point for Vienna on the map.

Importing the coordinates from WikiPedia into R:

vienna.N = char2dms("48d12'32\"N")
vienna.E = char2dms("16d22'21\"E")

Building a sp-object

vienna.coord = data.frame(
        name='Vienna', 
        x=as.numeric(vienna.E), 
        y=as.numeric(vienna.N)
)
coordinates(vienna.coord) = c('x', 'y')

The projection is taken from http://toolserver.org/~geohack/geohack.php?pagename=Vienna&params=48_12_32_N_16_22_21_E_type:city_region:AT where the WikiPedia article links to:

proj4string(vienna.coord) = '+proj=latlong +datum=WGS84'

Printing this object gives the numeric representation which equals exactly the data from mentioned website:

> vienna.coord
         coordinates   name
1 (16.3725, 48.2089) Vienna

The resulting

> class(vienna.coord)[1]
[1] "SpatialPointsDataFrame"

can be used with base- and sp-graphics. First we transform the projection to the one used in the Diva worldmap:

vienna.coord = spTransform(vienna.coord, CRS(proj4string(world)))

Then a plot for Europe can be produced:

spplot(world, 'COUNTRY', colorkey=FALSE, 
    xlim=c(-24.5, 69.5), ylim=c(26.6, 71.2),
    sp.layout = c('sp.points', vienna.coord, col='red', pch=16),
    par.settings = list(panel.background=list(col='lightblue')),
    col.regions = rep('white', length(levels(world$COUNTRY)))
)

world map with red dot for Vienna

FloE
  • 241
  • 2
  • 10
0

The part of my questions concerning 'grids' solved itself.

Obviously, the function gridlines generates a rectangular grid for the given projection. A distortion can only be seen after projecting the grid. This might help using 2 projections in one map.

FloE
  • 241
  • 2
  • 10