[GRASS5] Fixing up projection related code...

GRASS GRASS developer
1

I've been looking at some of the projection related code, and I'm
trying to come up with a course of action to clean it up in places.
Eventually, this should result in:

1. Having better defined projection files.
2. Creation of a set of predefined well known coordinate systems
    files.
3. A rewrite of g.setproj with these ideas in mind:
    - Types of coordinate system are XY, geographic, or projected
    - Ask users if they want to use a predefined coordinate
       system (state plane, UTM, others...) or if they want to
       define their own. Having a collection of predefined
       coordinate systems where the user can just pick it by
       "name", will simplify user's lives as well as insuring the
       info is correctly specified.
    - Both projected and geographic require datum/ellipsoid info.
       Ask for datum first, ellipsoid parameters can be looked up
       from it. Otherwise, just get ellipsoid parameters.
    - Remove state plane from projection list, it'll be covered
       by the predefined systems above (though still need units...)
    - Allow user defined datums??
    - Other...?
4. Have transparent handling for a variety of datum transformations
    incorporated into the library wrappers for proj4. Most of the
    necessary work for the transformations has been done by Andreas
    Lange, it's a matter of working in the strategies. Perhaps the
    programs can pass on a preferred stategy when specified by the
    user (for instance, NADCON vs. Bursa Wolf). Get rid of GRASS
    copy of proj4 and make a dependency on proj4 >= whatever version.
    It should be a fatal error if a datum transformation is impossible.
5. Miscellaneous:
    - Specify vertical datums? MSL, ellipsoid height, others... Mostly
       for documentation, since very few can be translated.
    - Add missing usfoot/usfeet <-> metric conversions. usinch? usmile?
    - Add prime meridian as datum parameter.
    - Other?

So, I hacked up a little perl script to convert the state27 and state83
files to PROJ_INFO files with all of the datum/ellipsoid parameters
present (current g.setproj fails on this). In thinking about this, I
wondered why PROJ_UNITS exists at all. Silly for a three line file.
That info should be in the PROJ_INFO file, IMHO. Anyway, I thought
we could have a Coordinate_System directory that g.setproj could
use. So far, with my state plane files, I have something like:

Coordinate_Systems/
   Projected/
      North_American/
         State_Plane_1927/
      State_Plane_Alabama_East
      State_Plane_Alabama_West
      ...
      State_Plane_Wyoming_West_Central
   State_Plane_1983/
      ...

This kind of thing is easy enough to add to, and g.setproj could
just use directory traversals until the user selects a file, then
just copy the file to the new location...

Example generated State Plane file:

# 5010:State_Plane_1927_Alaska_Zone_10
name: State_Plane_1927_Alaska_Zone_10
ellps: clark66
y_0: 0
es: 0.00676865799729109
a: 6378206.4
dx: -22
dy: 157
dz: 176
lat_0: 51
f: 294.9786982
lat_1: 53.8333333333333
lat_2: 51.8333333333333
proj: lcc
lon_0: -176
x_0: 914401.828803658

I would add the following for STP 1927:
unit: usfoot
units: usfeet
meters: 0.30480060960121920243

For STP 1983, it's case by case for the above or just meters...

Anyway, some of the fixing up will probably break programs all over
the place. I don't want to see switch statements like:

  switch (G_projection()) {
     case XY:
        ...
     case UTM:
        ...
     case STP:
        ...
     case LL:
        ...
     default:
        do_other_thing();
  }

I'd rather see something like:

  switch (G_projection_type()) {
     case PROJ_XY:
        ...
     case PROJ_GEOGRAPHIC:
        ...
     case PROJ_PROJECTED:
        ...
     default:
        some_error_func();
  }

The problems with the first are STP is poorly defined, what's so
special about UTM? and you still don't have very good info about
the projection if you really need it -- in which case, you have
to get the projection key/vals anyway... But, if you only need
to know if you're working with linear units vs. angular units,
then the second method gives you that (or does it?).

Likewise, the WIND files' idea of projection should change to
reflect this. Currently, the proj: <num> is next to useless
and zone: <num> isn't applicable in most cases. What are the
WIND files for? Just the region settings, north/south/east/west
and resolutions. But arguably, it should be possible to know
the XY vs. LL vs. PROJECTED from this file.

Anyway, I wanted to throw this out for discussion. I'm tired
of reading/responding to people confused about projections in
GRASS (and they are more confusing than need be).

--
Eric G. Miller <egm2@jps.net>

2

On Mon, 26 Aug 2002, Eric G. Miller wrote:

I've been looking at some of the projection related code, and I'm
trying to come up with a course of action to clean it up in places.
Eventually, this should result in:

...

4. Have transparent handling for a variety of datum transformations
    incorporated into the library wrappers for proj4. Most of the
    necessary work for the transformations has been done by Andreas
    Lange, it's a matter of working in the strategies. Perhaps the
    programs can pass on a preferred stategy when specified by the
    user (for instance, NADCON vs. Bursa Wolf). Get rid of GRASS
    copy of proj4 and make a dependency on proj4 >= whatever version.
    It should be a fatal error if a datum transformation is impossible.

I have already done half of this. Now to separate the GRASS dependency on
a local copy of the PROJ.4 source, I need to extend the GRASS datum.table
format to allow more than the simple 3-parameter dx dy dz datum
transformation parameters it stores at present.

