1 Jun 2010 19:42

### Spheroidal gnomonic projection

```The gnomonic map projection is a central projection of the sphere on a
tangent plane.  The key property of the projection is that all geodesics
map to straight lines in the projection.

This property cannot be preserved for an ellipsoid.  However, we can
obtain a projection where the projected geodesics are approximately
straight close to the center of the projection.  It is obtained as the
limit of a 2-point azimuthal projection as the two points approach one
another.  The method generalizes to any 2-dimensional surface; but I
only have the necessary formulas worked out for the ellipsoid.

This projection has the following properties:

(1) azimuthal, all azimuths from the center are correct;
(2) hence, all lines through the center are geodesics;
(3) all other lines are *approximately* geodesics.

To quantify point 3, consider this projection with some specified center
on the WGS84 ellipsoid.  Take an arbitrary pair of points within r =
1000 km of the center.  Draw a straight line between these on the map;
then back project this line onto the ellipsoid.  How close is this line
to a geodesic?  I find that *at worst*,

max error in initial/final azimuth = 1.0" = 0.39 * f*(r/a)^3
max deviation from geodesic = 1.66 m = 0.13 * f*(r/a)^3*r

(a = major radius, f = flattening, and the formulas show the scaling of
the errors).  For comparison, for a gnomonic projection obtained by
projecting from the center of the ellipsoid (i.e., approximating the
geodesic by a great ellipse), the equivalent figures are about 100 times
```

1 Jun 2010 20:18

### Re: Spheroidal gnomonic projection

Cole did the two-point azimuthal equidistant on the Helmert ellipsoid back in the early 1900's for the Survey of Egypt.  I do not have a copy, though.  It's pretty obscure, but the intent was to provide a graphic for determining the direction to face for daily prayers.

Clifford J. Mugnier, C.P., C.M.S.
Chief of Geodesy,
Center for GeoInformatics
Department of Civil Engineering
Patrick F. Taylor Hall 3223A
LOUISIANA STATE UNIVERSITY
Baton Rouge, LA  70803
Voice and Facsimile:  (225) 578-8536 [Academic]
Voice and Facsimile:  (225) 578-4578 [Research]
Cell: (225) 238-8975 [Academic & Research]
Honorary Life Member of the
Louisiana Society of Professional Surveyors
Fellow Emeritus of the ASPRS
Member of the Americas Petroleum Survey Group

From: proj-bounces <at> lists.maptools.org on behalf of Charles Karney
Sent: Tue 01-Jun-10 12:42
To: proj <at> lists.maptools.org
Cc: deakin <at> rmit.edu.au; william <at> karney.com; kevin <at> karney.com; knud.poder <at> mail.dk; Craig.M.Rollins <at> nga.mil
Subject: [Proj] Spheroidal gnomonic projection

The gnomonic map projection is a central projection of the sphere on a
tangent plane.  The key property of the projection is that all geodesics
map to straight lines in the projection.

This property cannot be preserved for an ellipsoid.  However, we can
obtain a projection where the projected geodesics are approximately
straight close to the center of the projection.  It is obtained as the
limit of a 2-point azimuthal projection as the two points approach one
another.  The method generalizes to any 2-dimensional surface; but I
only have the necessary formulas worked out for the ellipsoid.

This projection has the following properties:

(1) azimuthal, all azimuths from the center are correct;
(2) hence, all lines through the center are geodesics;
(3) all other lines are *approximately* geodesics.

To quantify point 3, consider this projection with some specified center
on the WGS84 ellipsoid.  Take an arbitrary pair of points within r =
1000 km of the center.  Draw a straight line between these on the map;
then back project this line onto the ellipsoid.  How close is this line
to a geodesic?  I find that *at worst*,

max error in initial/final azimuth = 1.0" = 0.39 * f*(r/a)^3
max deviation from geodesic = 1.66 m = 0.13 * f*(r/a)^3*r

(a = major radius, f = flattening, and the formulas show the scaling of
the errors).  For comparison, for a gnomonic projection obtained by
projecting from the center of the ellipsoid (i.e., approximating the
geodesic by a great ellipse), the equivalent figures are about 100 times
worse:

max error in initial/final azimuth = 108" = 1.0 * f*(r/a)
max deviation from geodesic = 265 m = 0.5 * f*(r/a)*r

The method of constructing the projection entails defining two
quantities,

