Saturday, 30 January 2016

vector - Loading a multipoint shapefile in R


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.


multipart shp file in qgis


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...