Saturday, 25 July 2015

polygon - Correctly setting up CRS systems for the NUTS (Eurostat) geographies and Doogal's data?


I'm interested in undertaking spatial matching in R with use of the 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:


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.





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)




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

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...