I intend to use the same keywords as PROJ.4, i.e. towgs84 for the
7-parameter transform (and optionally for the 3-parameter[1] ) and
nadgrids for the (NADCON?) American grid-shifting tables. I intend to
allow more than one set of transform parameters for each datum (it doesn't
seem right to delete anything from the table), but what I'm not sure about
is which should be used.

Should the default logic be to use the most accurate transformation (tables)
if it's available, then 7-parameter, falling back to 3-parameter if there
is nothing else in the datum.table?

Eric suggested above that programs could pass on a preferred type of
transform, but what if that's not available in the table? It would seem a
bit awkward to implement, making the system less transparent. Perhaps the
new improved g.setproj could allow the user to set which transformation
they wanted before running the re-projection command, and this could be
explained in the documentation for those commands?

That's what I'm leaning towards at the minute, but I haven't convinced
myself it's the best way. I was wondering does anyone know how it is
implemented in other GIS systems?

To summarise what I'll do if there are no better suggestions:
1. By default use the most accurate transformation available in the GRASS
datum.table
2. Explain this in the man pages / documentation for the re-projection
commands, and explain how the user can define his/her own datum
transformation parameters.

[1] In PROJ.4, 'towgs84=' should be followed by a comma-separated list of
numbers. If there are 3 numbers a 3-parameter transform is assumed; if 7
it is assumed a 7-parameter transform was intended.

Paul.

3

Paul Kelly wrote:

Should the default logic be to use the most accurate transformation (tables)
if it's available, then 7-parameter, falling back to 3-parameter if there
is nothing else in the datum.table?

Paul,

First an observation. It is hard in many cases to establish which is
the more accurate transformation. For instance if you have two three
parameter transformations, which is best? The answer is that likely
either is best for a particular area. That is, it is common to generate
an assortment of datum shifts, each appropriate for some subarea where
the datum is used. Even in cases like datum shift files where some
are high precision, but local to a region, it is hard to automatically
determine which is best.

This is not a problem I have dealt with well in PROJ, just punting and
leaving it up to the user to select an appropriate transformation.

In EPSG they list all the transformations, and the descriptions generally
contain information useful to the user in selecting the best for their
situation.

The most general solution then would likely be maintaining a list of
datum shift options for any given datum, each having some descriptive
text that can be offered to the user when setting up their location
(in g.setproj?).

A slightly easier approach would be to have a datum name for each
possible transformation. This is what they do in MapInfo for instance.
So you might have nad27-conus, nad27-ntv1, nad27-HPGN, each represening
a different approximation for shifting to WGS84. This is likely what
I would suggest in your case. The downside is that if a user select nad27
it may be difficult to automatically determine a list of all the datum
names that could apply.

Eric suggested above that programs could pass on a preferred type of
transform, but what if that's not available in the table?

I do think that import programs that report projection information, or
offer to setup new locations (like r.in.gdal) should offer the transformation
used if possible.

Best regards,

--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

4

Paul,

I have latest CVS GRASS and have found a couple of projection programs
are not quite working. I discovered the problem when using m.ll2db which
converts to and from GRASS database and Lat. Long. coordinates. When I
run the program I get the following error ...
cannot initialize pj
cause: major axis or radius = 0 or not given

The problem is with how the "ll" projection is initialized, namely ...
parms_in[0] = '\0';
pj_get_string(&info_in, parms_in);

If I change this to a better qualified initialization such as ...

sprintf(ll_buf, "+proj=ll +ellsp=wgs84"); /* +datum also works with
current setting */
pj_get_string(&info_in, ll_buf);

it works fine. This problem has implications with other programs such as
m.proj(2), and v.in.gshhs. There may be others with this type of
initialization.

You seem to be a few steps ahead on the PROJ stuff. Should the code be
fixed (m.proj seems to be on the chopping block) or do you have other
suggestions?

--
Bob Covill

Tekmap Consulting
P.O. Box 2016
Fall River, N.S.
B2T 1K6
Canada

E-Mail: bcovill@tekmap.ns.ca
Phone: 902-860-1496
Fax: 902-860-1498

5

Paul,

I just noticed on previous message that +ellps was incorrectly stated as
+ellsp.

Sorry about that.

--
Bob Covill

Tekmap Consulting
P.O. Box 2016
Fall River, N.S.
B2T 1K6
Canada

E-Mail: bcovill@tekmap.ns.ca
Phone: 902-860-1496
Fax: 902-860-1498

6

Hello Bob

On Mon, 3 Feb 2003, Bob Covill wrote:

Paul,

I have latest CVS GRASS and have found a couple of projection programs
are not quite working. I discovered the problem when using m.ll2db which
converts to and from GRASS database and Lat. Long. coordinates. When I
run the program I get the following error ...

BTW I couldn't find a man page for this program; could you give me an
example of a typical command line you might use with it?

cannot initialize pj
cause: major axis or radius = 0 or not given

Oops. It looks like we need to specify the ellipsoid for the input
projection; it doesn't just assume it is the same as the output if not
specified. I had been testing this behaviour using different combinations
of arguments given to cs2cs, and it seemed happy with a simple +proj=latlong,
but I realise now it must have been silently adding a default ellipsoid to
the input projection definition, to get round this problem.

The problem is with how the "ll" projection is initialized, namely ...
parms_in[0] = '\0';
pj_get_string(&info_in, parms_in);

If I change this to a better qualified initialization such as ...

sprintf(ll_buf, "+proj=ll +ellsp=wgs84"); /* +datum also works with
current setting */
pj_get_string(&info_in, ll_buf);

