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

[Jon Olav Skoien]

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

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


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