* m, the reduced length.  Take two geodesics of length s which start
from the center with azimuths differing by dalpha.  The endpoints
are separated by m * dalpha.

* M, the geodesic scale.  Take two geodesics of length s which are
parallel close to the center and initially separated by dt.  The
endpoints are separated by M * dt.

The reduced length was introduced by Gauss (1827) and Christoffel
(1868).  Helmert (1880), Sec. 6.5, gives an explicit formula for the
ellipsoid.  The geodesic scale was introduced (I think) by

G. V. Bagratuni,
Course in spheroidal geodesy,
FTD-MT-64-390 (US Air Force, 1967), Sec. 17
[Kurs sferoidicheskoi geodezii, (Moscow, 1962)].

(Bagratuni uses the symbol n.)  An explicit formula for the ellipsoid is
given in version 1.2 of GeographicLib by GeodesicLine::Scale.

To carry out the projection for an arbitrary point, solve the inverse
geodesic problem between the center and the point.  Let alpha be the
azimuth of the geodesic at the center and s be its length.  Compute

rho = m/M

Then the projected point is

x = rho * sin(alpha); y = rho * cos(alpha)

The radial scale is 1/M^2 and the transverse scale is 1/M.  For a
sphere, m = a * sin(s/a) and M = cos(s/a), thus

rho = a * tan(s/a)

which gives the spherical gnomonic projection.

In order to reverse the projection, compute

alpha = atan2(x, y); rho = sqrt(x^2 + y^2)

Use Newton's method to solve for s given rho; for this we need drho/ds
which is given by the radial scale = 1/M^2.  Solve the direct geodesic
problem from the center using length s and azimuth alpha to recover the
original point.

I just checked in a Gnomonic class for GeographicLib that implements the
forward and reverse projections.  It's available through svn from
SourceForge.

This projection minimizes the error in the straightness of geodesics
near the center.  If, instead, you want to minimize the error over some
region, then use a two-point azimuthal projection with two base points
judiciously chosen within the region.

Questions:

