[GIS] R using fortify to project shapefile on google map gives wrong results

ggmapggplot2rshapefile

I have a shapefile that looks like this

> str(FieldsMap)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 3 obs. of  2 variables:
  .. ..$ id  : int [1:3] 3 2 1
  .. ..$ Attr: Factor w/ 3 levels "A","B","C": 3 2 1
  ..@ polygons   :List of 3
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 8.6 47.7
  .. .. .. .. .. .. ..@ area   : num 5.88e-08
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 8.6 8.61 8.61 8.61 8.6 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 8.6 47.7
  .. .. .. ..@ ID       : chr "0"
  .. .. .. ..@ area     : num 5.88e-08
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 8.61 47.68
  .. .. .. .. .. .. ..@ area   : num 1.41e-06
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:15, 1:2] 8.61 8.61 8.61 8.61 8.61 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 8.61 47.68
  .. .. .. ..@ ID       : chr "1"
  .. .. .. ..@ area     : num 1.41e-06
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 8.61 47.68
  .. .. .. .. .. .. ..@ area   : num 1.11e-06
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:17, 1:2] 8.61 8.61 8.61 8.61 8.61 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 8.61 47.68
  .. .. .. ..@ ID       : chr "2"
  .. .. .. ..@ area     : num 1.11e-06
  ..@ plotOrder  : int [1:3] 2 3 1
  ..@ bbox       : num [1:2, 1:2] 8.6 47.68 8.61 47.68
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

No I want to plot the polygons with ggmap on a google map so I run:

MapExtent<-bbox(FieldsMap)
Map<-get_map(MapExtent,maptype="hybrid")
mapTransform<-fortify(FieldsMap,region="id")
p1<-ggmap(Map)+geom_polygon(data=mapTransform,aes(long,lat,group=as.factor(group)),fill="green",color="green")
p1

However, the polygons do not have the right shape after fortifying and are not exactly at the right place( a few hundred meters of). Does anyone has an idea on how to solve this? Plotting plot(FieldsMap) gives the right shape of the polygon. I also noticed that after fortifying there is sometimes no polygon plotted using another shapefile. Displaying both shapefiles in QGIS works as it should, so I supect fortify to cause the problem, because looking at the plot I get the impression that sometimes points are missing or are connected in the wrong order.

Best Answer

If your polygons do not overlap the map properly, you have a Coordinate Reference System (CRS) issue. If your map provider is Google Maps, the CRS is WGS84. Thus, you need to change the CRS of your SpatialPolygonDataFrame.

library(sp)
library(rgdal)
proj4string(FieldsMap) # Gives the CRS of the data 
FieldsMap=spTransform(FieldsMap,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")) # To convert it to WGS84

Then use ggplot2::fortify(FieldsMap,byid="your id var") to convert your SpatialPolygonDataFrame into a DataFrame, which you can use as input for ggplot. Beware, fortify() sometimes distort the polygons. If you want a simple graph, prefer the standard plot() function.

Related Question