1

I'm interested in undertaking spatial matching in R with use of the package on some of the publicly available data consisting of:

Problem

The challenge seems to be with setting up the correct CRS for both of those data sets.

Illustration

I've initially attempted to set up projections to those discussed here, via:

shpNUTS <- spTransform(
  shpNUTS,
  CRSobj = CRS(
    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
  )
)
proj4string(pointsDta) <- proj4string(shpNUTS)
res <- over(pointsDta, shpNUTS)

results in an unusable map:

unusable polygons

Similarly, the attempt to remove projection via:

proj4string(shpNUTS) <- CRS(as.character(NA))

Does not improve the map.

Matching

The matching is undertaken in a simple manner,

sp::over(pointsDta, sp::geometry(shpNUTS))

but I've also tried:

res2 <- spatialEco::point.in.poly(as(pointsDta, "SpatialPointsDataFrame"), shpNUTS)

The matching does not return overlapping points; this is obviously wrong as x/y coordinates sourced from Doogal'swebsite reflect UK postcodes so all of the postcodes should be matched against a set of UK NUTS geographies.


Remaining code

The is read via the readOGR function available through package, Doogal's data is read via readCSV in R. Remaining code reflects basic transformations:

pointsDta <- doogalDta[, c("x", "y")]
pointsDta <- na.omit(pointsDta)
coordinates(pointsDta) <- ~ x + y
# reproject nuts
shpNUTS <- spTransform(
  shpNUTS,
  CRSobj = CRS(
    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
  )
)

proj4string(pointsDta) <- proj4string(shpNUTS)
res <- over(pointsDta, shpNUTS)

Update

Updates in reply to provided comments.

Shapefiles

Concering the shapefiles, I've sourced this file, that is made available through the Eurostat website.

1

It is hard to say without seeing the shapefile data but I guess that both your datasets will most likely use WGS84 as CRS. At least the X/Y coordinates use that.

So maybe if you assign

wgs84_crs = CRS("+init=epsg:4326")
proj4string(pointsDta) <- wgs84_crs

it could solve your problem.

If not, then you need to figure out what CRS the layer shpNUTS uses, assign the correct CRS and then you can use spTransform() to transform it to the same CRS as pointsDta layer.

3
  • If the point data is Eastings/Northings then the CRS is British National Grid – Phil Dec 1 '16 at 9:56
  • Well, the dataset provides both Latitude/Longitude and Eastings/Northings coordinates. I used Latitude/Longitude as it is easier to estimate the CRS for those coordinates. – Jan Caha Dec 2 '16 at 8:33
  • Thanks very much for showing interest in my post, I will try those projections today. I was actually trying to find what NUTS uses but it's not clear. I came across this page that mentions urn:ogc:def:crs:EPSG::27700 as used coordinate system. – Konrad Dec 2 '16 at 9:18

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.