[GRASS5] Re: What proj library version is best?

GRASS GIS GRASS GIS developer
1

Hello Thierry,
(cc to Kergis, GRASS and PROJ lists and Maciek who was also interested in this)

First of all a summary of how I understand projections and co-ordinate systems and the C language structures used in the various PROJ derivatives to represent these.

'Pure' PROJ performs forward and inverse projections from latitude/longitude angles to co-ordinates in the projected system described by the 'PJ' structure. A PJ struct (created using the pj_init() function) includes a name (e.g. proj=tmerc), ellipsoid parameters (a, es or b etc.) and other projection-specific parameters such as latitude of origin, false easting and northing etc. Something like (from the original projects.h in GRASS CVS):
         /* base projection data structure */
typedef struct PJconsts {
         XY (*fwd)(LP, struct PJconsts *);
         LP (*inv)(XY, struct PJconsts *);
         void (*spc)(LP, struct PJconsts *, struct FACTORS *);
         void (*pfree)(struct PJconsts *);
         const char *descr;
         paralist *params; /* parameter list */
         int over; /* over-range flag */
         int geoc; /* geocentric latitude flag */
         double
                 a, /* major axis or radius if es==0 */
                 e, /* eccentricity */
                 es, /* e ^ 2 */
                 ra, /* 1/A */
                 one_es, /* 1 - e^2 */
                 rone_es, /* 1/one_es */
                 lam0, phi0, /* central longitude, latitude */
                 x0, y0, /* easting and northing */
                 k0, /* general scaling factor */
                 to_meter, fr_meter; /* cartesian scaling */
#ifdef PROJ_PARMS__
PROJ_PARMS__
#endif /* end of optional extensions */
} PJ;

Forward and inverse projection are two specific mathematical operations that form part of a co-ordinate system conversion but do not do this in themselves.

On Mon, 26 Apr 2004, Thierry Laronde wrote:

Hello Paul,

The original CERL sources included a version of the PROJ.4 library being
provided supplementary to another derived version used for SCS mapgen
(mapgen has not been updated to use the more recent (as of 1995) PROJ.4
and Dave Gerdes decided to provide the last PROJ.4 as an add-on).

This is interesting to me as I had thought the extensions to the PROJ version in GRASS (two extra files, get_proj.c and do_proj.c and some additions to projects.h) were GRASS-specific but it makes more sense that they were part of Mapgen and got integrated into GRASS along with other Soil Conservation Service enhancements. It is also consistent with an e-mail I remember reading from Gerald Evenden in the PROJ mailing list archive asking was Mapgen still part of GRASS... that suggested to me he was interested in Mapgen because it included some kind of extension to (some old USGS version of) PROJ.

This extension was to add an API (functions pj_get_kv() and pj_do_proj() primarily) for performing proper co-ordinate system conversions rather than just forward and inverse projections. A general co-ordinate system conversion might involve peforming an inverse projection into latitude/longitude and then a forward projection into a different system, at the same time accounting for perhaps different units used in the two projected systems. A co-ordinate system is described by the pj_info struct, which was added to projects.h:
struct pj_info {
       PJ *pj;
       double meters;
       int zone;
       char proj[100];
};

The string pj_info->proj is used to hold the PROJ name for the projection (e.g. tmerc) if it is a projected co-ordinate system, or proj=ll for latitude/longitude. This addition is needed because of course latitude/longitude can be a co-ordinate system in itself as well as an intermediate stage between two projected systems.

pj_info->meters holds a conversion factor for the projected units to meters (set by unfact=x in the parameter pairs passed to the initialising function pj_get_(string|kv)).

I'm really not sure what the zone is for and haven't seen it used anywhere. Finally pj_info->pj is a pointer to a PJ struct (as described above) for the projected co-ordinate system. If it is lat/long then this is uninitialised; zone and meters are set to 0 and 1.0 respectively. The pj_info struct is usually created by pj_get_kv(), pj_get_string() or pj_zero_proj(), but there were many examples in the GRASS source code where it was created manually for a simple lat/long co-ordinate system.

The author of PROJ.4 has, in 2003, extract the core libraries and made
some improvements and fixes in what he calls libproj4.

I am fairly sure if you take an old version of GRASS (before 5.0.0pre4) you could simply drop the files from Gerald's libproj4 in alongside the get_proj.c and do_proj.c and everything will work fine. I don't think the Mapgen version modified anything in PROJ, only added those two files and function prototypes and struct definition to projects.h.

