[R-sig-Geo] Problem with spTransform en +towgs84 ?
Roger Bivand
Roger.Bivand at nhh.no
Fri Sep 24 23:33:35 CEST 2010
On Fri, 24 Sep 2010, jgarcia at ija.csic.es wrote:
> Hi again,
> However,in my case, if instead of the previously mentioned version in my
> linux box, I use my windows version of R:
>
>> library(rgdal)
> Loading required package: sp
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.6.3, released 2009/11/19
> Path to GDAL shared files: C:/ARCHIV~1/R/R-29~1.2/library/rgdal/gdal
> Loaded PROJ.4 runtime: Rel. 4.6.1, 21 August 2008
> Path to PROJ.4 shared files: C:/ARCHIV~1/R/R-29~1.2/library/rgdal/proj
>>
>
> The mismatch problem is solved. So, could this be, in this case, a problem
> of the PROJ.4 Rel. 4.5.0, 22 Oct 2006, or GDAL 1.4.2.0 versions in the
> linux computer?
Older releases of the EPSG database for PROJ.4 and GDAL did not
necessarily include +towgs84= parameters for known +datum= - this has
improved with time. The burden is, though, on users to understand that
both the datum and the ellipsoid must be known, and in the case of
ambiguous "datum" definitions like ED50, what the +towgs84= should be.
See projInfo("datum") for known datum definitions.
Currently, the GDAL version is 1.7.2, and PROJ.4 is 4.7.0 (but declares
itself as 4.7.1). These are included in CRAN binaries for Windows and CRAN
extras binaries for OSX. If you install from source under Linux, these are
also what you should have, and Debian/Fedora packagers for GIS are usually
not far behind, and checking for updates in some fashion should be
possible. In R, do make a habit of running sessionInfo() and recording the
output if you need absolutely reproducible research; then
update.packages() to pick up bug fixes and enhancements. If you need to
step back, all older versions of packages are also available on CRAN.
It is known that datum specifications are a problem, which is why one
needs to be particularly careful about correct +towgs84= parameters.
Hope this helps,
Roger
>
> Thanks,
>
> Javier
> ---
>
>
>> Hi,
>> Just to add that I've also detected some problems in this sense. I've have
>> exported a shapefile from ArcGIS to KML, and also imported it into R, and
>> from there to KML. The final KML files differ in some 70m of latitude,
>> while the longitudes match.
>>
>> # import the shapefile
>>> knapdale.supgeo <- readOGR(dsn="test",layer="knapdale_supgeo")
>>> proj4string(knapdale.supgeo)
>> [1] " +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000
>> +ellps=airy +units=m +no_defs"
>>
>> # Note, the above is OSGB 1936 / British National Grid, which equals
>> # CRS("+init=epsg:27700")
>> # Then:
>>
>> knapdale.supgeoLL <- spTransform(knapdale.supgeo, CRS("+proj=longlat +
>> datum=WGS84"))
>> kmlf <- "knapdale_supgeo.kml"
>> writeOGR(knapdale.supgeoLL, dsn=kmlf, layer="supgeo", driver="KML")
>>
>> I'm not saying that the R part is responsible for the difference, but the
>> mismatch is there. Some details:
>>
>> ArcGIS 9.3 (Windows XP Pro)
>> R 2.10.1
>>> library(rgdal)
>> Loading required package: sp
>> Geospatial Data Abstraction Library extensions to R successfully loaded
>> Loaded GDAL runtime: GDAL 1.4.2.0, released 2007/06/27
>> Path to GDAL shared files: /usr/local/share/gdal
>> Loaded PROJ.4 runtime: Rel. 4.5.0, 22 Oct 2006
>> Path to PROJ.4 shared files: (autodetected)
>>
>> Best regards,
>> Javier
>>
>>
>>
>>> Roger,
>>>
>>> I apologize for not including version details etc, here they are:
>>> Loading required package: rgdal
>>> Geospatial Data Abstraction Library extensions to R successfully loaded
>>> Loaded GDAL runtime: GDAL 1.7.2, released 2010/04/23
>>> Path to GDAL shared files: /usr/share/gdal17
>>> Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
>>> Path to PROJ.4 shared files: (autodetected)
>>>
>>>> version
>>> _
>>> platform i486-pc-linux-gnu
>>> arch i486
>>> os linux-gnu
>>> system i486, linux-gnu
>>> status
>>> major 2
>>> minor 11.1
>>> year 2010
>>> month 05
>>> day 31
>>> svn rev 52157
>>> language R
>>> version.string R version 2.11.1 (2010-05-31)
>>>
>>> Also, I've put SPDF polsfl1_1000 here:
>>> https://sites.google.com/site/openfiles2/home/polsfl1_1000.rda?attredirects=0&d=1
>>>
>>> Finally, the towgs84 parameters come from the last paragraph in the pdf
>>> for The Kingdom of Spain (July) listed in the page you mention. (BTW,
>>> I recall having seen
>>> another list with slightly different parameters but cannot find it
>>> anymore)
>>> I had to put
>>> the same 3 parameters for the definition of the projection in QGIS
>>> (which qgis does not
>>> do by default, had to define a custom projection to get the datum
>>> shift working). The towgs
>>> parameters are not in the proj file of the shapefile set, at least not
>>> in the ones created
>>> by writeOGR.
>>>
>>> Hope this helps to clarify the problem. I'm still not certain that
>>> spTransform is doing anything wrong.
>>>
>>> Thanks for your help,
>>>
>>> Agus
>>>
>>>
>>> 2010/9/21 Roger Bivand <Roger.Bivand at nhh.no>:
>>>> On Tue, 21 Sep 2010, Agustin Lobo wrote:
>>>>
>>>>> Hi!
>>>>>
>>>>> Might it be that spTransform is not taking the +towgs84 string into
>>>>> account?
>>>>> I have:
>>>>>>
>>>>>> proj4string(polsfl1_1000)
>>>>>
>>>>> [1] " +proj=utm +zone=31 +ellps=intl +units=m +no_defs
>>>>> +towgs84=-87,-98,-121"
>>>>>
>>>>> and do:
>>>>>>
>>>>>> polsfl1_1000lonlat <- spTransform(polsfl1_1000, CRS("+proj=longlat
>>>>>> +datum=WGS84"))
>>>>>> writeOGR(polsfl1_1000,
>>>>>> dsn="/media/TRANSCEND/FLUXPYR/FLIGHT201009/polsfl1_1000",layer="polsfl1_1000",
>>>>>> driver="ESRI Shapefile")
>>>>>> writeOGR(polsfl1_1000lonlat,
>>>>>> dsn="/media/TRANSCEND/FLUXPYR/FLIGHT201009/polsfl1_1000lonlat.kml",layer="polsfl1_1000",
>>>>>> driver="KML")
>>>>>
>>>>> Then, in QGIS, with reprojection on the fly activated, the 2 vectors
>>>>> are not coincident.
>>>>> Instead, a polsfl1_1000lonlatv2 made from the original polsfl1_1000
>>>>> shapefile in QGIS, is coincident.
>>>>>
>>>>> As there are 2 different programs here, I'm not positive of my
>>>>> diagnostic.
>>>>> But
>>>>> I think there is a problem with spTransform
>>>>
>>>> Could you please provide links to screen shots and input and output
>>>> data
>>>> objects, as well as sessionInfo() output? Which OS is involved, and are
>>>> spTransform() and QGIS using the same PROJ.4 library and metadata. Is
>>>> your
>>>> input actually:
>>>>
>>>> CRS("+init=epsg:23031 +towgs84=-84,-107,-120")
>>>>
>>>> with the three-parameter datum shift from:
>>>>
>>>> http://www.asprs.org/resources/GRIDS/
>>>>
>>>> of July 2000? Where does yours come from? The way QGIS treats the *.prj
>>>> file
>>>> is also of interest.
>>>>
>>>> Roger
>>>>
>>>>>
>>>>> Thanks
>>>>>
>>>>> Agus
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at stat.math.ethz.ch
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>> --
>>>> Roger Bivand
>>>> Economic Geography Section, Department of Economics, Norwegian School
>>>> of
>>>> Economics and Business Administration, 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
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at stat.math.ethz.ch
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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
More information about the R-sig-Geo
mailing list