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?
Answer
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)

No comments:
Post a Comment