[GRASS5] More on Projections

GRASS GRASS developer
1

GRASS Luminaries,

A little more I need to know about projections:

o If the proj number in a raster header, or the WIND/DEFAULT_WIND file is
   99 (PROJECTION_OTHER) does that mean that the definition of the projection
   can only be found in the PROJ_INFO file, and that it will be one of the
   types other than state plane, utm, latlong or XY?

o I was looking at the r.proj command to see how it gets the projection
   definition of the source raster. It appears to depend on the PROJ_INFO
   file existing (via a call to G_get_projinfo()) even though this file
   does not seem to exist on most of the sample datasets on the GRASS Sample
   Dataset page. It seems to me that r.proj should fallback to looking at
   the proj: code of the raster, or the projection definition in the WIND
   file, but as far as I can tell it does not.

   How should a GRASS program safely determine what the projection of a raster
   is?

For my purposes, it is convenient to have a conventional PROJ.4 string as a
way of expressing the complete coordinate system of a raster, including
units conversions. For this purpose, I have implemented the following
functions:

/***************************************************************************
* char *G_get_cell_as_proj4( cell, mapset );
*
* char *cell Name of cell to query.
* char *mapset Name of mapset to query.
*
* Returns a PROJ.4 definition of the coordinate system. The returned
* string should be freed with G_free() when no longer needed. XY coordinate
* system, or other errors during processing will return in a return value
* of NULL.
**************************************************************************/

/***************************************************************************
* char *G_kv_to_proj4( in_proj_keys, in_units_keys );
*
* struct Key_Value *in_proj_keys
* List of projection keys, as returned by
* G_get_projinfo().
* struct Key_Value *in_units_keys
* List of units keys, as returned by
* G_get_projunits().
*
* Returns a PROJ.4 definition of the coordinate system. The returned
* string should be free with G_free() when no longer needed. Note that
* this function will return a PROJ.4.4.2 string, which includes support
* for the +proj=longlat.
**************************************************************************/

In retrospect, it seems having G_get_cell_as_proj4() is unnecessary, and that
I should instead have written a G_get_location_as_proj4() since everyone seems
to agree that the projection definition should be identical in all the
rasters in a location. Would it be useful to implement something like:

int G_get_location_coordsys( char *location /* can be NULL for current */,
                          Key_Value **proj_info, Key_Value **unit_info );

that would return the projection, and unit definition for a location (in
proj_info, and unit_info) in a form ready to use with pj_get_kv()? This
function would use the PROJ_INFO and PROJ_UNITS if available, otherwise it
would fallback to the WIND file? Then programs like r.proj could be updated
to use this function?

Finally, I am still in great need of examples of locations with "interesting"
PROJ_INFO and PROJ_UNITS files.

Best regards,

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

----------------------------------------
If you want to unsubscribe from GRASS Development Team mailing list write to:
minordomo@geog.uni-hannover.de with
subject 'unsubscribe grass5'

2

On Thu, 21 Sep 2000, Frank Warmerdam wrote:

A little more I need to know about projections:

o If the proj number in a raster header, or the WIND/DEFAULT_WIND file is
   99 (PROJECTION_OTHER) does that mean that the definition of the projection
   can only be found in the PROJ_INFO file, and that it will be one of the
   types other than state plane, utm, latlong or XY?

Definition of projection is only found in the PROJ_INFO file. I assume
that the number (1, 2, 3, 99) at the top of the cellhd file (and WIND,
DEFAULT_WIND) originally was a way to identify the projection of the map
or region, but as more projections were added to grass later they all got
the the number 99. As a consequence you cannot see what projection a 99
map is made in. If you copy a map and its associated cellhd file to
another location in a different projection, you may not get an error when
trying to display or process it, but the result will certainly not be what
you expect.

o I was looking at the r.proj command to see how it gets the projection
   definition of the source raster. It appears to depend on the PROJ_INFO
   file existing (via a call to G_get_projinfo()) even though this file
   does not seem to exist on most of the sample datasets on the GRASS Sample
   Dataset page.

I don't know this sample data, but it's probably in
XY-(non)projection. XY-locations do not have a PROJ_INFO file.

It seems to me that r.proj should fallback to looking at
   the proj: code of the raster, or the projection definition in the WIND
   file, but as far as I can tell it does not.

The WIND file does not contain projection info, it contains information on
the borders and the resolution. r.proj (and r.proj.new) uses
G_get_projinfo() to read in the parameteres from PROJ_INFO into a
memory structure for later use.

   How should a GRASS program safely determine what the projection of a raster
   is?

Something like
myprojkeys = G_get_projinfo();
myprojname = G_find_key_value("proj", myprojkeys);

will give you the name of the projection of the current location. Other
parameters should be extracted too of course. See
src/libes/proj/get_proj.c for examples.

[snip]

Currently I prepare a struct Cell_head, attempting to set the projection
to match that of the source file if possible.

