[R-sig-Geo] spTransform: datum syntax?

Agustin Lobo Agustin.Lobo at ija.csic.es
Fri Dec 7 10:09:43 CET 2007


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ó:
> On Thu, 6 Dec 2007, Agustin Lobo wrote:
> 
>> 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
> 
> 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
> 
>>
>> 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
>>
> 

-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster




More information about the R-sig-Geo mailing list