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

[Roger Bivand]

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

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