(1) Is this new?
(2) Is Bagratuni this first person to define M?
(3) Is GeographicLib the first place where the formula for M for an
ellipsoid is given?
(4) Does anyone have a better copy of Bagratuni than the one given
above.  (I'll even take a good copy in the Russian original.)

--
Charles Karney <ckarney <at> sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662
_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

```_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj```
1 Jun 2010 23:58

### Re: Spheroidal gnomonic projection

```> From: Clifford J Mugnier [cjmce <at> lsu.edu]
> Sent: Tuesday, June 01, 2010 14:18
>
> Cole did the two-point azimuthal equidistant on the Helmert ellipsoid
> back in the early 1900's for the Survey of Egypt.  I do not have a
> copy, though.  It's pretty obscure, but the intent was to provide a
> graphic for determining the direction to face for daily prayers.

Thanks, this gives me something to track down.  I assume that "two-point
azimuthal equidistant" means the azimuth from basepoint A is correct and
the distance from basepoint B is correct?  I had thought that Close (and
not Cole) was the originator of this projection:

C. F. Close, 1935, Two-Point Azimuthal-Equidistant projection,
Geographical Journal, v. 86, no. 5, p. 445-446.

But maybe Cole generalizes this the the ellipsoid?

The gnomonic projection is the limit of a "two-point azimuthal"
projection where the two azimuths from A and B are correct, and this was
given (for the sphere only?) by Maurer, Immler, and Close (again).  All
of these are in papers I haven't been able to get hold of yet.

--Charles
--
Charles Karney <ckarney <at> sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662
_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

```
7 Jun 2010 09:38

### TM06 for Portugal

```On Thu, 2010-02-18 at 16:40 +0100, Roger Oberholtzer wrote:

I now have some test data for the conversion from WGS84 lat/long to TM06
northing/eastings. The suggested proj conversion of:

"+proj=tmerc "
"+ellps=GRS80 "
"+k=1.0 "
"+x_0=0 "
"+y_0=0 "
"+lon_0=-8.133108333333334 "
"+lat_0=39.66825833333333 "
"+towgs84=0,0,0,0,0,0,0 "

seems in fact to work with ETRS89 lat/long, not WGS84 lat/long. When
provided with WGS84 lat/long, the resulting northings/eastings are a few
meters away. (Assuming my test data from the Portuguese Road Department
is correct. I have double checked, so it might very well be.)  In fact,
it is a small error, but I want to provide the best I can. I am guessing
that I need to change my +towgs84= to cause the WGS84 lat/longs to be
converted to ETRS89. Does anyone know the proper values?

This begs the question: what does
http://spatialreference.org/ref/epsg/3763/ assume is the datum of the
original lat/longs? Seems it assumes ETRS89, not WGS84.

Here is a sample location that I have been provided (all altitudes set
to 0.0):

Ref (WGS84):    Latitude:   41.6743830556 Longitude:   -8.3312663889
Ref (ETRS98):    Latitude:   41.6744052222 Longitude:   -8.3312324722
Ref (TM06):    Easting:    -16498.151 Northing:       222796.720

The proj conversion above properly converts the ETRS98 locations to
TM06. But I want to convert the WGS84 locations to TM06. When I try
that, I get:

Easting: -16500.981 Northing: 222794.265

which is 3.747 meters away.

--

--
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj```
7 Jun 2010 11:18

### Re: TM06 for Portugal

Roger

The test data cannot be correct. ETRS89 coincides with WGS84 (or better ITRF) in 1989 and the current discplacement is only a few decimeters. The lat-long differences between your coordinates are 0.080" and 0.122", which correspond to approximately 2.9 m in lat and 2.5 m in long, i. e the differences you get in the projected coordinates.

Yout test point belongs to the portuguese geodetic network. The ETRS89 coordinates are correct, as you can check in:

http://www.igeo.pt/produtos/geodesia/vg/rgn/ficha.asp?num=873
(open with Internet Explorer)

Where did you get the WGS84 coordinates? They must be wrong.

Regards

José A. Gonçalves
Univ. of Porto, Portugal

2010/6/7 Roger Oberholtzer
On Thu, 2010-02-18 at 16:40 +0100, Roger Oberholtzer wrote:

I now have some test data for the conversion from WGS84 lat/long to TM06
northing/eastings. The suggested proj conversion of:

"+proj=tmerc "
"+ellps=GRS80 "
"+k=1.0 "
"+x_0=0 "
"+y_0=0 "
"+lon_0=-8.133108333333334 "
"+lat_0=39.66825833333333 "
"+towgs84=0,0,0,0,0,0,0 "

seems in fact to work with ETRS89 lat/long, not WGS84 lat/long. When
provided with WGS84 lat/long, the resulting northings/eastings are a few
meters away. (Assuming my test data from the Portuguese Road Department
is correct. I have double checked, so it might very well be.)  In fact,
it is a small error, but I want to provide the best I can. I am guessing
that I need to change my +towgs84= to cause the WGS84 lat/longs to be
converted to ETRS89. Does anyone know the proper values?

This begs the question: what does
http://spatialreference.org/ref/epsg/3763/ assume is the datum of the
original lat/longs? Seems it assumes ETRS89, not WGS84.

Here is a sample location that I have been provided (all altitudes set
to 0.0):

Ref (WGS84):    Latitude:   41.6743830556 Longitude:   -8.3312663889
Ref (ETRS98):    Latitude:   41.6744052222 Longitude:   -8.3312324722
Ref (TM06):    Easting:    -16498.151 Northing:       222796.720

The proj conversion above properly converts the ETRS98 locations to
TM06. But I want to convert the WGS84 locations to TM06. When I try
that, I get:

Easting: -16500.981 Northing: 222794.265

which is 3.747 meters away.

--
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

```_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj```
7 Jun 2010 13:38

### Re: TM06 for Portugal

```On Mon, 2010-06-07 at 10:18 +0100, Jose Gonçalves wrote:
> Roger
>
> The test data cannot be correct. ETRS89 coincides with WGS84 (or
> better ITRF) in 1989 and the current discplacement is only a few
> decimeters. The lat-long differences between your coordinates are
> 0.080" and 0.122", which correspond to approximately 2.9 m in lat and
> 2.5 m in long, i. e the differences you get in the projected
> coordinates.
>
> Yout test point belongs to the portuguese geodetic network. The ETRS89
> coordinates are correct, as you can check in:
>
> http://www.igeo.pt/produtos/geodesia/vg/rgn/ficha.asp?num=873
> (open with Internet Explorer)
>
> Where did you get the WGS84 coordinates? They must be wrong.

(meaning that given the proj transformation I described, the expected
northing/easting values are derived). I was surprised when they provided
two sets of lat/long values. They explained them by saying that the two
sets of LAT/LONG describe the same physical location, one in WGS84 and
the other in ETRS89. I have to trust that they have not mislabeled the
data. Especially after I asked for a confirmation.

What I am trying to accomplish is transforming WGS84 LAT/LONG from a GPS
receiver to both TM06 and MPG73 northings/eastings. I have started with
TM06. I will not have any ETRS89 lat/long values. The only reason they
might exist would be if proj makes them internally as a step in
converting the source WGS84 lat/longs into the desired TM06
northing/eastings.

How would you convert WGS84 lat/long into TM06 northing/easting? I know
that WGS84 is so close to ETR89 here that one 'could' use the WGS84
values in the transform from ETRS89 to TM06. But how would I do it the
'proper' way with proj?

>
> Regards
>
> José A. Gonçalves
> Univ. of Porto, Portugal
>
>
>
> 2010/6/7 Roger Oberholtzer <roger <at> opq.se>
>         On Thu, 2010-02-18 at 16:40 +0100, Roger Oberholtzer wrote:
>
>         I now have some test data for the conversion from WGS84
>         lat/long to TM06
>         northing/eastings. The suggested proj conversion of:
>
>              "+proj=tmerc "
>              "+ellps=GRS80 "
>              "+k=1.0 "
>              "+x_0=0 "
>              "+y_0=0 "
>              "+lon_0=-8.133108333333334 "
>              "+lat_0=39.66825833333333 "
>              "+towgs84=0,0,0,0,0,0,0 "
>
>         seems in fact to work with ETRS89 lat/long, not WGS84
>         lat/long. When
>         provided with WGS84 lat/long, the resulting northings/eastings
>         are a few
>         meters away. (Assuming my test data from the Portuguese Road
>         Department
>         is correct. I have double checked, so it might very well be.)
>          In fact,
>         it is a small error, but I want to provide the best I can. I
>         am guessing
>         that I need to change my +towgs84= to cause the WGS84
>         lat/longs to be
>         converted to ETRS89. Does anyone know the proper values?
>
>         This begs the question: what does
>         http://spatialreference.org/ref/epsg/3763/ assume is the datum
>         of the
>         original lat/longs? Seems it assumes ETRS89, not WGS84.
>
>         Here is a sample location that I have been provided (all
>         altitudes set
>         to 0.0):
>
>          Ref (WGS84):    Latitude:   41.6743830556 Longitude:
>         -8.3312663889
>         Ref (ETRS98):    Latitude:   41.6744052222 Longitude:
>         -8.3312324722
>          Ref (TM06):    Easting:    -16498.151 Northing:
>         222796.720
>
>         The proj conversion above properly converts the ETRS98
>         locations to
>         TM06. But I want to convert the WGS84 locations to TM06. When
>         I try
>         that, I get:
>
>          Easting: -16500.981 Northing: 222794.265
>
>         which is 3.747 meters away.
>
>         --
>         Roger Oberholtzer
>
>         OPQ Systems / Ramböll RST
>
>         Ramböll Sverige AB
>         Krukmakargatan 21
>         P.O. Box 17009
>         SE-104 62 Stockholm, Sweden
>
>         Office: Int +46 10-615 60 20
>         Mobile: Int +46 70-815 1696
>
>         _______________________________________________
>         Proj mailing list
>         Proj <at> lists.maptools.org
>         http://lists.maptools.org/mailman/listinfo/proj
>
> _______________________________________________
> Proj mailing list
> Proj <at> lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj

--

--
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj```
7 Jun 2010 17:04

### Re: rotated pole ob_tran help needed!

```
Haileye wrote

> Hello All,
>
> I am trying to convert rotated pole coordinates to regular lat & lon coordinates. I have both coordinates,
but am trying to get the right proj4 definition (or at least close). The rotated pole projection
coordinates are defined as yc & xc, the projection description is described below:
>         double yc(yc) ;
>                 yc:axis = "Y" ;
>                 yc:standard_name = "grid_latitude" ;
>                 yc:long_name = "latitude in rotated grid" ;
>                 yc:units = "degrees" ;
>         double xc(xc) ;
>                 xc:axis = "X" ;
>                 xc:standard_name = "grid_longitude" ;
>                 xc:long_name = "longitude in rotated pole grid" ;
>                 xc:units = "degrees" ;
>         char rotated_pole ;
>                 rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
>                 rotated_pole:grid_north_pole_latitude = 42.5 ;
>                 rotated_pole:grid_north_pole_longitude = 83. ;

Hi, I am quite new to proj and trying to use it probably for the same
purpose as indicated here (processing of Numerical Weather Prediction
model data on rotated grids), so I come back to this old thread. This is
the link to the archive, in case it has been forgotten:

http://osgeo-org.1803224.n2.nabble.com/rotated-pole-ob-tran-help-needed-td4257959.html

I can reproduce haileye's result with the following command line
parameters to proj:

invproj -f "%f" -m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lon_p=83
+o_lat_p=42.5 +lon_0=180
in> -127.163970947 12.5568675995
out> -33.880011	-28.379995

then the opposite transformation:

proj -f "%f" -m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lon_p=83
+o_lat_p=42.5 +lon_0=180
in> -33.880011	-28.379995
out> -127.163971	12.556868

-m 57.295779506 is needed for getting output or providing input in
+lon_0=180 is required otherwise there is a shift of 180 degrees,
probably because the transformation is underdetermined with only the
given parameters?!

From what I could understand, the point here is that the interpretation
of pole coordinates is the opposite in the GIS/Cartography community and
in the Numerical Weather Prediction community: in the first case they
are the coordinates of the Earth NP in the map system, while in the
second case they are the geographical coordinates of the rotated NP. For
this reason the transformations are used in the opposite way as designed.
It would be nice to hear the opinion from some proj guru on whether this
is a correct use, or a different approach should be used. I enclose the
definition of "rotated grid" from the World Meteorological Organization"
document, here Southern Pole is used and point (c), angle of rotation
should be 0 in the examples examined here.

=== WMO specification ===
(6) Three parameters define a general latitude/longitude coordinate
system, formed by a general rotation of the sphere. One
choice for these parameters is:
(a)  The geographic latitude in degrees of the southern pole of the
coordinate system, θp for example;

(b)  The geographic longitude in degrees of the southern pole of
the coordinate system, λp for example;

(c)  The angle of rotation in degrees about the new polar axis
(measured clockwise when looking from the southern to
the northern pole) of the coordinate system, assuming the new
axis to have been obtained by first rotating the
sphere through λp degrees about the geographic polar axis, and
then rotating through (90 + θp) degrees so that
the southern pole moved along the (previously rotated)
Greenwich meridian.
=== end WMO specification ===

Thank you very much for your attention, bye, Davide

--

--
============================= Davide Cesari ============================
Servizio IdroMeteorologico ARPA Emilia Romagna
Phone/Fax: +39 051525926/+39 0516497501
E-mail:    dcesari <at> arpa.emr.it
Address:   ARPA-SIM, Viale Silvani 6, 40122 Bologna, Italy
========================================================================

This message has been scanned for malware by Websense. www.websense.com
_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj```
8 Jun 2010 00:07

### Fw: postgis-users Digest, Vol 98, Issue 6

```Hi All,

How can I determine the proper PROJ4 parameters considering I have the
following projection parameters for my data files (see PARAM for the North
American CRCM grid herein below my signature)? Is it documented somewhere
within POSTGIS documentation? Or somewhere within PROJ4 documentation?

Here's the best hack I can come up with for the PROJ4 definition:
Proj4js.defs["SR-ORG:6728"]="+proj=stere +lat_0=49.65698 +lat_ts=60
+lon_0=-96.99443 +k=1 +x_0=-1300000 +y_0=100000 +ellps=WGS84 +datum=WGS84
+units=m +no_defs"; But it is not valid. Need some help in adjusting these
parameters. Are the basics PROJ4 parameters documented somewhere?.

Stephane

--------------------------------------------------------------------
Stephane Poirier, M.Sc. Optical Physics
Remote Sensing Application Software Developer
Geography Department, Burnside Hall
http://geog.mcgill.ca
Stephane.Poirier2 <at> mcgill.ca
Cell. (514) 994-3532
Fax: (514) 509-8833
--------------------------------------------------------------------

PARAM for the North American Canadian Regional Climate Model (CRCM) grid:
"
float yc(yc) ;
yc:units = "m" ;
yc:long_name = "y-coordinate of polar-stereographic
projection" ;
yc:standard_name = "projection_y_coordinate" ;
yc:axis = "y" ;
yc:coordinate_defines = "point" ;
yc:actual_range = 22500.f, 7807500.f ;
float xc(xc) ;
xc:units = "m" ;
xc:long_name = "x-coordinate of polar-stereographic
projection" ;
xc:standard_name = "projection_x_coordinate" ;
xc:axis = "x" ;
xc:coordinate_defines = "point" ;
xc:actual_range = 22500.f, 8167500.f ;
float lon(yc, xc) ;
lon:units = "degrees_east" ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
lon:actual_range = 160.4761f, 332.2443f ;
float lat(yc, xc) ;
lat:units = "degrees_north" ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
lat:actual_range = 11.56135f, 87.47531f ;
char polar_stereographic ;
polar_stereographic:straight_vertical_longitude_from_pole =
245.f ;
polar_stereographic:standard_parallel = 60.f ;
polar_stereographic:false_easting = -2790000.f ;
polar_stereographic:false_northing = -8113500.f ;
polar_stereographic:hemisphere_of_standard_parallel = 1.f ;
polar_stereographic:resolution_at_standard_parallel = 45000.f
;
"

Stephane

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

```
8 Jun 2010 08:15

### Re: Fw: postgis-users Digest, Vol 98, Issue 6

```On 06/08/10 00:07, Stephane Poirier wrote:
> Hi All,
>
> How can I determine the proper PROJ4 parameters considering I have the
> following projection parameters for my data files (see PARAM for the North
> American CRCM grid herein below my signature)? Is it documented somewhere
> within POSTGIS documentation? Or somewhere within PROJ4 documentation?
>
> Stephane

Bonjour Stéphane,

You can find the PROJ4 parameters for this projection on the
spatialreference.org web site :
http://spatialreference.org/ref/sr-org/6728/proj4/
I am not sure they are valid, because they miss the +towgs84 parameter.

The documentation of the PROJ4 parameters is at
http://trac.osgeo.org/proj/wiki/GenParms

Jean-Claude
_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

```
8 Jun 2010 09:40

### Re: TM06 for Portugal

```Roger Oberholtzer wrote:

> How would you convert WGS84 lat/long into TM06
> northing/easting? I know that WGS84 is so close
> to ETR89 here that one 'could' use the WGS84 values
> in the transform from ETRS89 to TM06. But how
> would I do it the 'proper' way with proj?

For this purpose, there are seven-parameter
transformations that are time-dependent (since
the ETRS89 follows the contintental drift of
Eurasia, while WGS84 does not).

The British Ordnance Survey cites the following
datum shift from ITRS2005 to ETRS89:

tX =  0.056 m
tY =  0.048 m
tZ = -0.037 m
rX =  0.000054 sec * deltatime
rY =  0.000518 sec * deltatime
rZ = -0.000781 sec * deltatime
s  =  0 ppm

where deltatime is the difference, in years, between
the time when the ITRS2005 coordinates were observed,
and the fixed epoch of ETRS89, which is 1989.0.
http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesyst
emsinfo/guidecontents/guide6.html

For example, today its June 8, 2010, which is the
159th day of the year, so expressed in decimal years
it is today 2010 + 159/365 = 2010.44.  So deltatime
would be 2010.44 - 1989.0 = 21.44.

I believe ITRS2005 can be regarded as equivalent to
WGS84 with centimeter accuracy (although I am a bit
out of my depth here).  If so, the datum shift can
be interpreted as going from from WGS84 to ETRS89,
and in Proj.4 we want to express it in the opposite
direction.  So we reverse all signs.  (The Ordnance
Survey uses the same rotation sign convention as
Proj.4, so that is one worry less.)

So, as an accurate datum shift for ETRS89, we get

+towgs84=-0.056,-0.048,0.037, < -0.000054 * 21.44 >, < -0.000518 *
21.44 >, < 0.000781 * 21.44 >, 0

or

+towgs84=-0.056,-0.048,0.037,-0.00115776,-0.01110592,0.01674464,0

At least today.

transformation.  I believe they are known as
"14 parameter transformations" in general.
Also, there may be regional versions that are
more accurate for Portugal.

Disclaimer: I am not a geodesist, and I haven't
tested this datum shift on any real data.

Mikael Rittri
Carmenta AB
Sweden
www.carmenta.com

-----Original Message-----
From: proj-bounces <at> lists.maptools.org
[mailto:proj-bounces <at> lists.maptools.org] On Behalf Of Roger Oberholtzer
Sent: Monday, June 07, 2010 1:38 PM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] TM06 for Portugal

(full text at
http://lists.maptools.org/pipermail/proj/2010-June/005233.html )

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

```

Gmane