The problem with that of course is that it won't do datum transformation, which as far as I can make out from discussions on the GRASS mailing lists is a requirement for most of the uses people want to put the re-projection modules to (i.e. [rsv].proj). There are also other modules in GRASS that use PROJ but don't need datum transformation, just conversions to/from latitude and longitude angles, e.g. r.sun etc.

Remotesensing.org has taken (old) PROJ.4 maintainance.

Not just maintenance but the extremely useful addition of general co-ordinate system conversions incorporating datum transformation. The interesting thing is comparing the approach Frank Warmerdam took in extending PROJ to handle general co-ordinate system conversions to that chosen by the authors in the SCS. I feel Mapgen kept more to the spirit of the PJ struct by retaining it as description of only a projected system, and adding their own 'wrapper' struct for incorporating details of the units etc. In remotesensing.org proj the PJ struct itself was extended from describing purely a projection, to
describing a complete co-ordinate system. proj=longlat or latlong is used in
place of the Mapgen ll, and to_meter in place of unfact for describing the
unit conversion factor. In fact I think this to_meter may have been in place
in later versions of USGS PROJ than the one GRASS/Mapgen was based on.

The datum transformation parameters are considered an attribute of the
co-ordinate system and if the datum transformation parameters are different
then two co-ordinate systems are considered different, even if the
projection parameters and ellipsoid settings are the same. This can
sometimes result in projecting a point between two identical co-ordinate
systems resulting in an (erroneous) change because they have different
datum parameters. This is now confirmed as a feature of the remotesensing.org
PROJ.4 implementation (see PROJ bug no. 368
http://bugzilla.remotesensing.org/show_bug.cgi?id=368 )

I tend to want to put libproj4 in place of old PROJ.4, but since this is
a field you know far better than me I'd like to have your advice.
libproj4 is the library minus the programs. But I'm not fear anymore of
writing new things or updating old ones...

Yes well as I said above, I think libproj4 should slot in easily enough to
any version of GRASS prior to 5.0.0pre4, so that certainly applies to the 4.x
version KerGIS is based on. Then in the future perhaps datum transformation
could be implemented with the datum information as attributes of the pj_info
(i.e. co-ordinate system) struct rather than the PJ (i.e. projection) struct.

You (or me or whoever ends up making the change) could possibly even re-use
Frank's datum transformation code from remotesensing.org PROJ, in a
re-arranged format.

TIA.

Yes, that was well in advance :slight_smile: Sorry for taking so long to reply but I wanted
to try and make my answer clear and useful to other people and a lot of work
intervened.

PS: I have read your page about datum conversions for France. If you
need to read documentations in french and find parts you are
uncomfortable with, tell me I will translate them for you/answer your
questions (about the meaning of the french writing of course! Not about
what the text is talking about :wink:

I still have not had any feedback on how useful this file is so I will give
the web address here in case anyone from France is interested:
http://www.stjohnspoint.co.uk/gis/france.htm

Good luck in the continuing efforts with KerGIS. I know there is a lot of
background work in tidying up the code structure and preparing it all for
one day when it will be used by many other developers in diverse GIS projects;
hopefully this will not be too far away.

Paul

2

Hello Paul,

On Mon, Jul 05, 2004 at 10:55:15AM +0100, Paul Kelly wrote:

Hello Thierry,
(cc to Kergis, GRASS and PROJ lists and Maciek who was also interested in
this)

[snipped]

While the documentation is still under construction, if you agree I will
put a copy of your message under the "documentation" part of the KerGIS
site, so that people can have access to your explanations (and,
hopefully, so that french people at least can read your page about
projection studies about french systems [1]).
.
This is meant to be a temporary solution since I do believe that I (and
others) need to provide an up-to-date and summarizing documentation
pointing to more in-depth discussion of important points [and once 1.0
is obtained I will focus on the documentation].

And indeed projection management is an essential part of a GIS software
so that is why the proj library has to be provided by KerGIS (and not
put as an external dependency). So from your explanations and my
feelings :

>
> I tend to want to put libproj4 in place of old PROJ.4, but since this is
> a field you know far better than me I'd like to have your advice.
> libproj4 is the library minus the programs. But I'm not fear anymore of
> writing new things or updating old ones...

Yes well as I said above, I think libproj4 should slot in easily enough to
any version of GRASS prior to 5.0.0pre4, so that certainly applies to the 4.x
version KerGIS is based on. Then in the future perhaps datum transformation
could be implemented with the datum information as attributes of the pj_info
(i.e. co-ordinate system) struct rather than the PJ (i.e. projection) struct.

You (or me or whoever ends up making the change) could possibly even re-use
Frank's datum transformation code from remotesensing.org PROJ, in a
re-arranged format.

I do prefer this solution of handling datum through attributes via
pj_info since it will allow:

1) a more fine grained handling of conversions without the side-effect
of the PJ integration in the present PROJ.4

2) an extension that will be compatible enough so it can be scheduled
for the 1.x serie of KerGIS.

