This issue has only been raised in the community once, two years ago, to the best of my knowledge. I think it's high time to raise the issue again!
readOGR
, from the rgdal package, cannot handle multipoint shapefiles, even though they load fine in other GIS packages like QGIS. Of course, you can save the file as a single part object (see image below) in an external program, but ideally it would also be possible to do it in R. Please see worked example below.
download.file("http://www.personal.leeds.ac.uk/~georl/egs/lnd-stations-multi.zip",
"lnd-stations.zip") # download multi part polygon
unzip("lnd-stations.zip") # unzip
lndS <- readOGR(".", "lnd-stations", p4s = "+init=epsg:27700") # load
OGR data source with driver: ESRI Shapefile
Source: ".", layer: "lnd-stations"
with 2532 features and 27 fields
Feature type: wkbMultiPoint with 2 dimensions
Error in readOGR(".", "lnd-stations", p4s = "+init=epsg:27700") :
Incompatible geometry: 4
After performing "Multipart to singleparts" in the QGIS Vector menu (illustrated above) and saving the file, it loads fine in R. Try it:
download.file("http://www.personal.leeds.ac.uk/~georl/egs/lnd-stns.zip",
"lnd-stns.zip")
unzip("lnd-stns.zip")
lndS <- readOGR(".", "lnd-stns", p4s = "+init=epsg:27700")
plot(lndS)
Answer
Seems to be working with maptools. Loading the fixed shapefile as per your code:
download.file("http://www.personal.leeds.ac.uk/~georl/egs/lnd-stns.zip", "lnd-stns.zip")
unzip("lnd-stns.zip")
lndS.fix <- readOGR(".", "lnd-stns", p4s = "+init=epsg:27700")
plot(lndS.fix)
OGR data source with driver: ESRI Shapefile
Source: ".", layer: "lnd-stns"
with 2532 features and 27 fields
Feature type: wkbPoint with 2 dimensions
What does its size and structure look like?
str(lndS.fix)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
..@ data :'data.frame': 2532 obs. of 27 variables:
.. ..$ CODE : int [1:2532] 5520 5520 5520 5520 5520 5520 5520 5520 5520 5520 ...
Compare to the original multipoint data, but loaded with readShapePoint (which loads with no errors):
library(maptools)
download.file("http://www.personal.leeds.ac.uk/~georl/egs/lnd-stations-multi.zip", "lnd-stations.zip") # download multi part polygon
unzip("lnd-stations.zip") # unzip
lndS <- readShapePoints("lnd-stations.shp")
#Look at it
str(lndS)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
..@ data :'data.frame': 2532 obs. of 27 variables:
.. ..$ CODE : int [1:2532] 5520 5520 5520 5520 5520 5520 5520 5520 5520 5520 ...
Seems to have loaded OK - but are they really identical?
identical(lndS, lndS.fix)
FALSE
No ... but they have the same variables and the same set of XY coordinates:
identical(lndS@coordinates[[2]], lndS.fix@coordinates[[2]]
TRUE
No comments:
Post a Comment