[GRASS5] More on Projections

[Frank Warmerdam]

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 at 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 at geog.uni-hannover.de with
subject 'unsubscribe grass5'