In this case, and for most of the programs that use the proj library,
there is no datum transformation, so the input ellipsoid should be the
same as the output ellipsoid? Following that logic, could you please test
the following patch for me and see if it gives you your expected results.

Index: src/misc/m.ll2db/main.c

RCS file: /grassrepository/grass/src/misc/m.ll2db/main.c,v
retrieving revision 1.3
diff -u -r1.3 main.c
--- src/misc/m.ll2db/main.c 26 Jan 2003 17:47:01 -0000 1.3
+++ src/misc/m.ll2db/main.c 3 Feb 2003 16:21:35 -0000
@@ -148,11 +148,6 @@
exit(1);
}

-/* In Info */
-parms_in[0] = '\0';
-pj_get_string(&info_in, parms_in);
-
-
/* Out Info */
out_proj_keys = G_get_projinfo();
out_unit_keys = G_get_projunits();
@@ -160,6 +155,12 @@
exit (0);
}

+/* In Info */
+if( G_find_key_value("ellps", out_proj_keys) != NULL )
+ sprintf(parms_in, "proj=ll ellps=%s", G_find_key_value("ellps", out_proj_keys) );
+else
+ sprintf(parms_in, "proj=ll ellps=wgs84");
+pj_get_string(&info_in, parms_in);

     if (isatty(0))
     {

If it works I'll add something similar to the other modules that use proj
with no datum transformation.

Paul

7

Paul,

I tried the patch listed below, and it does the trick. The other two
places I know where the projection is initialized this way are in
m.proj(2), and v.in.gshhs. It may occur in other places as well.

Paul Kelly wrote:

Hello Bob

On Mon, 3 Feb 2003, Bob Covill wrote:

> Paul,
>
> I have latest CVS GRASS and have found a couple of projection programs
> are not quite working. I discovered the problem when using m.ll2db which
> converts to and from GRASS database and Lat. Long. coordinates. When I
> run the program I get the following error ...

BTW I couldn't find a man page for this program; could you give me an
example of a typical command line you might use with it?

The m.ll2db program is derived from m.ll2u (and m.u2ll) and accepts
pretty much the same format. Instead of the user supplying projection
parameters, they are obtained from PROJ_* files. It accepts forward and
reverse projection to and from Lat. Long. coordinates. In other words it
should work in all GRASS Projections with the exception of XY and
LatLong (ll). Also, the program uses the PROJ library instead of the
coorcnv library.

For example a table of Lat. Long coordinates can be converted to the
users current projection with the following ....
m.ll2db in=coords_file.LL out=coords_file.UTM
There are flags for reverse coordinate order, and inverse projection.
Data can also be piped through the program.

> cannot initialize pj
> cause: major axis or radius = 0 or not given

Oops. It looks like we need to specify the ellipsoid for the input
projection; it doesn't just assume it is the same as the output if not
specified. I had been testing this behaviour using different combinations
of arguments given to cs2cs, and it seemed happy with a simple +proj=latlong,
but I realise now it must have been silently adding a default ellipsoid to
the input projection definition, to get round this problem.

> The problem is with how the "ll" projection is initialized, namely ...
> parms_in[0] = '\0';
> pj_get_string(&info_in, parms_in);
>
> If I change this to a better qualified initialization such as ...
>
> sprintf(ll_buf, "+proj=ll +ellsp=wgs84"); /* +datum also works with
> current setting */
> pj_get_string(&info_in, ll_buf);

In this case, and for most of the programs that use the proj library,
there is no datum transformation, so the input ellipsoid should be the
same as the output ellipsoid? Following that logic, could you please test
the following patch for me and see if it gives you your expected results.

Index: src/misc/m.ll2db/main.c

RCS file: /grassrepository/grass/src/misc/m.ll2db/main.c,v
retrieving revision 1.3
diff -u -r1.3 main.c
--- src/misc/m.ll2db/main.c 26 Jan 2003 17:47:01 -0000 1.3
+++ src/misc/m.ll2db/main.c 3 Feb 2003 16:21:35 -0000
@@ -148,11 +148,6 @@
exit(1);
}

-/* In Info */
-parms_in[0] = '\0';
-pj_get_string(&info_in, parms_in);
-
-
/* Out Info */
out_proj_keys = G_get_projinfo();
out_unit_keys = G_get_projunits();
@@ -160,6 +155,12 @@
exit (0);
}

