0

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.

1
  • When posting a question please do not imply that the software/library is producing incorrect results unless you have a reproducible example that illustrates a potential bug. Commonly, ~99% of the time an OP implies that the software is incorrect it is user error. Developers work hard, and often for free, and do not need users intimating that their code is producing incorrect results unless there is ample proof of a bug that can be addressed. – Jeffrey Evans Jul 1 '15 at 15:15
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.