[R-sig-Geo] Problem with spTransform en +towgs84 ?
jgarcia at ija.csic.es
jgarcia at ija.csic.es
Fri Sep 24 20:59:15 CEST 2010
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?
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
>
More information about the R-sig-Geo
mailing list