Since I have already discovered that some of my initial feelings about
future implementations (for example of the textual attributes
management) have to be revised with the understanding I have gained of
the logic of the historical GRASS code, the goal is to experiment,
use intensively and extend to the limits _inside its present logic_ the
4.0 kernel in order to obtain clear indications about what are the
bottlenecks and the impossibilities and how some changes will impact the
whole.

Handling of the datum will be a neat progress, but will not impose a
complete revolution :wink: so I'm for it. And there is a lot of good stuff
that can fit easily in the actual base of the code.

Good luck in the continuing efforts with KerGIS. I know there is a lot of
background work in tidying up the code structure and preparing it all for
one day when it will be used by many other developers in diverse GIS projects;
hopefully this will not be too far away.

Thanks! We are at a couple of months from the whole sources fixed
(alpha-5) but I want to be more conservative with the 1.0 since the
directory tree, the function naming scheme and the organization of the
libraries shall be stable for the 1.x serie and I will try to minimize
the mistakes (the dramatic changes, not in the logic but in the
organization must be done now so that people can rely on stable features
while improvements continue to come in).

Thanks for taking time to answer thoroughly.

Cheers,

1: The projection stuff should be study carefully by all people dealing
with georeferencing, and the comparison and distorsions studies for
example in France with different projection systems ---as is illustrated
by your paper--- should convince people to use libre GIS software to
visualize them, simply because these have practical consequences.
I once have seen people puzzled because the
establishment on the ground of a building more than 400 meters long was
giving, when verifying the coordinates of the corners, "wrong" results.
The explanation was that the coordinates were in Lambert II, so with a
conformal conic projection i.e. with a projection that keeps the angles
but not the distances. And the establishment on the ground by the
construction industry was using true distance (local coordinates system)
while the surveyor was putting everything in Lambert II. With 400 meters
or more there were several cm of difference...
And going form Lambert II étendu to Lambert-93 should make clear that
great care has to be taken...
--
Thierry Laronde (Alceste) <tlaronde@polynum.com>
http://www.kergis.org/ | http://www.kergis.com/
Key fingerprint = 0FF7 E906 FBAF FE95 FD89 250D 52B1 AE95 6006 F40C

3

This discussion reminds me of the email interaction Frank Warmerdam and I
had about what should be included in a projection library. My position was,
and still is, that *only* the projection (conversion between geographic and
cartesian space) is the domain of the library. All non-projection elements like
datum shifts and x-y-z transformations are separate and unrelated problems
and belong in their own library systems.

The reason for this is twofold: 1) projections have uses that do not involve
datums and their is no reason to burden these uses with elements that are
not involved and 2) developers often only have expertise in only one of the
areas and should not have to burden themselves with library maintenance
involving unfamiliar elements.

I may be partly to blame for the tendency to combine projection and datum
manipulation software by including the 'nad2nad' program with the initial
PROJ.4 distributions. But this was included as an example of how to use
the projection library and used by our local shop for datum operations.

I believe Frank and I did reach a consensus but he has not had time to
redo Remotesensing software to conform to the more recent libproj4
that I support. The last I checked, libproj4 has several more projections
(~135 total) than the Remotesensing version including some newly supported
grid systems.

I believe a great deal of the apparent current chaos could have been
avoided by maintaining pristine PROJ.4 libraries and not indulging in
parochial additions that diminish the ease of the importing of later versions.

_____________________________________
Jerry and the low riders: Daisy Mae and Joshua

4

Hello Gerald,

On Mon, Jul 05, 2004 at 08:00:11PM -0400, Gerald Evenden wrote:

This discussion reminds me of the email interaction Frank Warmerdam and
I
had about what should be included in a projection library. My position
was,
and still is, that *only* the projection (conversion between geographic
and
cartesian space) is the domain of the library. All non-projection
elements like
datum shifts and x-y-z transformations are separate and unrelated
problems
and belong in their own library systems.

The reason for this is twofold: 1) projections have uses that do not
involve
datums and their is no reason to burden these uses with elements that
are
not involved and 2) developers often only have expertise in only one of
the
areas and should not have to burden themselves with library maintenance
involving unfamiliar elements.

