[PATCH] New projections for grass5.0beta2

Sun Aug 1 18:52:30 EDT 1999

Hello list

I needed the Mollweide and Goode projections in GRASS5.0beta2 for some
global maps I am working on, so I added them. It was not too much work,
because the PROJ4 package that everything is based on is already part of
GRASS, but has never been fully implemented in the g.setproj, r.proj,
v.proj or m.proj modules. I added support for the rest of included
projections as well, because it was not that different. 

The week-end is over and I have to go back to normal work. But I would
like to release my (incomplete) patch to the interested now, because I do
not have time to work with this on a regular basis and don't know when/if
I will be able to continue. I hope releasing it will help encourange
testing and clarify some points on the work that has to be done to include
full support for all PROJ4 projections in a future GRASS version. There
are still a lot of gotchas in sanity checks and making things to work
flawlesly with other GRASS modules.

The file patch-grass5.0beta2-mh1 does the following:

	-enables creating and editing locations in any of the projections of the
	proj version that is currently included in grass5.0beta2. With this 
	the total number of projections supported by GRASS is 76.

	-enables reprojection of rasters and vectors with v.proj and r.proj 
	into the new projections.

WARNING AND DISCLAIMER

	This is experimental, incomplete, unofficial and probably buggy code. 
	I do not take any responsibility for what may happen to your data (or 
	even to your hardware :v) if you apply this patch. But if you know what 
	you are doing you will get additional projections to play around with 
	in GRASS. But be prepared to revert to the official version at its next 
	release. This patch should only be applied to a clean grass5.0beta2 source 
	tree. (But changes you have made to files not mentioned in the CHANGES 
	section does not matter.)

WHERE TO GET IT

	You can download the patch from 
	http://www.ngb.se/~morten/grass/patch-grass5.0beta2-mh1.tar.gz

HOW TO APPLY

	1) Untar the patch-grass5.0beta2-mh1.tar.gz file into the top level of 
		the source directory of grass5.0beta2. 
	2) Run the command
		patch -p1 < patch-grass5.0beta2-mh1 
	3) Recompile the sources.
		'make install' after removing the file in
		src/CMD/next_step  
		(You don't have to clean your source tree. 'Make' will know
		which files need to be recompiled.)

CHANGES

	The following files have been changed:

	src/libes/gis/projection.h
		Extended list of supported projections
	src/general/g.setproj/local_proto.h
		Added prototypes for new functions in get_num.c 
	src/general/g.setproj/get_num.c
		Added new functions to get values for new options
	src/general/g.setproj/table.h
		Added headers for new projections and options 
	src/general/g.setproj/table.c
		Added new projection and parameter default values
	src/general/g.setproj/main.c
		Updated code to build PROJ_INFO and added some sanity 
		checks for some of the new projections.

	I have tried to change as little as possible, though some
	of the code looks perculiar to me and I would probably have
	done things differently. (But what do I know?)

LIMITATIONS AND BUGS

	The input map for r.proj _must_ be in a 'unprojected' (LL) or cylindical 
	projection. The reason is that only in a cylindical projection can 
	the N,S,E and W values from cellhd be combined and used directly as 
	corner coorinates (NW NE SE SW), which is what r.proj does. Of course,
	trying to project rasters from a spherical projection map in nonsense, 
	because the corner co-ordinates of the rectangular area may not be part of 
	the map at all. This does not affect v.proj; v.proj can project vectors using 
	any projection as input.	

	The 'north/east must be larger than south/west' restriction prevents
	projections and display of partial maps in some areas of conical and
	azimuthal projections. It seems to be a very basic assumpion in 
	GRASS, so I have not even though seriously about by-passing it.
	But it means that only one half of a 360-degree polar cap map can be 
	projected into or displayed from. You _can_ project and display a 
	complete map (covering the whole default region), but if the 
	resolution is too high may not have enough memory to do that 
	projection. r.proj wants to read the whole input map into memory 
	becore projecting.

	Sanity checks and default values of some of the new projections may 
	be inclomplete or incorrect. I don't know enough about projections to 
	be able to provide reasonable defaults or checks for all cases.
	(Since PROJ modules _do_ error checking and give diagnosis back,
	it would be a better solution to just pass input through to a test
	exec and have setup fail if PROJ complains.)

	The omerc projection needs special handling of input because of the
	alternative ways to use parameters. Haven't finished this. 

	Some projections have fixed parameters (alsk, gs50). g.setproj will not
	ask the user for input for these but will write the fixed parameters directly
	into PROJ_INFO. PROJ_INFO should always	reflect these values and should never 
	be changed manually. (The PROJ routines don't care but other GRASS
	modules may fail). 

	In global spherical projections v.proj produces vector maps with artifact 
	lines accoss the map, connecting (i guess) opposite points on the perimeter, 
	which are actually part of the same line or polygon.

	r.proj sometimes produces empty maps with projecting rasters into
	some of the spherical projections. Don't know why.

	---------------------------------------------------------------------------

	I found two BUGS in the original g.setproj:

	In src/general/g.setproj/main.c line 180
	          buf = G_find_key_value("radius", out_proj_keys);
	should read
	          buf = G_find_key_value("a", out_proj_keys);

	otherwise, if the ellipsoid is a 'sphere', the radius 'a' will be
	reset to 0.0000000000 if g.setproj is run a second time. This bug
	is fixed in this patch, but if you do not want to apply the patch
	you may at least want to change line 180 manually.

	The second bug is that old projection parameters will not be cleared
	from PROJ_INFO if the projection is changed with g.setproj, and the
	old projection uses parameters that are not used by the new one. I 
	haven't fixed that. I doubt the usefulness of allowing users to
	change the projections of existing locations anyway, because their 
	data will become unusable.

NOTES
	Maximum cartesian N/S and E/W values are of course calculated differently for 
	different projections. For example, Orthographic global earth has 
	N/S/E/W_max=earth radius, while Sinusoidal has N/S_max = earth perimeter/4 
	and E/W_max = eart perimeter. I use my externally installed proj program 
	to calculate default region values in meters before setting up a new location.  

TO DO

	It was my intention to update m.proj as well to support the new
	projections, but I did not have the time. Code and headers are
	very similar to g.setproj, in fact, at least the headers and 
	table.c could be merged. Maybe I will give it a try later.
	But since I have the PROJ4.3.3 package installed externally I have
	had no need for m.proj.

	Except for trying to get around the limitations mentioned above, there
	are some 45 projections in the PROJ4.3.3 package that are not yet
	inlcuded in the GRASS sources. But I hesitate to get into this
	deeper than necessary for my own purposes. 

AUTHOR
	Morten Hulden
	on linux 2.2.9, libc6, patched RH6 and Mandrake, i686