I don't it's a good idea to try to change the projection of an existing
location (even though it's possible, e.g with g.setproj) because the
user's old data will be rendered useless if the original projection cannot
be restored afterwards. Or do you mean setting the projection of the just
output map? In that case, the approach is wrong, because in GRASS the
projection is a property of the location, which all maps in it will
inherit. The only place where you find that info is in PROJ_INFO.
Individual maps do not have separate projection definitions.

Am I responsible for doing validation to ensure the projection of
my new raster matches the location?

Yes, unless you instruct the user that a new location matching the source
map must be created manually before attempting to import the map.

If the coordinate system doesn't match the existing location I would
like it easy for the user to create a new location with a coordinate
system matching the source file. Is this something that raster import
programs ever do? Can you point me to an example that does this?

Good idea, but it may be hard to accomplish, because there are so many
paramenters to check, and all of them must somehow be included in the
source file, and your program should be able to parse them out.

Finally, as a user how do I change my current location? Do I need to
exit out of grass and restart it specifying a new location?

Yes

regards
Morten Hulden

----------------------------------------
If you want to unsubscribe from GRASS Development Team mailing list write to:
minordomo@geog.uni-hannover.de with
subject 'unsubscribe grass5'

3

Morten Hulden wrote:

> o I was looking at the r.proj command to see how it gets the projection
> definition of the source raster. It appears to depend on the PROJ_INFO
> file existing (via a call to G_get_projinfo()) even though this file
> does not seem to exist on most of the sample datasets on the GRASS Sample
> Dataset page.

I don't know this sample data, but it's probably in
XY-(non)projection. XY-locations do not have a PROJ_INFO file.

Morten,

In particular I have downloaded the following datasets from
http://www.geog.uni-hannover.de/grass/data.html

o GRASS 5.x Spearfish: No PROJ_INFO, everything is in UTM.

o IMAGERY dataset: No PROJ_INFO, but everything is projection 0, unreferenced.

o ETOPO5 elevation mode: Has a PROJ_INFO, though it's just in lat/long.

o Global Dataset: No PROJ_INFO, everything is in lat/long.

It seems to me that it is traditional for GRASS programs to work, even with
relatively old grass databases. Based on this assumption I think there needs
to be transparent support for datasets without a PROJ_INFO, but that do have
a well defined projection.

Rather than adding the g_get_location_coordsys() I suggested before, I now
believe
we should extend G_get_projinfo() to create a key/value list based on the
current
window setting in the WIND file, if no PROJ_INFO file exists. This would mean
that applications can safely just call G_get_projinfo() without having to worry
whether the dataset is a modern one with a PROJ_INFO file, or an older one
without.
I would be happy to make this change, if folks support it.

The WIND file does not contain projection info, it contains information on
the borders and the resolution. r.proj (and r.proj.new) uses
G_get_projinfo() to read in the parameteres from PROJ_INFO into a
memory structure for later use.

Well, it does contain projection information via the proj: keyword, but you
are correct that it can't define the full set of projections.

> How should a GRASS program safely determine what the projection of a raster
> is?

Something like
myprojkeys = G_get_projinfo();
myprojname = G_find_key_value("proj", myprojkeys);

will give you the name of the projection of the current location. Other
parameters should be extracted too of course. See
src/libes/proj/get_proj.c for examples.

On odd thing about get_proj() is this stuff about looking at the name keyword.
The code seems to be organized as if the proj keyword doesn't always exist
and if it isn't it will check to see if name start with "Lat". Is this
generally necessary?

[snip]

> Currently I prepare a struct Cell_head, attempting to set the projection
> to match that of the source file if possible.

I don't it's a good idea to try to change the projection of an existing
location (even though it's possible, e.g with g.setproj) because the
user's old data will be rendered useless if the original projection cannot
be restored afterwards.

The G_set_window() does not update the WIND file, but does affect the working
window of the current process, and affects what gets written out to the cellhd
of newly created rasters. It is necessary to do this get the right extents,
but at the same time it is necessary to set something for the projection, and
zone to put into the Cell_head that will be used to set the current window.

In retrospect, I think I should always just set those by calling G_projection()
and G_zone() to get the projection and zone of the current projection (which
normally comes from the WIND file). Checking whether the file r.in.gdal is
reading actually matches the current location should be a separation activity.

Or do you mean setting the projection of the just
output map? In that case, the approach is wrong, because in GRASS the
projection is a property of the location, which all maps in it will
inherit. The only place where you find that info is in PROJ_INFO.
Individual maps do not have separate projection definitions.

Is there any utility that will check a location or database for validity?
Check stuff like all the projection information in a location making sense?

> Am I responsible for doing validation to ensure the projection of
> my new raster matches the location?

Yes, unless you instruct the user that a new location matching the source
map must be created manually before attempting to import the map.

I really don't feel it is appropriate to make the user check that the
projections check if it is possible for the program to do so. I suppose
there might be cases where a user should be able to override what appear
to be differences in coordinate system when they feel they are unimportant
(such as slightly different ellipsoids on data at a scale where that doesn't
matter).

> If the coordinate system doesn't match the existing location I would
> like it easy for the user to create a new location with a coordinate
> system matching the source file. Is this something that raster import
> programs ever do? Can you point me to an example that does this?

Good idea, but it may be hard to accomplish, because there are so many
paramenters to check, and all of them must somehow be included in the
source file, and your program should be able to parse them out.

Well, I think it behooves us to find a way to make this easy. I firmly
believe that GRASS needs to provide better library facilities making it
easy for application programmers to do projection validation and to import
data with proper coordinate systems. For instance, in a translator there
should be an easy way for me to return a Key_Value list defining the
coordinate system of a file, and verify that it matches the current location.
If it doesn't I should have an easy way for the user create a new location
based on that information.

Should I propose an API for doing this, and implement it once it is agreed
upon?

To summarize I want to do the following this:

o Upgrade G_get_projinfo() to "generate" appropriate information if only
   a WIND file is available, and there is no PROJ_INFO.

o I want to propose some functions that applications can use to easily
   compare their datasources projections with the current PROJ_INFO to ensure
   they are consistent.

o I want to initially deploy use of these functions in r.in.gdal, but I hope
   they will be retrofit to other coordinate system aware importers in the
   future.

Does anyone object to make doing these things? Note the only portion of this
that risks introducing instability is the G_get_projinfo() change, but I can't
really see much risk assuming it is done properly.

Best regards,

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

----------------------------------------
If you want to unsubscribe from GRASS Development Team mailing list write to:
minordomo@geog.uni-hannover.de with
subject 'unsubscribe grass5'