Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'd like to plot the distribution of a species on a map of Switzerland, following this tutorial.

library(maps)
library(mapdata)

map('worldHires', 'Switzerland', xlim=c(5, 12), ylim=c(45,48))

From this website I download the shapefile for e.g. Pinus sylvestris

psylvestris <- readShapePoly("Pinus sylvestris.shp")

I want to add the range to the map:

plot(psylvestris, add=TRUE, xlim=c(5, 12), ylim=c(45,48))

But I do not see anything (other than the map).. It appears to me the coordinates might be given in a different system (?)

> psylvestris
class       : SpatialPolygonsDataFrame 
features    : 52 
extent      : -1202142, 4710481, -925131.2, 6334053  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 1
names       : Id 
min values  :  0 
max values  :  0 

Unfortunately I do not know anything about spatial analysis and thus I have no idea how to fix this. Can you help me please?

share|improve this question

1 Answer 1

up vote 2 down vote accepted

You have to re-project the shapefile to LatLong/WGS84 (spTransform()) :

library("maps")
library("mapdata")
library("rgdal")

psylvestris <- readOGR("Pinus Sylvestris.shp", "Pinus Sylvestris")
proj4string(psylvestris)
# [1] "+proj=laea +lat_0=48 +lon_0=9 +x_0=0 +y_0=0 +a=6378388 +b=6378388 +units=m +no_defs"

psylvestris <- spTransform(psylvestris, CRS("+proj=longlat +datum=WGS84"))
proj4string(psylvestris) 
# [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

plot(psylvestris, xlim=c(5, 12), ylim=c(45,48), axes=TRUE)
map("worldHires", "Switzerland", add=TRUE, col="red", lwd=3)

map

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.