4

I have a dataset with lat/longs mostly in the UK and want to add a column indicating parliamentary constituency. The whole dataset is ~50k rows by 30 columns. A short sample is:

    lat      long Age city_id

1  50.78839 -3.648326  26    1721
2  51.47254 -2.500891  19    1724
3  52.76901 -1.202122  27    1727
4  53.57782 -1.002811  27    1730
5  54.66645 -5.913355  45    1733
6  51.06004 -1.309250  27    1736
7  51.50796 -2.606150  23    1724
8  51.47989 -1.442205  44    1739
9  53.72372 -1.862787  55    1742
10 50.72398 -1.788323  31    1745

I've downloaded shapefiles for Westminster parliamentary constituencies from the ONS: https://data.gov.uk/dataset/westminster-parliamentary-constituencies.

I've then used over in the sp package in R to check which points fall within which polygons. The code looks like:

WMCs<-readOGR(file.choose())
data_SPDF<-SpatialPointsDataFrame(data[,c("long","lat")],data,CRS("+proj=longlat +datum=WGS84"))
data_SPDF<-spTransform(data_SPDF,CRS(proj4string(WMCs)))
data_over_WMCs<-sp::over(data_SPDF,WMCs)
data$WMC<-data_over_WMCs$NAME

This works for most of the data, but around 8k points just get NA. For those first ten lines:

        lat      long Age city_id              WMC
1  50.78839 -3.648326  26    1721             <NA>
2  51.47254 -2.500891  19    1724        Kingswood
3  52.76901 -1.202122  27    1727     Loughborough
4  53.57782 -1.002811  27    1730             <NA>
5  54.66645 -5.913355  45    1733             <NA>
6  51.06004 -1.309250  27    1736             <NA>
7  51.50796 -2.606150  23    1724             <NA>
8  51.47989 -1.442205  44    1739          Newbury
9  53.72372 -1.862787  55    1742          Halifax
10 50.72398 -1.788323  31    1745 Bournemouth East

These points do fall within the UK and within parliamentary constituencies. The ones that fail don't generally fall at the boundaries between constituencies - they're mostly inside, and seem randomly interspersed with points that succeed. Identical code with different shapefiles has allowed me to successfully add a column with MLSOAs.

Any ideas what could be going wrong?

4

There is no need to reproject your data with spTransform as, both datasets are in the same projection. The WGS84 datum will assume a WGS84 spheroid/ellipse and there are no additional parameters defined. Just use proj4string(data_SPDF) <- proj4string(WMCs) to match the projections so over will not throw an error.

Based on your abbreviated example data, you do, in fact, have points that fall outside your polygons. Using over, to return an index of intersecting polygons, indicates that the fifth point does not have a matching polygon.

sp::over(data_SPDF, sp::geometry(WMCs))    
    1   97
    2   252 
    3   270 
    4   145  
    5   NA 
    6   516  
    7   71 
    8   306 
    9   197  
    10  57

Plotting the data illustrates the same thing.

plot(WMCs)
  plot(data_SPDF, pch=20, col="red", cex=1.5, add=TRUE)

I am just not seeing the inconsistency that you are implying and it looks like error in your code in piping mismatched data. I would recommend an approach like this which will return the joined data as a SpatialPointsDataFrame object and drop the NA observations where no points intersect polygons.

data_SPDF <- data_SPDF[!is.na(sp::over(data_SPDF, sp::geometry(WMCs))), ]
( data_SPDF@data <- data.frame(data_SPDF@data, stats::na.omit(sp::over(data_SPDF, WMCs))) )
  • Thank you, that's really helpful! I hadn't realised I didn't need to reproject - not doing that does indeed seem to make it work (i.e., it only gives NAs for points which really are outside the shapes). Do you know what it was about reprojecting the data that made some points appear to be outside the shapes? Anyway, thank you! – tzirtzi May 8 '17 at 18:53

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.