+/* In Info */
+if( G_find_key_value("ellps", out_proj_keys) != NULL )
+ sprintf(parms_in, "proj=ll ellps=%s", G_find_key_value("ellps", out_proj_keys) );
+else
+ sprintf(parms_in, "proj=ll ellps=wgs84");
+pj_get_string(&info_in, parms_in);

     if (isatty(0))
     {

If it works I'll add something similar to the other modules that use proj
with no datum transformation.

Paul

Thanks for your help.

--
Bob Covill

Tekmap Consulting
P.O. Box 2016
Fall River, N.S.
B2T 1K6
Canada

E-Mail: bcovill@tekmap.ns.ca
Phone: 902-860-1496
Fax: 902-860-1498

8

Hi Paul,

I just discovered another program that produces the same PROJ error
reported earlier. When g.region -l is executed suffers the same fate ...
cannot initialize pj
cause: major axis or radius = 0 or not given

Probably we need to search the GRASS code to catch all of the projection
initialization problems.

--
Bob Covill

Tekmap Consulting
P.O. Box 2016
Fall River, N.S.
B2T 1K6
Canada

E-Mail: bcovill@tekmap.ns.ca
Phone: 902-860-1496
Fax: 902-860-1498

9

Bob Covill wrote:

Probably we need to search the GRASS code to catch all of the projection
initialization problems.

The following programs are linked against libproj:

etc/bin/cmd/g.region
etc/bin/cmd/m.proj2
etc/bin/cmd/r.in.gdal
etc/bin/cmd/r.proj
etc/bin/cmd/r.sun
etc/bin/cmd/r.sunmask
etc/bin/cmd/s.proj
etc/bin/cmd/v.build.polylines
etc/bin/cmd/v.in.gshhs
etc/bin/cmd/v.in.tig.lndmk
etc/bin/cmd/v.mkquads
etc/bin/cmd/v.proj
etc/bin/inter/i.points3
etc/bin/inter/m.proj
etc/bin/inter/v.digit

--
Glynn Clements <glynn.clements@virgin.net>

10

Hello again Bob
Thanks for all the testing

On Mon, 3 Feb 2003, Bob Covill wrote:

Paul,

I tried the patch listed below, and it does the trick. The other two
places I know where the projection is initialized this way are in
m.proj(2), and v.in.gshhs. It may occur in other places as well.

Yes, see http://grass.itc.it/pipermail/grass5/2003-January/004562.html for
a longer list. I've committed a fix similar to the patch we discussed to
all those modules.

I suppose it is good practice to always specify an ellipsoid as
well as the projection type anyway; it doesn't really make sense
otherwise. Two GRASS modules already did this, v.mkquads
(src/mapdev/v.mkquads/convert.c) and i.points3
(src/libes/image3/convert_ll.c), so that was very far-sighted of their
authors :-).

I would be grateful if somebody familiar with malloc and free etc. could
have a quick look at the changes I made to v.mkquads and i.points3 (the
two files mentioned above) just to check my implementation.

I also fixed m.proj2 interactive mode. It wouldn't let the user enter an
ellipsoid if the projection was lat/long, which was definitely a bug. It
still doesn't ask for a datum in interactive mode, but I have changed the
man page to mention that.

v.in.tig.lndmk and v.build.polylines don't seem to use any of the proj
wrapper functions, so maybe their dependency on libproj is wrong. They
probably need their projection code checked anyway....

Paul

11

On Mon, Feb 03, 2003 at 10:05:05PM +0000, Paul Kelly wrote:

I would be grateful if somebody familiar with malloc and free etc. could
have a quick look at the changes I made to v.mkquads and i.points3 (the
two files mentioned above) just to check my implementation.

If you're asking about convert.c in v.mkquads, your usage is erroneous.

