[R-sig-Geo] Projection problem (epsg 3857 - maybe not R-related)

[Jon Olav Skoien]

Hi,

We are having a problem with projections that we dont know the source of 
(bugs in software or lack of understanding). The problem might not be 
related to R, but I hope someone can still help. We have a set of data 
(polygons and rasters) that are all supposed to be in projection 
epsg:3857 (WGS 84 / Pseudo-Mercator, often referred to as web mercator 
as it is used in web services such as Google Maps). However, they do not 
appear to be the same when using them in R. Below is an example that 
starts with a hypothetical polygon somewhere in Europe, and then 
transforming it with the three different proj4strings we got from the 
different objects :

cords = matrix(c( 650000, 700000, 700000, 650000, 650000, 5600000, 
5600000, 5700000, 5700000, 5600000), ncol = 2)
pol = 
SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(cords)), 
ID = "1")),
     proj4string = CRS("+init=epsg:3857")), data = data.frame(test = 1))
proj4string(pol)
# [1] " +init=epsg:3857 +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"
# This is the proj4string that R creates from this epsg code

# Below are the proj4strings that we got from the imported objects
# 1 is from an export and transformation from a postgis database of 
polygons (latlong) with ogr2ogr and target crs epsg:3857
# 2 is an export from the same data base with psql2shape
# the same proj4strings were also generated from the .prj-files when 
transforming with QGIS.
# 3 is from a geoTIFF, and practically identical to the proj4string that 
R generated from epsg:3857
projs = list()
projs[[1]] = " +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 
+units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
projs[[2]] = " +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 
+b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
projs[[3]] = " +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"

# Here we print the y-coordinate of the original polygon (centroid) 
together with the y-coordinate
# of the reprojected polygon. In addition we write the result to a 
shapefile:
pols = list()
for (i in 1:3) {
   pols[[i]] = spTransform(pol, CRS(projs[[i]]))
#  writeOGR(pols[[i]], "e:/temp/Rtmp/polygons", paste("pol", i, sep = 
""), driver = "ESRI Shapefile")
   print(paste(i, coordinates(pol)[,2], coordinates(pols[[i]])[,2]))
}
#[1] "1 5650000 5619680.01465088"
#[1] "2 5650000 5619663.0232111"
#[1] "3 5650000 5650000"

The writeOGR call above is to also be able to examine the .prj files. 
The last two of these are:
pol2.prj
PROJCS["Mercator",GEOGCS["GCS_unnamed ellipse",DATUM["D_unknown", 
SPHEROID["Unknown",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]], 
PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0], 
PARAMETER["false_northing",0],PARAMETER["standard_parallel_1",0.0],UNIT["Meter",1]] 


pol3.prj
PROJCS["Mercator_2SP",GEOGCS["GCS_unnamed ellipse",DATUM["D_unknown",  
SPHEROID["Unknown",6378137,0]],PRIMEM["Greenwich",0], 
UNIT["Degree",0.017453292519943295]], 
PROJECTION["Mercator_2SP"],PARAMETER["standard_parallel_1",0], 
PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

The only differences between these two are the PROJCS and PROJECTION 
parameters, which are "Mercator" and "Mercator_2SP", respectively, and 
the order of some parameters; the last one puts 
PARAMETER["standard_parallel_1",0] before meridian and false easting and 
northing, whereas the first one puts it after. The parameters are 
otherwise the same.

The questions are therefore:
Are the differences due to the use of "Mercator" and "Mercator_2SP"?
Can it be correct that we get different coordinates for the same epsg-code?
If not, which version should we trust, or how can we make sure that the 
same epsg-code refers to the same projection?

Thanks,
Jon

BTW:
 > sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United 
States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgdal_0.7-8 sp_0.9-98

loaded via a namespace (and not attached):
[1] grid_2.15.0    lattice_0.20-6
 > getPROJ4VersionInfo()
[1] "Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]"

-- 
Jon Olav Skøien
Joint Research Centre - European Commission
Institute for Environment and Sustainability (IES)
Land Resource Management Unit

Via Fermi 2749, TP 440,  I-21027 Ispra (VA), ITALY

jon.skoien at jrc.ec.europa.eu
Tel:  +39 0332 789206

Disclaimer: Views expressed in this email are those of the individual and do not necessarily represent official views of the European Commission.