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

Roger Bivand Roger.Bivand at nhh.no
Fri Dec 7 11:34:17 CET 2007


On Fri, 7 Dec 2007, Agustin Lobo wrote:

> 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.

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.

>
> 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]]

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:

>       X-axis translation: -87 m
>       Y-axis translation: -98 m
>       Z-axis translation: -121 m

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

>
> 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)
>

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

> 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
>> >
>> 
>
>

-- 
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