G_find_key_value() returns a pointer to the value in the key_value
struct array or *NULL*. First, you don't need to allocate any memory;
second, you're clobbering it; Third, if ellps is not found and a NULL
pointer is returned you try to copy "wgs84" like sprintf (NULL, "wgs84);
Oops!

I think what you want is something like:

   const char *ellps_default = "wgs84";
   char *ellps;

   ....
   ellps = G_find_key_value ("ellps", in_proj_keys);
   if (ellps == NULL)
       ellps = ellps_default;
   G_set_key_value ("ellps", ellps, out_proj_keys);
   ....

The G_set_key_value function makes copies of the strings (if needed) for
the key_value array.

I haven't looked at i.points3.

--
echo ">gra.fcw@2ztr< eryyvZ .T pveR" | rot13 | reverse

12

OK, v.mkquads and i.points3 should be more robust now.

On Mon, 3 Feb 2003, Eric G. Miller wrote:

If you're asking about convert.c in v.mkquads, your usage is erroneous.

G_find_key_value() returns a pointer to the value in the key_value
struct array or *NULL*. First, you don't need to allocate any memory;

This was the way it was used in src/image3/convert_ll.c (i.points3), which
usage I copied over to v.mkquads. I suppose its not a good idea to infer a
function's correct usage by the way other modules use it!

second, you're clobbering it; Third, if ellps is not found and a NULL
pointer is returned you try to copy "wgs84" like sprintf (NULL, "wgs84);
Oops!

Yes this was just me; I didn't look over it for long enough. Was just
trying to get the bugfix out too quickly.

Thank you for your help. I've corrected the usage now. It should only be
an extremely unlikely case that a PROJ_INFO file wouldn't contain an
ellps: value anyway. I'll try and tidy up the whole idea of using default
values at a future date sometime.

Paul

13

Hello everyone
I've been looking at improving and tidying the datum transformation support
for a while now and made a few changes in my local copy of GRASS. I would
just like a little more discussion and / or suggestions and hints before
commiting the changes to CVS.

On Sun, 2 Feb 2003, Frank Warmerdam wrote:

The most general solution then would likely be maintaining a list of
datum shift options for any given datum, each having some descriptive
text that can be offered to the user when setting up their location
(in g.setproj?).

This is what I have chosen to go for, implementing a new file
('datumtransform.table') which can contain multiple entries for each
datum. There are 4 fields: datum name, transformation parameters in PROJ.4
syntax, "Where Used" and "Comment". The last two fields should be
descriptive and help the user to pick the appropriate parameters.

It doesn't really make sense the way it is partially currently implemented
in GRASS, to look up a 3-parameter transform from the datum.table and
just use it, because as Frank said each datum may have many different
tranformations associated with it and GRASS shouldn't pretend to be able to
choose the best one for the user when really it needs a human making the
decision. E.g. someone might want to use an old or custom set of parameters
simply for compatibility with existing data.

So after selection, the datum transformation parameters will be stored in
the PROJ_INFO file under the new key 'datumparams'.
E.g. 'datumparams: towgs84=a,b,c' will be equivalent to the old
'dx: a' 'dy: b' and 'dz: c'. The pj_get_kv() wrapper function will still
interpret the old style properly but g.setproj will only write new-style
PROJ_INFO files.

Regarding passing ellipsoid and datum names onto PROJ.4, this will only be
done if GRASS doesn't recognise the names from its own ellipse.table and
datum.table files. Otherwise the underlying parameters will be passed on
but not the names. This will probably solve bug 1047 with the problem with
the PROJ_INFO files r.in.gdal automatically creates but I will have a look
at it soon as well.

A slightly easier approach would be to have a datum name for each
possible transformation. This is what they do in MapInfo for instance.
So you might have nad27-conus, nad27-ntv1, nad27-HPGN, each represening
a different approximation for shifting to WGS84. This is likely what
I would suggest in your case. The downside is that if a user select nad27
it may be difficult to automatically determine a list of all the datum
names that could apply.

Looking back now I realise this method was sort of partially implemented
in GRASS, e.g. there is an 'a-can' datum for Alaska and Canada, separate
from NAD27. These can't really be purged from the table as it would break
compatibility with existing PROJ_INFO files but we should be careful about
adding new datums to the file in the future. There is no problem with
changing entries in the datumtransform.table file though, as this file is
only looked up when creating projection information.

I do think that import programs that report projection information, or
offer to setup new locations (like r.in.gdal) should offer the transformation
used if possible.

I don't really think that programs that import data should try to
reproject it. This is getting too complicated and then there would be too
many modules to try to keep up to date. The conventional GRASS way seems
to be to import the data into a location with the correct projection, then
reproject it to your target location using [rsv].proj.

Currently g.setproj is only interactive. When we create a cmd version it
can have a flag to only change the datum transformation parameters and not
change anything else. There could be a library function to do this as
well, which could maybe be called by the re-projection programs

Some points I would appreciate if people could comment on:
As far as I can make out a datum has a one-to-many relationship with datum
transformations, and one-to-one with ellipsoid and prime meridian. Neither
the datum lists in GRASS nor PROJ.4 include space for an associated prime
meridian so I'm wondering is this correct. I would like to make the GRASS
format as correct as possible. Probably I should investigate how PROJ.4
uses prime meridians but I don't really have that much spare time at the
minute and would appreciate a hint.

If we were to extend the datum.table format we would need to change the
function which reads it.
At present it is like (in src/libes/gis/datum.c)

if (sscanf(buf, "%99s \"%255[^\"]\" \"%255[^\"]\" \"%255[^\"]\"",
            name, params, where_used, comment) != 4)
...
It would seem to be very hard to add extra fields to this without it
causing problems when reading old-style files, but maybe I'm missing
something obvious. Could we just leave out the check for the return value
being 4?

In GRASS 5.1 there can be proper support for co-ordinate systems other
than State Plane:

On Mon, 26 Aug 2002, Eric G. Miller wrote:

That info should be in the PROJ_INFO file, IMHO. Anyway, I thought
we could have a Coordinate_System directory that g.setproj could
use. So far, with my state plane files, I have something like:

Coordinate_Systems/
   Projected/
      North_American/
         State_Plane_1927/
            State_Plane_Alabama_East
            State_Plane_Alabama_West
            ...
            State_Plane_Wyoming_West_Central
         State_Plane_1983/
            ...

This kind of thing is easy enough to add to, and g.setproj could
just use directory traversals until the user selects a file, then
just copy the file to the new location...

This sounds like a good idea but what are directory traversals and how
should they be implemented in g.setproj? There should be a GUI version as
well to make it look user-friendly and nice as g.setproj is a lot of
people's first encounter with a GRASS command (mine anyway and I spent
days over it) and it is really a terrible piece of software and very
off-putting. Is it all right to use Tcl/Tk in GRASS 5.1 or will there be
a different GUI?

If there is no consensus I can go ahead without the prime meridians; I'm
sure it can probably be added later.

Paul

14

Paul Kelly wrote:

This is what I have chosen to go for, implementing a new file
('datumtransform.table') which can contain multiple entries for each
datum. There are 4 fields: datum name, transformation parameters in PROJ.4
syntax, "Where Used" and "Comment". The last two fields should be
descriptive and help the user to pick the appropriate parameters.

It doesn't really make sense the way it is partially currently implemented
in GRASS, to look up a 3-parameter transform from the datum.table and
just use it, because as Frank said each datum may have many different
tranformations associated with it and GRASS shouldn't pretend to be able to
choose the best one for the user when really it needs a human making the
decision. E.g. someone might want to use an old or custom set of parameters
simply for compatibility with existing data.

So after selection, the datum transformation parameters will be stored in
the PROJ_INFO file under the new key 'datumparams'.
E.g. 'datumparams: towgs84=a,b,c' will be equivalent to the old
'dx: a' 'dy: b' and 'dz: c'. The pj_get_kv() wrapper function will still
interpret the old style properly but g.setproj will only write new-style
PROJ_INFO files.

Paul,

Do correctly understand then that each named datum in the datumtransform
table will have a *default* transformation to use with the datum, but that
by explicitly including the transformation in the PROJ_INFO file it is
intended to be easy for users to alter what transformation values are used?

Regarding passing ellipsoid and datum names onto PROJ.4, this will only be
done if GRASS doesn't recognise the names from its own ellipse.table and
datum.table files. Otherwise the underlying parameters will be passed on
but not the names. This will probably solve bug 1047 with the problem with
the PROJ_INFO files r.in.gdal automatically creates but I will have a look
at it soon as well.

Right. I would add that there is a pj_get_def() function in PROJ.4 which can
be used to get an "expanded" definition of a PROJ.4 coordinate system string.
For instance, it will expand a +datum= into the actual transformation values
that would be used with it. This is intended to be useful for applications
that want to find out the magic in certain coordinate system definitions.
This also could be used to fetch back details from lookups from the various
"cookbook" coordinate system files, like the state plane ones, and the EPSG
database.

I do think that import programs that report projection information, or
offer to setup new locations (like r.in.gdal) should offer the transformation
used if possible.

I don't really think that programs that import data should try to
reproject it. This is getting too complicated and then there would be too
many modules to try to keep up to date. The conventional GRASS way seems
to be to import the data into a location with the correct projection, then
reproject it to your target location using [rsv].proj.

I meant that an import program like r.in.gdal should be able to associate
particular datum shift transformation information with the dataset. That is
if it knows that a non-default datum shift should be used it could populate
the PROJ_INFO file with it. This would seem to be very easy with your
approach.

Some points I would appreciate if people could comment on:
As far as I can make out a datum has a one-to-many relationship with datum
transformations, and one-to-one with ellipsoid and prime meridian. Neither
the datum lists in GRASS nor PROJ.4 include space for an associated prime
meridian so I'm wondering is this correct. I would like to make the GRASS
format as correct as possible. Probably I should investigate how PROJ.4
uses prime meridians but I don't really have that much spare time at the
minute and would appreciate a hint.

I don't see the prime meridian as an attribute of the datum, though it is
an attribute of a geographic coordinate system based on the datum. That is,
I think the prime meridian handling is orthagonal to the datum shift info.

That info should be in the PROJ_INFO file, IMHO. Anyway, I thought
we could have a Coordinate_System directory that g.setproj could
use. So far, with my state plane files, I have something like:

Coordinate_Systems/
  Projected/
     North_American/
        State_Plane_1927/
           State_Plane_Alabama_East
           State_Plane_Alabama_West
           ...
           State_Plane_Wyoming_West_Central
        State_Plane_1983/
           ...

This kind of thing is easy enough to add to, and g.setproj could
just use directory traversals until the user selects a file, then
just copy the file to the new location...

This sounds like a good idea but what are directory traversals and how
should they be implemented in g.setproj? There should be a GUI version as
well to make it look user-friendly and nice as g.setproj is a lot of
people's first encounter with a GRASS command (mine anyway and I spent
days over it) and it is really a terrible piece of software and very
off-putting. Is it all right to use Tcl/Tk in GRASS 5.1 or will there be
a different GUI?

