[R-sig-Geo] what Manifold does that Rgdal fails? Projection transformation R

Slawomir Kozielec kozielec at gmail.com
Mon May 22 16:07:41 CEST 2017


Dear friends,

pretty new to the topic so let me know if i do not follow and reproducible
examples convention. I am OK with R, but spatial data is still terra
incognita

I got a Mapinfo file and i had to convert it to ESRI shp and then within r
into spatial lines data frame and use with leaflet+OS maps epsg:3857

whatever i tried within R with spTransform (rgdal) it did not work - either
was not visible or was visible with offset (approx 2 meters). Proj4 did not
want to work with SLDF.
I took the file into Manifold, converted it to 3857, returned to R,
transformed it to epsg:4269 and it works perfectly OK with leaflet + OS map
3857. But i have to automate this process to work within R, without any
Manifold intervention.

I checked how the files differed:
i selected one item (the same of course) from the file and to make it easy
to present created a centroid for it. Then I checked projection, xy, and
latlong.

the item created from the file without using Manifold:
proj4string

    "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
+y_0=-100000 +ellps=airy +units=m +no_defs"

centroid:

    528685 181836.2

centroid after spTransform to 4269

    -0.1450149 51.52028

the item from the file pre-processed by Manifold:
proj4string

    "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137
+units=m +no_defs"

centroid:

    -16321.64 6713938

centroid after spTransform to 4269:

    -0.1466198 51.52079

the latter is correct of course. Just to emphasize - it is the same file
and the same item

Few questions:

1. any idea what happens in Manifold that fails in rgdal?

2. is there any mechanism in R that I can apply to fix the error (to match
the 'with Manifold' output)? Maybe a kind of offset or work on prj files?

I tried to use the latter CRS for transformation

    P5 <-spTransform(ptp.points1P, CRS( "+proj=merc +lon_0=0 +lat_ts=0
+x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"))

and it returned still incorrect values
the xy is t least inverted but shifted

     -16142.98 6713847

the latlong is of course incorrect, but the same as without any attempt to
transform

    -0.1450149 51.52028


content of prj file not transformed - failing to correctly show

    PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy
1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

content of prj manifold transformed file - good to go

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

thanks in advance

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list