I'm interested in undertaking spatial matching in R with use of the sp package on some of the publicly available data consisting of:
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:
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.
The shapefile is read via the readOGR
function available through rgdal 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)
Updates in reply to provided comments.
Shapefiles
Concering the shapefiles, I've sourced this file, that is made available through the Eurostat website.
Answer
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.
No comments:
Post a Comment