I agree that some degree of categorization of coordinate system definitions
would be desirable. In many systems I have looked at there is really just
one level of categorization, but each pre-defined coordinate system can be
part of more than one group. So the groups might include:

   Common
   NAD27 State Plane
   NAD83 State Plane
   European
   North American
   EPSG Projected Coordinate Systems
   EPSG Geographic Coordinate Systems

and so on.

Some coordinate systems might appear in several groupings. If something like
this is built, I would ideally like to see it implemented within PROJ.4, so
that other systems could easily take advantage of it. If done within PROJ.4,
one approach would be to use each "initialization" file as a group, and add
a mechanism to supply some of the auxilary information with each coordinate
system. The existing EPSG, and state plane files would be a good start.

I am cc:ing this to the PROJ.4 list for possible additional comment.

Best regards,

--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

15

Hello Frank

On Wed, 12 Mar 2003, Frank Warmerdam wrote:

Paul Kelly wrote:
>
> So after selection, the datum transformation parameters will be stored in
> the PROJ_INFO file under the new key 'datumparams'.

Do correctly understand then that each named datum in the datumtransform
table will have a *default* transformation to use with the datum, but that
by explicitly including the transformation in the PROJ_INFO file it is
intended to be easy for users to alter what transformation values are used?

Yes I intended that if, for example, additional datum transform keywords
(other than towgs84 or nadgrids) are added to PROJ.4 in the future only
the PROJ_INFO file will need to be modified to take advantage of this.

I suppose it is implicit that there is a default transformation (as
defined by dx= dy= dz= in datum.table); this is really just for legacy
compatibility with the way datums are currently implemented in GRASS.
Both the current g.setproj and my changed version explicitly write the
transformation parameters into the PROJ_INFO file and these are what are
used. So only in the unlikely situation of the datum keyword being there
on its own with no parameters will the default transformation be used.

>>I do think that import programs that report projection information, or
>>offer to setup new locations (like r.in.gdal) should offer the transformation
>>used if possible.
>
> I don't really think that programs that import data should try to
> reproject it. This is getting too complicated and then there would be too
> many modules to try to keep up to date. The conventional GRASS way seems
> to be to import the data into a location with the correct projection, then
> reproject it to your target location using [rsv].proj.

I meant that an import program like r.in.gdal should be able to associate
particular datum shift transformation information with the dataset. That is

