Skip to content

towgs snd prj files

5 messages · Roger Bivand, Agustin Lobo

Agustin Lobo
# 
This might be a bit off-topic.
I do:
CRS arguments:
 +init=epsg:23031 +proj=utm +zone=31 +ellps=intl
+towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs
But the created prj file:
PROJCS["UTM_Zone_31_Northern_Hemisphere",GEOGCS["GCS_International
1909 (Hayford)",DATUM["D_unknown",SPHEROID["intl",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

lacks the +towgs84=-87,-98,-121,0,0,0,0 information, which is a
problem if I do not remember I have to add it for a GIS display
with reprojection on the fly.
Is there any vector format that can  be written by rgdal in which the
CRS info would include the +towgs84=-87,-98,-121,0,0,0,0 ?

Thanks!

Agus
Roger Bivand
#
On Wed, 7 Nov 2012, Agustin Lobo wrote:

            

      
      
      
    

        
There is no guarantee of any kind that GDAL/OGR will reproduce the same 
projection parameters on re-reading, and there is nothing in sp/rgdal that 
can be done about this. The thread at:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016581.html

did permit an extra argument to be passed to readGDAL() to preserve datum 
in addition to equivalent ellps and towgs84 tags.

The underlying difficulties are that none of the various spatial reference 
systems are unequivocal, and changing between PROJ.4, WKT, and others 
(EPSG codes in different versions) may drop tags, especially where they 
are not uniquely defined.

On R-forge (revision 383) I have added a morphToESRI= logical argument to 
writeOGR(), NULL default, set automatically TRUE if the driver is "ESRI 
Shapefile" and FALSE for all other drivers (this is the same as internal 
code has always been). The user may try to see whether this makes a 
difference, I didn't see any difference in the exported WKT written in the 
*.prj file.

The main issue is that for any given specification, the towgs84 tag may 
not be unique (or other tags), so that the PROJ.4/GDAL policy is not to do 
anything that might not be accurate in conversion. There are no easy 
solutions here.

In your case, for the copy og EPSG distributed with rgdal in the binary 
Windows and OSX packages, and for my PROJ.4/GDAL combination:

            

      
      
      
    

        
CRS arguments:
  +init=epsg:23031 +proj=utm +zone=31 +ellps=intl
+towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs

but http://spatialreference.org/ref/epsg/23031/ does not include a 
+towgs84 term. Which is credible? This is an ED50 projection/datum, so I'd 
actually trust dropping towgs84 more than keeping it. The web page 
suggests that the ED50 datum should be defined by country with respect to 
towgs84; I often look at Grids & Datums to cross-check:

http://www.asprs.org/Grids-Datums.html

Hope this helps,

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
Just to complement the info I provided, QGIS writes, for the same
proj4 info, the following prj file:
PROJCS["ED50_UTM_zone_31N",GEOGCS["GCS_European_1950",DATUM["D_European_1950",SPHEROID["International_1924",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

which is a bit different than the prj written by writeOGR() (see former message)

The point is that when QGIS opens a shapefile with the prj written by
writeOGR(), the towgs84 coefficients are not assumed, while
when it reads a shapefile with the prj written by QGIS (above), the
towgs84 coefficients are assumed. And this has implications
in the onthefly reprojection.

I do not know if this is a QGIS idiosyncratic behaviour or it is
generally assumed that ["ED50_UTM_zone_31N" includes the towgs84
coefficients

Agus

        
On Wed, Nov 7, 2012 at 11:44 AM, Agustin Lobo <alobolistas at gmail.com> wrote:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
ok, thanks, just got your message while I was typing my second one.
It seems we have to live with this issue. The best is probably make a version
of the object in both datums via spTransform() and then write 2
shapefiles as well.

Or perhaps using spatialite objects instead of shapfiles would be a solution?
http://www.gaia-gis.it/spatialite-2.0/SpatiaLite2-tutorial.html#t5

Agus

        
On Wed, Nov 7, 2012 at 1:21 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 7 Nov 2012, Agustin Lobo wrote:

        

            

      
      
      
    

        
The issue is that where the +towgs84 coefficients are not known with 
complete certainty (the case for ED50, where they differ a good deal 
depending on where you are), they are dropped when converting from PROJ.4 
to the internal representation used by OGR:

library(rgdal)
pt <- matrix(c(10000, 3000000), ncol=2)
spdf <- SpatialPointsDataFrame(SpatialPoints(pt), data=data.frame(ID=1))
proj4string(spdf) <- CRS("+init=epsg:23031")
proj4string(spdf)
td <- tempdir()
writeOGR(spdf, td, "tmp", driver="ESRI Shapefile") 
readLines(paste(td, "tmp.prj", sep="/"))
proj4string(readOGR(td, "tmp"))

If we use OSGB instead, with D_OSGB_1936 rather than D_unknown for its 
datum, the +towgs84 values are preserved:

spdf <- SpatialPointsDataFrame(SpatialPoints(pt), data=data.frame(ID=1))
proj4string(spdf) <- CRS("+init=epsg:27700")
proj4string(spdf)
writeOGR(spdf, td, "tmposgb", driver="ESRI Shapefile")
readLines(paste(td, "tmposgb.prj", sep="/"))
proj4string(readOGR(td, "tmposgb"))

So the logic in OGR is trying to ensure that only information that is 
authoritative is encoded in the output file. QGIS has a "custom" way of 
handling spatial reference, trying to make things "easy" for the user, 
whereas OGR prefers not to guess when it doesn't have authority. ED50 is 
not a standard at all, it is a collection of national (sometimes 
sub-national) standards that have also changed over time. It is only since 
GRS80/WGS84 that we have satellite-based measurements on which to base 
spheroid and datum definitions, hence the importance of having authority 
for transformations to and from WGS84. It doesn't make matters easier that 
the Earth also changes, but these changes are small in comparison with the 
differences induced by poor datum transformation. The datum values 
recognised by PROJ.4 are:

projInfo("datum")

which includes OSGB36, but not ED50. This:

http://www.ogp.org.uk/pubs/373-10.pdf

gives a flavour of the difficulties involved in oil exploration in the 
North Sea when also transforming between ED50 and WGS84 - note that these 
coefficients are very different from those shown in another source:

http://earth-info.nga.mil/GandG/coordsys/onlinedatum/CountryEuropeTable.html

referred to in this thread:

http://lists.osgeo.org/pipermail/gdal-dev/2009-May/020656.html

This is a never-ending story ...

Roger

            

      
      
      
    

        
No, nothing to do with this.