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

[Roger Bivand]

On Thu, 3 May 2012, Jon Olav Skoien wrote:

> On 02-May-12 18:39, Roger Bivand wrote:
>> On Wed, 2 May 2012, Jon Olav Skoien wrote:
>> 
>>> 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?
>> 
>> As http://spatialreference.org/ref/sr-org/6864/ points out, this is a 
>> rather broken configuration on a sphere. The definitions seem to be 
>> identical in your cases, as a @null grid or zero towgs84 coefficients will 
>> have the same effect. You could ask on the proj list, but my intuition is 
>> that this is useful for display only, not for data transfer.
>> 
>> Hope this clarifies,
>> 
>> Roger
>
> Roger,
>
> If I understand the reference correct, this projection could cause errors in 
> the order of 800 meters, while the ones we encountered are around 30 km. You 
> are probably right that we should be careful with this projection and not use 
> it without manual checking. I will also ask on the proj list and see if they 
> can give some more information.
>
> BTW, the proj4string from spatialreference.org also gives a y-transformation 
> of points defined with epsg:3857 in R:
> p1 = SpatialPoints(data.frame(x = 700000, y = 5700000), proj4string = 
> CRS("+init=epsg:3857"))
> p2 = spTransform(p1, CRS("+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"))
> p1
> p2

Yes, it does, which implies that +towgs84=0,0,0,0,0,0,0 and 
+nadgrids=@null are not equivalent in this case. Is it the minor axis 
difference in going from WGS84 to sphere that is changing y? I think that 
the proj list is the right place to follow this up, unless the final 
point on:

http://proj.maptools.org/faq.html

helps. As I read the FAQ, the definition including +towgs84= is probably 
wrong, and is forcing an understanding of +datum=WGS84. Has this changed 
between PROJ.4 releases?

Roger

>
> Thanks a lot for your help!
> Jon
>
>> 
>>> 
>>> 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]"
>>> 
>>> 
>> 
>
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no