Yes, sorry didn't read your comment there properly at all. r.in.gdal does
conform to that convention exactly. I was actually thinking more of
v.in.gshhs which I was looking at earlier today. It offers to reproject
from lat/long into the location it is being imported into, but it is
complicated by the fact that GSHHS data doesn't have an associated ellipsoid
or datum. The program must decide which ellipsoid to associate the
lat/long data with. (Obvious choice: same one as the location the data is
being imported into, but still it should be the user making that decision
and not hidden inside the program.)

if it knows that a non-default datum shift should be used it could populate
the PROJ_INFO file with it. This would seem to be very easy with your
approach.

Yes definitely that would be easily possible.

I don't see the prime meridian as an attribute of the datum, though it is
an attribute of a geographic coordinate system based on the datum. That is,
I think the prime meridian handling is orthagonal to the datum shift info.

Well it looks like that can wait until co-ordinate systems are implemented
then. Good.

Thank you for the feedback

Paul

16

A little bit more of this is done now...

On Mon, 26 Aug 2002, Eric G. Miller wrote:

I've been looking at some of the projection related code, and I'm
trying to come up with a course of action to clean it up in places.
Eventually, this should result in:

1. Having better defined projection files.
2. Creation of a set of predefined well known coordinate systems
    files.
3. A rewrite of g.setproj with these ideas in mind:

Not re-written yet, just fixed up a bit. Re-write with a GUI can go in
GRASS 5.1

    - Types of coordinate system are XY, geographic, or projected
    - Ask users if they want to use a predefined coordinate
       system (state plane, UTM, others...) or if they want to
       define their own. Having a collection of predefined
       coordinate systems where the user can just pick it by
       "name", will simplify user's lives as well as insuring the
       info is correctly specified.
    - Both projected and geographic require datum/ellipsoid info.
       Ask for datum first, ellipsoid parameters can be looked up
       from it. Otherwise, just get ellipsoid parameters.

This is now mostly sorted out in the correct order

    - Remove state plane from projection list, it'll be covered
       by the predefined systems above (though still need units...)
    - Allow user defined datums??

This is possible; the PROJ.4 datum string (e.g. "towgs84=...") needs to be
known

    - Other...?
4. Have transparent handling for a variety of datum transformations
    incorporated into the library wrappers for proj4. Most of the
    necessary work for the transformations has been done by Andreas
    Lange, it's a matter of working in the strategies. Perhaps the
    programs can pass on a preferred stategy when specified by the
    user (for instance, NADCON vs. Bursa Wolf).

Anything possible with a parameter in the string passed to PROJ.4 can now
be used by including it in the 'datumparams: ...' line in the PROJ_INFO
file (or specify a custom datum in g.setproj). g.setproj can be re-run to
change the transformation parameters, though it would be wise to make a
copy of the current projection parameters by running g.projinfo first;
probably easier just to edit the PROJ_INFO file directly if you know the
parameters you want.

    Get rid of GRASS
    copy of proj4 and make a dependency on proj4 >= whatever version.

Today I have done this (well not get rid of the GRASS copy; I just made it
optional to use an external PROJ.4 library). This can be invoked by using
the --with-proj=yes option to the configure script and the usual
--with-proj-includes=... --with-proj-libs=... also apply if needed. It
needs to be PROJ.4.4.6 or later as the earlier versions had some bugs in the
datum transformations.

The changes required to the modules that use proj were:
1. Include gprojects.h instead of projects.h
2. Include $(PROJINC) in the EXTRA_CFLAGS in the Gmakefile

I had to split the projects.h up into PROJ.4 internal stuff and bits
specific to the GRASS proj wrapper functions, and rename it so it wouldn't
conflict with an externally-installed copy of PROJ.4.

As far as I can see, under some circumstances a 'make clean' might be
needed, or maybe just delete the old libproj.a in src/libes/LIB.arch. The
library containing the GRASS proj wrapper functions is now called libgproj
instead of libproj, again so as not to conflict with the external version.

    It should be a fatal error if a datum transformation is impossible.
5. Miscellaneous:
    - Specify vertical datums? MSL, ellipsoid height, others... Mostly
       for documentation, since very few can be translated.
    - Add missing usfoot/usfeet <-> metric conversions. usinch? usmile?
    - Add prime meridian as datum parameter.

Frank said he didn't see this as an attribute of the datum, but of the
co-ordinate system, so I have left it for now.

Coordinate_Systems/
   Projected/
      North_American/
         State_Plane_1927/
      State_Plane_Alabama_East
      State_Plane_Alabama_West
      ...
      State_Plane_Wyoming_West_Central
   State_Plane_1983/
      ...

This kind of thing is easy enough to add to, and g.setproj could
just use directory traversals until the user selects a file, then
just copy the file to the new location...

Frank thought it would be a good idea to put this into PROJ.4, also I
don't know how to do directory traversals, so I will forget about this
for now.

One point that will need to be emphasised (for American users I suppose)
when this work is incorporated into a release of GRASS: The last few releases
(5.0.0pre4 to 5.0.1 approx.) have, for the specific case of transforming NAD27
to NAD83 and vice versa, ignored the dx dy dz-style transformation parameters
given in the PROJ_INFO file and used the nadcon conversion tables. Under
the new 'improvements' the 3-parameter shift will be used if it is
specified in the PROJ_INFO file, and g.setproj will need to be re-run on
the location (or the PROJ_INFO file edited manually) to remove the dx dy
dz and add a 'datumparams: nadgrids=conus' line.

Hoping for some testing and feedback,

Paul

17

Paul Kelly wrote:

> Get rid of GRASS
> copy of proj4 and make a dependency on proj4 >= whatever version.