I do agree with your points. More generally ---and that is the goal of
the KerGIS project [restarted from last CERL sources of GRASS]--- I'm
fanatic about size of the code (no duplication), orthogonality and Unix
spirit: small programs doing one thing but well.

So your libproj4 will be shipped with KerGIS (due credit, licence notes
being specified on the web site and in the sources as for every portion
of code) and no modifications will be done to libproj4 to "customize" it.
The datum handling, as suggested by Paul and you, will be handled
supplementary to this.
libproj4 will be shipped with KerGIS just in order for people retrieving
it to have a standalone GIS (no need to access external ressources not
included in the tarball). But since it will be used as is, it should not
be difficult for us to keep in sync.

I will send you a message once the library is "officially" included in
the KerGIS sources (should be done in october 2004, perhaps sooner since
this does not represent at the moment huge modifications of the
surrounding KerGIS code).

Best regards,
--
Thierry Laronde (Alceste) <tlaronde@polynum.com>
http://www.kergis.org/ | http://www.kergis.com/
Key fingerprint = 0FF7 E906 FBAF FE95 FD89 250D 52B1 AE95 6006 F40C

5

Gerald Evenden wrote:

This discussion reminds me of the email interaction Frank Warmerdam and I
had about what should be included in a projection library. My position was,
and still is, that *only* the projection (conversion between geographic and
cartesian space) is the domain of the library. All non-projection elements like
datum shifts and x-y-z transformations are separate and unrelated problems
and belong in their own library systems.

The reason for this is twofold: 1) projections have uses that do not involve
datums and their is no reason to burden these uses with elements that are
not involved and 2) developers often only have expertise in only one of the
areas and should not have to burden themselves with library maintenance
involving unfamiliar elements.

[...]

I don't know how to interpret this : Datum is an important part of coordinate's
transformation. The same projection on two diferent datums can lead to errors
of 500 m (at least in France, between the NTF and WGS84 Datum). Projection
and Datum transformation are two different steps (and hence may be done in
different functions), but you cannot forget Datum transformation. For exemple,
the correct way to tranform UTM/WGS84 to NTF/Lambert coordinate is :

(E,N) UTM -> lat/long WGS84 -> X,Y,Z (tridimentional) -> Datum shift, rotation
and scale from WGS84 to NTF -> lat/long NTF -> (E,N) Lambert.

Even with this formula you have an accuracy of 5 m in some places : You have to
use a correction matrix to improve the precision to a geodetic one (5 cm).

You can only forget datum when you work at very small scale (1:2,000,000 and under)

--
Michel Wurtz - Auzeville-Tolosane

6

Sorry, but I have been through this issue several times but I see it is still
necessary to reiterate my position: cartographic projections are a subject
area unto themselves and have *no* relationship to datums nor other
aspects of common cadastral and geodetic activities. One can spend
their entire life studying and developing cartographic projections and
not have one whit of knowledge of datums nor how to manipulate them.

I take John Snyder as an excellent example: he wrote many articles and
books related to cartographic projections. He created new projections.
He was recognized as a world class expert on cartographic projections.
However, you do not need to take your shoes off to count the number
of paragraphs he wrote on datum operations. In fact, I cannot think
of one case other than his listing of the projection parameters for
both NAD27 and NAD83. Surely, Snyder may have been aware of
datum shifting and how to do it but he did not seem to include it in
his research activities.

Surely, projections are required in the practical world of cadastral and
geodetic activities, but so are sine and cosine and a bunch of other
transcendental functions. The software development business
normally compartmentalizes solutions by individuals specializing
in solving different parts of a larger problem with the result of a
collection of libraries that are linked together to produce a more
general application program that utilizes a wide area of disciplines.

Lastly, the argument that because the major axis and eccentricity are
major details of a datum and that they are required by a projection, then the
projection should be included as part of a datum library just obfuscates
the issue. This whole subject area uses the sine function in many
places---does that mean that we have to include the sine function in a
cadastral library? I think not.

I am sorry to bore y'all but I think this is an important issue and needs
to be understood.
_____________________________________
Jerry and the low riders: Daisy Mae and Joshua

7

I just wanted to say that you cannot propose a GIS without offering a
*correct* function to transform a map from one coordinate system to
another one, even if the Datum involved is not the same. The problem
is not "should we develop this as internal or external library" : that
doesn't matter. The problem is "I have many sources of data, with an
accuracy of a few meters. I want to display/compare/mix/analyse them
together. My GIS should manage that".

So the coordinate transformation functions must take care of datums.

