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

Slawomir Kozielec kozielec at gmail.com
Tue May 23 10:48:31 CEST 2017


all, Thank you for your assistance,
I resolved the issue, though I am not sure how to 'call off' my request for
advice.
Basically the MapInfo was projected to BNG (27700), so rgdal tried to
transform already projected file. On Manifold I most likely and purely by
chance set it to initial BNG and this is how it was produced correctly (i
have no practical experience with Manifold - I got it once upon a time with
intention to learn but never had time to). Probably Proj4 could have worked
as well (as you provide the initial projection), but it refuses to work
with spatial dataframes, especially lines (as points you can just extract
as coords and transform).
What I did was basically to replace by substitution (not transformation)
the proj4file attribute to BNG and then transformed it - everything went
smoothly from then on

Mike, the file format is not my choice - I am preparing an application for
my team and shp is our internal standard. Alas our GIS team leader took few
months off, some project to clean oceans off plastic, and his deputees did
not know how to arrange the file for me - so i took a 'raw' file from our
vendor. The file was in TAB format. I am not disputing whatever is better
as I am a proper layman in the order of datums, elipsums and all other
Latin words I would love to know:)


On Tue, May 23, 2017 at 8:40 AM, Michael Sumner <mdsumner at gmail.com> wrote:

> Quickly, there is no need to convert to shapefile. MapInfo TAB or MIF are
> both vastly superior to SHP and are well supported by rgdal and sf.
>
> That takes out one weak link. I'm happy to explore with Manifold if you
> provide reproducible examples. Also I would ask on georeference.org too,
> though you'll get a very disparaging (of open source) perspective from some
> of their staff. It works both ways.
>
> Cheers, Mike
>
> On Tue, 23 May 2017, 17:13 Roger Bivand <Roger.Bivand at nhh.no> wrote:
>
>> On Mon, 22 May 2017, Slawomir Kozielec wrote:
>>
>> > 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
>>
>> Please find out about ellipsoids and datums (sorry, not plural data). They
>> describe the assumed shape of the world, and the relative distance of a
>> point on the surface to an assumed 3D origin (roughly).
>>
>> In one of your CRS, the ellipsoid is Airy, in another where a=b it is a
>> sphere. In addition, EPSG 4269 presupposes NAD 1983, equivalent to GRS
>> 1980 and effectively WGS84. What you do not have are +towgs84=
>> definitions, but if these are not given, they are assumed by proj in rgdal
>> to be WGS84.
>>
>> None of these are "correct", all are based on chosen assumptions. Why
>> Manifold is assuming a spheroid is unknown. Define your datums correctly,
>> and things should clear up.
>>
>> Hope this clarifies,
>>
>> Roger
>>
>> >
>> > 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]]
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo at r-project.org
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55 <+47%2055%2095%2093%2055>; e-mail:
>> Roger.Bivand at nhh.no
>> Editor-in-Chief of The R Journal, https://journal.r-project.org/
>> index.html
>> http://orcid.org/0000-0003-2392-6140
>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list