Today I have done this (well not get rid of the GRASS copy; I just made it
optional to use an external PROJ.4 library). This can be invoked by using
the --with-proj=yes option to the configure script and the usual
--with-proj-includes=... --with-proj-libs=... also apply if needed. It
needs to be PROJ.4.4.6 or later as the earlier versions had some bugs in the
datum transformations.

As soon as this is known to work, the private copy of PROJ should be
removed, IMHO.

BTW, what are $(GRASSINT) and $(PROJINT) for? And where do they get
their values?

--
Glynn Clements <glynn.clements@virgin.net>

18

On Fri, 11 Apr 2003, Glynn Clements wrote:

BTW, what are $(GRASSINT) and $(PROJINT) for? And where do they get
their values?

These are defined in the Gmakefile for the GRASS proj library
src/libes/proj/Gmakefile. $(PROJINT) includes all the object files from
PROJ.4 and $(GRASSINT) contains get_proj.o and do_proj.o, the two files
that contain the GRASS wrapper functions for PROJ.4. So by default
$(GRASSINT) and $(PROJINT) are both compiled, i.e. the grass libgproj.a
contains all the PROJ.4 functions and the GRASS wrapper functions. If
--with-proj=yes is specified then only $(GRASSINT) is compiled, i.e.
libgproj.a contains only the wrapper functions and is only about 40-50kB
in size.

I had never used autoconf before but after a bit of thinking it seemed
an elegant enough way of doing it. In GRASS 5.1 I have made an external
PROJ.4 an absolute dependency so this bit isn't needed there. It seems
there is quite a consensus to remove the local copy of PROJ.4 for GRASS
5.0 as well, but I will still wait a while and hope for further discussion
before doing this. The required PROJ.4.4.6 was only released recently so a
lot of people may have an older version installed, which will have had
a small bug resulting in errors at the 10s of centimetres level when doing
7-parameter datum transformations (this is now fixed in the local GRASS
copy but I suppose it is possible there are other bugs there nobody has
noticed).

An aside: the Gmakefile / Makefile for GRASS libgproj also converts the
ASCII nadcon datum shift files into binary using the nad2bin program from
the PROJ.4 distribution. So this program must be installed in the path,
e.g. /usr/local/bin or somewhere for it to work. I couldn't see any
obvious convention for specifying the path to the nad2bin binary in the
configure script, but it should be fine as long as PROJ.4 is properly
installed on the system.

Paul

19

Paul Kelly wrote:

> BTW, what are $(GRASSINT) and $(PROJINT) for? And where do they get
> their values?

These are defined in the Gmakefile for the GRASS proj library
src/libes/proj/Gmakefile. $(PROJINT) includes all the object files from
PROJ.4 and $(GRASSINT) contains get_proj.o and do_proj.o, the two files
that contain the GRASS wrapper functions for PROJ.4. So by default
$(GRASSINT) and $(PROJINT) are both compiled, i.e. the grass libgproj.a
contains all the PROJ.4 functions and the GRASS wrapper functions. If
--with-proj=yes is specified then only $(GRASSINT) is compiled, i.e.
libgproj.a contains only the wrapper functions and is only about 40-50kB
in size.

OK. You need to use single quotes; bash interprets $(...) in the same
way as `...` (command substitution), resulting in errors from
configure:

/usr/src/grass/configure: PROJINT: command not found
/usr/src/grass/configure: GRASSINT: command not found

and GPROJLIBOBJS being empty.

An aside: the Gmakefile / Makefile for GRASS libgproj also converts the
ASCII nadcon datum shift files into binary using the nad2bin program from
the PROJ.4 distribution. So this program must be installed in the path,
e.g. /usr/local/bin or somewhere for it to work. I couldn't see any
obvious convention for specifying the path to the nad2bin binary in the
configure script, but it should be fine as long as PROJ.4 is properly
installed on the system.

The usual approach is to just run the program by name, and assume that
it's in the path. If you want to use a full path, the autoconf macro
AC_PATH_PROG can be used.

However, there's an issue regarding the datum files which we should
really deal with at some point. Compiling the datum files during the
build process doesn't work when cross-compiling, as the generated
etc/nad2bin program will be compiled for the target platform rather
than the build platform, so it won't run on the build host.

Using the installed nad2bin should solve that particular problem, but
create another: when we last discussed this, I was told that the
compiled datum files weren't portable so, while the installed nad2bin
would run, the resulting files may not work on the target platform.

The only reliable solution (unless the PROJ maintainers can be
persuaded to make the binary format portable) is to include a suitable
nad2bin program in the final package, and compile the datum files on
the target platform as part of the installation process.

BTW, a similar problem exists for the display fonts in src/fonts.

--
Glynn Clements <glynn.clements@virgin.net>

20

I have the following errors building the current 5.1 code on Solaris 8 sparc

configure:
The mysql checks fail.
Line 9519 - add "-lsocket"
Line 9559 - add "-lsocket -lm"

The prototype for G_ask_datum_name in gisdefs.h is different from the
function in get_datum_name.c

/lib/gis/get_datum_name.c has: "int G_ask_datum_name(char *datum)"

Change /include/gisdefs.h line 520
from: int G_ask_datum_name(char *datum, char *datum)
to: int G_ask_datum_name(char *datum)

This one I have not figured out:
make is looking for the non-existent file "datumtransform.table" in
/lib/gis. There is a "datum.table", but no "datumtransform.table" Where
should I look?