Wednesday, 22 February 2017

coordinate system - Why did spTransform fail here?


I'm trying to convert a map I found here into lat/lon compatible with another data set I created from Google Maps lat/lons.


From the directory with that .shp file, I run:


library(rgdal)
SG = readOGR('../data', 'MP14_PLNG_AREA_NO_SEA_PL')


A subset of the data I'm trying to co-plot is:


dim = c(5L, 2L)
# polygon coordinate matrices
poly_mats = list(structure(c(103.843, 103.843, 103.854, 103.854, 103.843,
1.318, 1.324, 1.324, 1.318, 1.318), .Dim = dim),
structure(c(103.843, 103.843, 103.854, 103.854, 103.843,
1.324, 1.329, 1.329, 1.324, 1.324), .Dim = dim),
structure(c(103.854, 103.854, 103.865, 103.865, 103.854,
1.318, 1.324, 1.324, 1.318, 1.318), .Dim = dim),

structure(c(103.854, 103.854, 103.865, 103.865, 103.854,
1.324, 1.329, 1.329, 1.324, 1.324), .Dim = dim))
# nest properly to create SP object
polySP = SpatialPolygons(lapply(seq_along(poly_mats), function(ii) {
Polygons(list(Polygon(poly_mats[[ii]])), ID = ii)
# assign lat/lon CRS
}), proj4string = CRS('+init=epsg:3857'))

Both of these objects look fine when plotted separately:


png('separate.png')

par(mfrow = c(1, 2))
plot(SG)
plot(polySP)
dev.off()

polygons separately


Not that it should matter, but for certain, the objects are currently in different CRS:


identical(proj4string(SG), proj4string(polySP))
# [1] FALSE


So we should be able to transform SG like so:


SG = spTransform(SG, proj4string(polySP))


But this fails:


png('together_no.png')
plot(SG)
plot(polySP, col = 'red', add = TRUE)
dev.off()

not together


It seems spTransform has failed, as the coordinates in the output are non-sensical as lat/lons:



coordinates(SG)[1:5, ]
# [,1] [,2]
# 0 11559649 153646.0
# 1 11569258 147405.4
# 2 11559462 150848.3
# 3 11544079 146374.9
# 4 11549925 150925.5

What's going wrong/how can this be rectified?



Answer




You've created your polygons with epsg:3857:


proj4string = CRS('+init=epsg:3857'))

which is Google's Web Mercator, not lat-long. You mean epsg:4326, WGS84 lat-long, most likely.


Some more detailed reading:



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