If Grass can't do that, Gras is unusable for a vast majority of application,
and we will not have any usable open-source GIS that can stand in front of
ArcView, MapInfo, GeoConcept and other commercial GIS softwares...

Gerald Evenden wrote:

Sorry, but I have been through this issue several times but I see it is still
necessary to reiterate my position: cartographic projections are a subject
area unto themselves and have *no* relationship to datums nor other
aspects of common cadastral and geodetic activities. One can spend
their entire life studying and developing cartographic projections and
not have one whit of knowledge of datums nor how to manipulate them.

I take John Snyder as an excellent example: he wrote many articles and
books related to cartographic projections. He created new projections.
He was recognized as a world class expert on cartographic projections.
However, you do not need to take your shoes off to count the number
of paragraphs he wrote on datum operations. In fact, I cannot think
of one case other than his listing of the projection parameters for
both NAD27 and NAD83. Surely, Snyder may have been aware of
datum shifting and how to do it but he did not seem to include it in
his research activities.

Surely, projections are required in the practical world of cadastral and
geodetic activities, but so are sine and cosine and a bunch of other
transcendental functions. The software development business
normally compartmentalizes solutions by individuals specializing
in solving different parts of a larger problem with the result of a
collection of libraries that are linked together to produce a more
general application program that utilizes a wide area of disciplines.

Lastly, the argument that because the major axis and eccentricity are
major details of a datum and that they are required by a projection, then the projection should be included as part of a datum library just obfuscates
the issue. This whole subject area uses the sine function in many
places---does that mean that we have to include the sine function in a
cadastral library? I think not.

I am sorry to bore y'all but I think this is an important issue and needs
to be understood.
_____________________________________
Jerry and the low riders: Daisy Mae and Joshua

--
Michel Wurtz - Auzeville-Tolosane

8

Hello Michel

On Thu, 8 Jul 2004, Michel Wurtz wrote:

I just wanted to say that you cannot propose a GIS without offering a
*correct* function to transform a map from one coordinate system to
another one, even if the Datum involved is not the same. The problem
is not "should we develop this as internal or external library" : that
doesn't matter. The problem is "I have many sources of data, with an
accuracy of a few meters. I want to display/compare/mix/analyse them
together. My GIS should manage that".

Gerald is talking about a projection library (PROJ), not a GIS (GRASS). Of course GRASS can do datum transformations but it is perfectly feasible to only use the PROJ library for projection calculations and to implement datum shifting separately. This is what Thierry is doing with KerGIS (another
GIS, based on GRASS).

I would also add that within GRASS there are examples of the utility of projection software that do not involve datums, e.g. in calculating the latitude and longitude angles of a position on the earth for use in sun angle calculations (r.sun etc.).

So the coordinate transformation functions must take care of datums.

If Grass can't do that, Gras is unusable for a vast majority of application,
and we will not have any usable open-source GIS that can stand in front of
ArcView, MapInfo, GeoConcept and other commercial GIS softwares...

Just to clarify for other people reading, GRASS *can* do that. The point of this discussion was more a low-level technical one about how this support is organised and programmed within the GIS and not a discussion of the high-level functionality visible to the user.

Paul

9

Hello,

On Thu, Jul 08, 2004 at 12:03:28AM +0200, Michel Wurtz wrote:
[..]

the correct way to tranform UTM/WGS84 to NTF/Lambert coordinate is :

(E,N) UTM -> lat/long WGS84 -> X,Y,Z (tridimentional) -> Datum shift, rotation
and scale from WGS84 to NTF -> lat/long NTF -> (E,N) Lambert.

As already explained by Gerald Evenden and Paul Kelly, your
"transformation formula" is a typical pipe i.e. a combination of
different tools. Since the beginning and end
are reverse functions, they belong to the very same library: Gerald
Evenden's PROJ.4 or libproj4 (the new version).

The datum handling deserves its proper tools (library). The discussion
was not about user visible result, but how the code shall be organized.

GPL GRASS is already providing such tools (minus some problems in some
cases linked precisely to the fact that proj and datum are merged in a
library).
Since KerGIS is new I have the possibility (and the duty) to choose what
I think is the best way to implement features. This best is for me to
separate clearly pieces that are logically distinct, allowing to rectify
the curve of needs (approaching almost exactly a "spline" via the
smallest segments (programs)) and to insulate problems.
--
Thierry Laronde (Alceste) <tlaronde@polynum.com>
http://www.kergis.org/ | http://www.kergis.com/
Key fingerprint = 0FF7 E906 FBAF FE95 FD89 250D 52B1 AE95 6006 F40C