Skip to content

spTransform: datum syntax?

4 messages · Roger Bivand, Agustin Lobo

Original
Agustin Lobo
# 
Dear list,

I have to transform the datum of an spatial object.
First, I create the object with:
pressp <- SpatialPointsDataFrame(pres[,c(10,6)], data=pres[,c(1,2,11)],
+                proj4string = CRS("+proj=longlat +ellps=WGS84"))

Then I have to transform to UTM but using the International or Hayford 
datum.

The following works:
presspUTM <- spTransform(pressp, CRS("+proj=utm +zone=31))

but then I keep WGS84.

The following also works:
presspUTM <- spTransform(pressp, CRS("+proj=utm +zone=31 +datum=NAD27"))

but this is not the datum I want

Then I've tried:
presspUTM <- spTransform(pressp, CRS("+proj=utm +zone=31 +datum=ED50"))
presspUTM <- spTransform(pressp, CRS("+proj=utm +zone=31 
+datum=International"))
presspUTM <- spTransform(pressp, CRS("+proj=utm +zone=31 +datum=Int"))
presspUTM <- spTransform(pressp, CRS("+proj=utm +zone=31 +datum=Hayford"))

and many other similar names, also tried with ellps=

The doc refers to the proj.4 doc, but cannot find any list of datums 
there, including the pdf docs of proj.4

Does anyone know where I could find a list of supported datums in proj.4
and their correct names for the +datum= argument?

Thanks!

Agus
Roger Bivand
#
On Thu, 6 Dec 2007, Agustin Lobo wrote:

            

      
      
      
    

        
projInfo("datum")

but the problem is to find the actual datum for the ellipse:

projInfo("ellps")

you need. You can move to "+ellps=intl", but finding the three or seven 
parameter +towgs84= values is more difficult. The Grids & Datums column 
at: http://www.asprs.org/resources/grids/ is often a good place to start.

You can search in the output data frame from:

EPSG <- make_EPSG()

too, but you'd need to know more about the target coordinate reference 
system.

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
I'm surprised that are there only 10 available datums. I understand
that this is a problem of proj4, not R. But if R has OGR-based
functionality, there should be somewhere a list with more datums.

What I'm doing now is importing other georeferenced layers of the same 
area to copy the proj4string from them (I remind you that I'm making
a sp polygon file out of lat-lon vertex coordinates and
then reprojecting to UTM). But I'm surprised by the following:

For a shp file of the same region, the proj file is:
PROJCS["ED_1950_UTM_Zone_31N",GEOGCS["GCS_European_1950",DATUM["D_European_1950",SPHEROID["International_1924",6378388.0,297.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",3.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

But once imported, I get only:
 > delme <- 
readOGR("C:/ALOBO/dipu2006/espaisoberts_P037/FBARO_OUTPUTS/STLLOR","STLL_EOA_FB3")
OGR data source with driver: ESRI Shapefile
Source: "C:/ALOBO/dipu2006/espaisoberts_P037/FBARO_OUTPUTS/STLLOR", 
layer: "STLL_EOA_FB3"
with  144  rows and  6  columns
 > delme at proj4string
An object of class ?CRS?
Slot "projargs":
[1] " +proj=utm +zone=31 +ellps=intl +units=m +no_defs"

Instead, I have another sp object of the same region (but do not
find any notes on how I managed to import it) with a much richer pro4string:

 > limitGFO at proj4string
CRS arguments:
  +proj=tmerc +lat_0=0 +lon_0=2.999999982811267 +k=0.999600 +x_0=500000 
+y_0=0
+ellps=intl +units=m +no_defs

Is there any way of making a correct proj4string out of a *.prj file?
(or from a *.txg file for raster layers)

For completeness, I also provide the following information:

If I read the gdal info of a geotif file of the same area:
 > GDALinfo("C:/ALOBO/dipu2006/foc2003/fcov2005.tif")
Closing GDAL dataset handle 0x021bab40...  destroyed ... done.
rows        1950
columns     2221
bands       1
ll.x        404629.1
ll.y        4624285
res.x       10
res.y       10
oblique.x   0
oblique.y   0
driver      GTiff
projection  +proj=utm +zone=31 +ellps=intl +units=m +no_defs
file        C:/ALOBO/dipu2006/foc2003/fcov2005.tif

While its companion txg file is much more comprehensive (not sure if the
actual info stored in the tif file itself is as comprehensive, don't 
have any way of looking at the geotiff header directly):

Format: GeoTIFF
Description: DELTA proporcional a PVI05
Unknown text 'impexpnames.GeoTIFF'
Image type: 32-bit floating-point
Image size: 2221 columns by 1950 lines
Resolution: 10.000000 (column) by 10.000000 (line) meters per cell
Units: meters
Corners:
    Upper Left:  404629.090000 E  4624285.300000 N
   Upper Right:  426839.090000 E  4624285.300000 N
    Lower Left:  404629.090000 E  4604785.300000 N
   Lower Right:  426839.090000 E  4604785.300000 N

Coordinate Reference System:
Name: ED50 / UTM zone 31N (CM 3E)
Projection: UTM zone 31N (CM 3E)
    Method: Transverse Mercator
    Latitude of natural origin: N 0 00 00.000
    Longitude of natural origin: E 3 00 00.000
    Scale factor at natural origin: 0.9996
    False easting: 500000 m
    False northing: 0 m
Datum: European Datum 1950 (ED50)
    Type: Geodetic
    Epoch: 1950
    Ellipsoid: International 1924
       Semi-major axis 6378388
       Inverse flattening: 297
    Prime Meridian: Greenwich
    Valid Area
       Europe - ED50
Coordinate System
    Type: Projected
    Axis 1: Easting
       Unit: meter
       Symbol: E
    Axis 2: Northing
       Unit: meter
       Symbol: N
Datum Transformations
    To: World Geodetic System 1984 (WGS84)
       Method: Geocentric translation
       X-axis translation: -87 m
       Y-axis translation: -98 m
       Z-axis translation: -121 m
       Valid Area
          Europe - west

   Note: Georeference control points are at outside edges of image.
         This is the upper left corner of the upper left cell,
         the upper right corner of the upper right cell, etc.






Roger Bivand escribi?:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 7 Dec 2007, Agustin Lobo wrote:

        

            

      
      
      
    

        
Well, that is why I suggested looking at the excellent Grids & Datums 
columns. Their author is a highly experienced cartographer and surveyor, 
and each column (monthly for many years) demonstrates that the only common 
factor in coordinate reference systems is surprise.

Briefly, your original states that the datum is ED50. However, ED50 only 
is a framework that differs not only from country to country, but often 
also within countries. It is only from WGS84 that surveyors have had 
access to a satellite measurement based "global standard" to link the 
"local standards" to.

Consequently EPSG and PROJ.4 do not "assume" a single ED50, which may be 
correct in some places but not in others - the datums used are those that 
are "know" to be invariant.

            

      
      
      
    

        
PROJCS["ED_1950_UTM_Zone_31N",GEOGCS["GCS_European_1950",
DATUM["D_European_1950",SPHEROID["International_1924",6378388.0,297.0]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",3.0],
PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]

If you can identify the EPSG code number for your specific zone 31 
setting, you can get maybe to the +towgs84= argument. Since the datum is 
not known precisely, you have to provide a conversion to WGS84, with 
either a table look-up (North America NAD27->NAD83), or three or seven 
parameter sets defining the datum shift.

The July 2000 G&D is for the Kingdom of Spain, and its last paragraph 
gives an approximate three parameter conversion from ED50 to WGS84. If you 
are willing to get into fiction, "La carta esf?rica" ("The Nautical Chart" 
I believe) by Arturo P?rez-Reverte or translations may amuse.

The string you need to convert to may be:

CRS("+proj=utm +zone=31 +ellps=intl +towgs84=-84,-107,-120,0,0,0,0")

or with your data source:

            

      
      
      
    

        
CRS("+proj=utm +zone=31 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0")

but I'd likely trust G&D here unless ground truth (a known point in both 
reference systems) said otherwise. If my signs or values are wrong, 
someone else will get to the treasure first.

Roger

            

      
      
      
    

        
PS: these are exactly equivalent, utm is a tmerc with specified central 
longitude.