[GRASS5] datum transforms...
Paul Kelly
paul-grass at stjohnspoint.co.uk
Wed Jan 15 08:49:53 EST 2003
Hello everyone
Well I've had a bit more time to look at this now and with the
suggestions from Eric and Frank I think I've come up with something
usable. See below for my comments on what I've changed and the attached
diff file for the changes.
grass/src/libes/gis/datum.table
grass/src/libes/proj/pj_datums.c
grass/src/libes/proj/pj_ellps.c
Add support for Ireland 65 Datum, correct missing underscores in
ellipsoid and datum names. Add support for sphere in pj_ellps.c
grass/src/libes/proj/do_proj.c
Add two new functions written by Eric Miller, pj_do_proj, a direct
replacement for the existing pj_do_proj which uses the PROJ.4
pj_transform() function instead of pj_fwd() and pj_inv(), which takes
advantage of datum transformation. Second new function pj_do_transform()
is for future use and allows a height value to be tranformed as well
(from ellipsoidal to orthometric height?) and arrays of co-ordinates to
be transformed in one pass.
grass/src/libes/proj/get_proj.c
get_proj.c contains three functions for setting up co-ordinate systems:
pj_get_kv() (which uses PROJ.4 parameters read from a PROJ_INFO file),
pj_get_string() (which uses PROJ.4 parameters passed to it in a string)
and pj_zero_proj() (which initialises a pj_info struct to 'no data').
'latlong' is now treated the same as a projected co-ordinate system and
a parameter 'proj=latlong' must be passed to PROJ.4. pj_get_kv and
pj_get_string have been changed to do this.
The parameter-passing logic in pj_get_kv has been changed to treat
pj_ellps.c and pj_datums.c as the definitive source of datum and
ellipsoid parameters, i.e. ellipsoid parameters are only passed on if
there is no ellipsoid name, and ellipsoid name is only passed on if
there is no datum name. When using pj_get_string() the programmer can
override this behaviour.
To use pj_transform, there must be an instance of the PJ struct (which
holds the parameters passed to PROJ.4) within the pj_info struct (which
GRASS programs use to communicate with these wrapper functions). At its
simplest, the PJ struct will contain one parameter, 'proj=latlong', but
this still must be set up now.
A PJ struct is created when pj_init is called. The first two functions
(pj_get_kv and pj_get_string) call pj_init and create a PJ struct within
the pj_info struct; the last one (pj_zero_proj) doesn't so it should not
be used any more.
Some programs set up a latlong projection system as follows:
pj_zero_proj(&oproj);
sprintf(oproj.proj,"%s","ll");
This won't work anymore as pj_zero_proj() does not create a PJ struct;
some programs used the following syntax which will work:
pj_zero_proj(&info_in);
parms_in[0] = '\0';
pj_get_string(&info_in, parms_in);
as pj_get_string creates calls pj_init and creates a PJ struct (passing
a null string creates a latlong projection). However the call to
pj_zero_proj() is superfluous as the 3 lines in that function are also
contained in pj_get_string. So I removed that call from these programs.
And the other programs that used the first method to set up a basic
latlong projection system have been changed to use the second method
(see details below). pj_get_string can also take NULL for the second
argument; is this the same as a string with its first character '\0'? As
I wasn't able to test most of these changes I kept the syntax exactly
the same as in other GRASS programs so that there should be no problems.
I would suggest removing pj_zero_proj() from the GRASS API now as it
doesn't seem to have any valid purpose and I removed all calls to it in
the user programs.
Also in get_proj.c, a new function added:
set_proj_lib() is a finder function for use with the PROJ.4 function
pj_set_finder(), and sets the path name where PROJ.4 should look for the
datum conversion tables. pj_set_finder is called with this function as
an argument, in pj_get_kv and pj_get_string, before pj_init is called.
(Deleted:
grass/src/libes/proj/datum_shift.c
grass/src/libes/proj/do_proj_nad.c)
These contained the set_datumshift() and do_proj_nad() functions which
were only necessary for the previous implementation of NAD/NAD datum
transformation, which used lower level PROJ.4 functions than
pj_transform().
grass/src/include/projects.h
Add prototypes for new functions and remove deleted ones
grass/src/libes/proj/Gmakefile
Remove datum_shift.o and do_proj_nad.o
Other GRASS programs which use PROJ.4 library:
1. grass/src/general/g.region/cmd/printwindow.c
Change from
pj_zero_proj(&oproj);
sprintf(oproj.proj,"%s","ll");
style initialisation to
parms_in[0] = '\0';
pj_get_string(&info_in, parms_in);
style.
2. grass/src/mapdev/v.in.gshhs/main.c
No problems; already used pj_get_string as above. Removed superfluous
call to pj_zero_proj().
3. grass/src/mapdev/v.proj/local_proto.h
grass/src/mapdev/v.proj/main.c
Revert program from using separate NAD/NAD datum conversion
4. grass/src/misc/m.ll2db/main.c
No problems; already used pj_get_string. Removed superfluous call to
pj_zero_proj().
5. grass/src/raster/r.proj/cmd/bordwalk.c
grass/src/raster/r.proj/cmd/main.c
Revert program from using separate NAD/NAD datum conversion
6. grass/src/raster/r.sun/main.c
Same as 1.
7. grass/src/raster/r.sunmask/g_solposition.c
Same as 1.
8. grass/src/sites/s.datum.shift/main.c
Same as 1.
9. grass/src/sites/s.proj/main.c
Revert program from using separate NAD/NAD datum conversion
10. grass/src/mapdev/v.digit/map_init_new.c
No problems; already used pj_get_string.
11. grass/src/libes/image3/convert_ll.c
12. grass/src/mapdev/v.mkquads/convert.c
13. grass/src/raster/r.in.gdal/main.c
No problems; all already used pj_get_kv(). Only v.mkquads and r.in.gdal
appear to read the PROJ_INFO file for two separate locations, and will
thus be able to take advantage of datum conversion straight away, as
will [rsv].proj.
14. grass/src/misc/m.proj/main.c
15. grass/src/misc/m.proj2/main.c
inter and cmd versions of this program will have no problems with the
new pj_do_transform() but from what I can see it looks like they won't
be able to take advantage of datum transformation either as they
manually form the projection parameters for PROJ.4 and don't pass on the
datum, even if it existed in the PROJ_INFO file. Would need to be looked
into more though.
More general notes:
PROJ.4 will perform datum conversion if a datum= parameter is passed to
it, or a towgs84= or nadgrids= parameter, for BOTH the input AND output
projections. All the datum= parameter really does is expand the
corresponding definition in the pj_datums.c file, which includes the
correct towgs84= or nadgrids= value, together with the ellipsoid name.
If the PROJ_INFO file contains a datum: entry, and the program passes
all the entries in this file onto PROJ.4 for both projections (as
[rsv].proj do and r.in.gdal and v.mkquads appear to do, as far as I can
see), then datum conversion will automatically happen. If only one
projection has a datum specified then datum conversion won't take place.
For a datum that isn't included in the pj_datums.c list, you could
manually add the necessary parameters to your PROJ_INFO file. Other
parts of GRASS would probably ignore them but I can't be sure. E.g.
dx:
dy:
dz:
or
towgs84: dx,dy,dz
should do the same thing,
or
nadgrids: alaska,hawaii
The last example might work and would be useful if your location is one
of the NAD conversion areas not covered by the conus and ntv1_can.dat
files.
For more details on how the changes work look at the diff file.
If nobody has any objections I'd really like to have this put into the
CVS for testing. I wonder is the attached file enough for that; the CVS
server seems to have an access problem (I get a Connection refused error
upon 'cvs login') so I took the latest CVS snapshot from the website and
diffed against that. I hope it will do.
Paul
-------------- next part --------------
diff -ru grass50_exp_2003_01_10/src/general/g.region/cmd/printwindow.c grass/src/general/g.region/cmd/printwindow.c
--- grass50_exp_2003_01_10/src/general/g.region/cmd/printwindow.c 2002-05-17 04:12:34.000000000 +0100
+++ grass/src/general/g.region/cmd/printwindow.c 2003-01-15 12:39:32.712333000 +0000
@@ -89,8 +89,8 @@
G_fatal_error("Can't get projection key values of current location");
/* set output projection to lat/long */
- pj_zero_proj(&oproj);
- sprintf(oproj.proj, "%s", "ll");
+ buf[0]='\0';
+ pj_get_string(&oproj, buf);
/* do the transform
* syntax: pj_do_proj(outx, outy, in_info, out_info)
diff -ru grass50_exp_2003_01_10/src/include/projects.h grass/src/include/projects.h
--- grass50_exp_2003_01_10/src/include/projects.h 2002-04-26 04:12:31.000000000 +0100
+++ grass/src/include/projects.h 2003-01-15 13:02:06.345311000 +0000
@@ -365,14 +365,12 @@
/* do_proj.c */
int pj_do_proj PROTO((double *, double *, struct pj_info *, struct pj_info *));
-/* do_proj_nad.c */
-int pj_do_proj_nad PROTO((double *, double *, struct pj_info *, struct pj_info *));
+int pj_do_transform PROTO((int, double *, double *, double *,
+ struct pj_info *, struct pj_info *));
/* get_proj.c */
int pj_get_string PROTO((struct pj_info *, char *));
int pj_zero_proj PROTO((struct pj_info *));
-
-/* pointer for datum conversion */
-int (*proj_f)(double *, double *, struct pj_info *, struct pj_info *);
+const char * set_proj_lib PROTO((const char *name));
#ifdef GRASS_GIS_H
int pj_get_kv PROTO((struct pj_info *, struct Key_Value *, struct Key_Value *));
diff -ru grass50_exp_2003_01_10/src/libes/gis/datum.table grass/src/libes/gis/datum.table
--- grass50_exp_2003_01_10/src/libes/gis/datum.table 2002-09-13 04:12:41.000000000 +0100
+++ grass/src/libes/gis/datum.table 2003-01-15 12:32:37.720817000 +0000
@@ -63,3 +63,6 @@
carthage "Carthage 1934 Tunisia" clark80 dx=-263.0 dy=6.0 dz=431.0
# Hermannskogel
hermannskogel "Hermannskogel" bessel dx=653 dy=-212 dz=449
+# Ireland 1965 - From NIMA document - but more accurate 7-parameter transform
+# from Ordnance Survey Ireland is specified in src/libes/proj/pj_datums.c
+ire65 "Ireland 1965" modif_airy dx=506 dy=-122 dz=611
diff -ru grass50_exp_2003_01_10/src/libes/proj/Gmakefile grass/src/libes/proj/Gmakefile
--- grass50_exp_2003_01_10/src/libes/proj/Gmakefile 2002-05-31 04:12:41.000000000 +0100
+++ grass/src/libes/proj/Gmakefile 2003-01-15 12:37:39.086416000 +0000
@@ -51,7 +51,7 @@
RLIB = pj_release.o
-GRASSINT = get_proj.o do_proj.o do_proj_nad.o datum_shift.o
+GRASSINT = get_proj.o do_proj.o
NAD_DIR = $(GISBASE)/etc/nad
NAD2BIN = $(GISBASE)/etc/nad2bin
Only in grass50_exp_2003_01_10/src/libes/proj: datum_shift.c
diff -ru grass50_exp_2003_01_10/src/libes/proj/do_proj.c grass/src/libes/proj/do_proj.c
--- grass50_exp_2003_01_10/src/libes/proj/do_proj.c 2002-04-26 04:12:44.000000000 +0100
+++ grass/src/libes/proj/do_proj.c 2003-01-15 12:36:29.710319000 +0000
@@ -6,52 +6,116 @@
/********************* PARMS FOR PROJ **********************/
-#define STP_FT 0.30480060960121920243
+/* a couple defines to simplify reading the function */
+#define MULTIPLY_LOOP(x,y,c,m) \
+do {\
+ int i; \
+ for (i = 0; i < c; ++i) {\
+ x[i] *= m; \
+ y[i] *= m; \
+ }\
+} while (0)
+
+#define DIVIDE_LOOP(x,y,c,m) \
+do {\
+ int i; \
+ for (i = 0; i < c; ++i) {\
+ x[i] /= m; \
+ y[i] /= m; \
+ }\
+} while (0)
-static double STP = 1.0;
static double METERS_in = 1.0, METERS_out = 1.0;
-/* typedef struct { double u, v; } UV; */
+/* Improved version of pj_do_proj uses the pj_transform function instead
+ * of pj_fwd and pj_inv, to take advantage of datum conversions */
int pj_do_proj(x,y,info_in,info_out)
double *x,*y;
struct pj_info *info_in,*info_out;
{
- int inverse;
- projUV data;
+ int ok;
+ double u,v;
+ double h = 0.0;
+
+ METERS_in = info_in->meters;
+ METERS_out = info_out->meters;
+
+ if (strncmp(info_in->proj,"ll",2) == 0 ) {
+ if (strncmp(info_out->proj,"ll",2) == 0 )
+ ok = pj_transform (info_in->pj, info_out->pj, 1, 0, x, y, &h);
+ else {
+ u = (*x) / RAD_TO_DEG;
+ v = (*y) / RAD_TO_DEG;
+ ok = pj_transform(info_in->pj,info_out->pj, 1, 0, &u, &v, &h);
+ *x = u / METERS_out;
+ *y = v / METERS_out;
+ }
+ }
+ else {
+ if (strncmp(info_out->proj,"ll",2) == 0 ) {
+ u = *x * METERS_in;
+ v = *y * METERS_in;
+ ok = pj_transform(info_in->pj,info_out->pj, 1, 0, &u, &v, &h);
+ *x = u * RAD_TO_DEG;
+ *y = v * RAD_TO_DEG;
+ }
+ else {
+ u = *x * METERS_in;
+ v = *y * METERS_in;
+ ok = pj_transform(info_in->pj,info_out->pj, 1, 0, &u, &v, &h);
+ *x = u / METERS_out;
+ *y = v / METERS_out;
+ }
+ }
+ return ok;
+}
- METERS_in = info_in->meters;
- METERS_out = info_out->meters;
+/* pj_do_transform()
+ * New function to take advantage of height transformations and
+ * transformations of multiple co-ordinates stored in arrays (as afforded by
+ * PROJ.4 function pj_transform().
+ * Not used yet (1/2003). */
- if (strncmp(info_in->proj,"ll",2) == 0 ) {
- if (strncmp(info_out->proj,"ll",2) == 0 )
- /* LL->LL is a no-op */ ;
- else {
- data.u = (*x) / RAD_TO_DEG;
- data.v = (*y) / RAD_TO_DEG;
- data = pj_fwd(data,info_out->pj);
- *x = data.u / METERS_out;
- *y = data.v / METERS_out;
- }
- }
- else {
- if (strncmp(info_out->proj,"ll",2) == 0 ) {
- data.u = *x * METERS_in;
- data.v = *y * METERS_in;
- data = pj_inv(data,info_in->pj);
- *x = data.u * RAD_TO_DEG;
- *y = data.v * RAD_TO_DEG;
- }
- else {
- data.u = *x * METERS_in;
- data.v = *y * METERS_in;
- data = pj_inv(data,info_in->pj);
- data = pj_fwd(data,info_out->pj);
- *x = data.u / METERS_out;
- *y = data.v / METERS_out;
- }
+int pj_do_transform (int count, double *x, double *y, double *h,
+ struct pj_info *info_in, struct pj_info *info_out)
+{
+ int ok;
+ int has_h = 1;
+
+ METERS_in = info_in->meters;
+ METERS_out = info_out->meters;
+
+ if (h == NULL) {
+ int i;
+ h = G_malloc (sizeof *h * count);
+ /* they say memset is only guaranteed for chars ;-( */
+ for (i = 0; i < count; ++i)
+ h[i] = 0.0;
+ has_h = 0;
+ }
+ if (strncmp(info_in->proj,"ll",2) == 0 ) {
+ if (strncmp(info_out->proj,"ll",2) == 0 )
+ ok = pj_transform (info_in->pj, info_out->pj, count, 1, x, y, h);
+ else {
+ DIVIDE_LOOP(x,y,count,RAD_TO_DEG);
+ ok = pj_transform(info_in->pj,info_out->pj, count, 1, x, y, h);
+ DIVIDE_LOOP(x,y,count,METERS_out);
+ }
+ }
+ else {
+ if (strncmp(info_out->proj,"ll",2) == 0 ) {
+ MULTIPLY_LOOP(x,y,count,METERS_in);
+ ok = pj_transform(info_in->pj,info_out->pj, count, 1, x, y, h);
+ MULTIPLY_LOOP(x,y,count,RAD_TO_DEG);
+ }
+ else {
+ MULTIPLY_LOOP(x,y,count,METERS_in);
+ ok = pj_transform(info_in->pj,info_out->pj, count, 1, x, y, h);
+ DIVIDE_LOOP(x,y,count,METERS_out);
+ }
}
- return 1;
+ if (!has_h) free (h);
+
+ return ok;
}
-
-
Only in grass50_exp_2003_01_10/src/libes/proj: do_proj_nad.c
diff -ru grass50_exp_2003_01_10/src/libes/proj/get_proj.c grass/src/libes/proj/get_proj.c
--- grass50_exp_2003_01_10/src/libes/proj/get_proj.c 2002-04-26 04:12:44.000000000 +0100
+++ grass/src/libes/proj/get_proj.c 2003-01-15 13:07:41.521202000 +0000
@@ -7,6 +7,10 @@
#include "projects.h"
#define MAIN
+/* GRASS relative location of datum conversion lookup tables */
+#define GRIDDIR "/etc/nad"
+/* Finder function for datum conversion lookup tables */
+#define FINDERFUNC set_proj_lib
#define PERMANENT "PERMANENT"
#define MAX_PARGS 100
@@ -53,30 +57,60 @@
if (strlen(info->proj) <= 0)
sprintf(info->proj, "ll");
- if (strncmp(proj_in, "Lat", 3) == 0
- || strcmp(info->proj,"ll") == 0) {
- return 1;
- }
nopt1 = 0;
- south = 0;
+ south = 0;
/* if zone is negative, write abs(zone) and define south */
for (i = 0; i < in_proj_keys->nitems; i++) {
/* the name parameter is just for grasses use */
if( strncmp(in_proj_keys->key[i],"name",4) == 0 ) {
continue;
- /* Don't pass through any ellps= settings if we have the
- underlying parameters. */
+ /* Don't need dx dy dz if datum= is specified; PROJ.4 will
+ * get these from pj_datums.c file; otherwise dx dy dz used to
+ * form +towgs84 parameter after end of this loop PK */
+
+ } else if( strncmp(in_proj_keys->key[i], "dx", 2) == 0
+ || strncmp(in_proj_keys->key[i], "dy", 2) == 0
+ || strncmp(in_proj_keys->key[i], "dz", 2) == 0 ) {
+ continue;
+
+
+ /* Don't pass through underlying parameters if we have ellps=
+ * settings; don't pass through ellps= settings if we have
+ * datum= settings. Changed from old logic by PK to use only
+ * top-level data sources if possible (i.e. within PROJ.4.)
+ *
+ * If there are many underlying ellipsoid parameters pass
+ * them all through and let proj sort it out. */
+
+ } else if( strncmp(in_proj_keys->key[i], "a", 1) == 0
+ || strncmp(in_proj_keys->key[i], "b", 1) == 0
+ || strncmp(in_proj_keys->key[i], "es", 2) == 0
+ || strncmp(in_proj_keys->key[i], "f", 1) == 0
+ || strncmp(in_proj_keys->key[i], "rf", 2) == 0 ) {
+ if( G_find_key_value("ellps", in_proj_keys) != NULL
+ || G_find_key_value("datum", in_proj_keys) != NULL)
+ continue;
+ else
+ sprintf(buffa, "%s=%s",
+ in_proj_keys->key[i], in_proj_keys->value[i]);
} else if( strncmp(in_proj_keys->key[i], "ellps", 5) == 0 ) {
- if( G_find_key_value("a", in_proj_keys) != NULL
- && (G_find_key_value("b", in_proj_keys) != NULL
- || G_find_key_value("es", in_proj_keys) != NULL
- || G_find_key_value("rf", in_proj_keys) != NULL))
+ if( G_find_key_value("datum", in_proj_keys) != NULL)
continue;
else
sprintf(buffa, "%s=%s",
in_proj_keys->key[i], in_proj_keys->value[i]);
+
+ /* PROJ.4 uses latlong instead of ll as 'projection name' */
+
+ } else if( strncmp(in_proj_keys->key[i], "proj", 4) == 0 ) {
+ if( strncmp(G_find_key_value("proj", in_proj_keys), "ll", 2) == 0)
+ sprintf(buffa, "%s=latlong", in_proj_keys->key[i]);
+ else
+ sprintf(buffa, "%s=%s",
+ in_proj_keys->key[i], in_proj_keys->value[i]);
+
} else if (strncmp(in_proj_keys->key[i], "south", 5) == 0) {
sprintf(buffa, "south");
south = 1;
@@ -122,6 +156,25 @@
sprintf(opt_in[nopt1-1], buffa);
}
+ /* If datum= is not specified but dx dy dz are, we can use them to
+ * form the towgs84 parameter for datum shift PK */
+ if( G_find_key_value("datum", in_proj_keys) == NULL
+ && G_find_key_value("dx", in_proj_keys) != NULL
+ && G_find_key_value("dy", in_proj_keys) != NULL
+ && G_find_key_value("dz", in_proj_keys) != NULL ) {
+ sprintf(buffa, "towgs84=%s,%s,%s",
+ G_find_key_value("dx", in_proj_keys),
+ G_find_key_value("dy", in_proj_keys),
+ G_find_key_value("dz", in_proj_keys) );
+ nsize = strlen(buffa);
+ if (!(opt_in[nopt1++] = (char *) malloc(nsize + 1)))
+ G_fatal_error("cannot allocate options\n");
+ sprintf(opt_in[nopt1-1], buffa);
+ }
+
+ /* Set finder function for locating datum conversion tables PK */
+ pj_set_finder( FINDERFUNC );
+
if (!(pj = pj_init(nopt1, opt_in))) {
fprintf(stderr, "cannot initialize pj\ncause: ");
fprintf(stderr,"%s\n",pj_strerrno(pj_errno));
@@ -141,7 +194,7 @@
char *s;
int nopt = 0;
int nsize;
- char zonebuff[50];
+ char zonebuff[50], buffa[300];
PJ *pj;
info->zone = 0;
@@ -149,53 +202,71 @@
info->meters = 1.0;
if ((str == NULL) || (str[0] == '\0')) {
+ /* Null string is supplied for parameters, implying latlong
+ * projection; just need to set proj param and call pj_init PK */
sprintf(info->proj, "ll");
- return 1;
- }
- s = str;
- while (s = strtok(s, " \t\n")) {
- if (strncmp(s, "+unfact=", 8) == 0) {
- s = s + 8;
- info->meters = atof(s);
- } else {
- if (strncmp(s, "+", 1) == 0)
- ++s;
- if (nsize = strlen(s)) {
- if (nopt >= MAX_PARGS) {
- fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
- G_fatal_error("Option input overflowed option table");
- }
- if (!(opt_in[nopt] = (char *) malloc(nsize + 1)))
- G_fatal_error("Option input memory failure");
- sprintf(opt_in[nopt++], s);
-
- if (strncmp("proj=", s, 5) == 0) {
- sprintf(info->proj, "%s", s + 5);
- if (strncmp(info->proj, "ll", 2) == 0) {
- sprintf(info->proj, "ll");
- return 1;
+ sprintf(buffa, "proj=latlong");
+ nsize = strlen(buffa);
+ if (!(opt_in[nopt] = (char *) malloc(nsize + 1)))
+ G_fatal_error("Option input memory failure");
+ sprintf(opt_in[nopt++], buffa);
+ } else {
+ /* Parameters have been provided; parse through them but don't
+ * bother with most of the checks in pj_get_kv; assume the
+ * programmer knows what he / she is doing when using this
+ * function rather than reading a PROJ_INFO file PK */
+ s = str;
+ while (s = strtok(s, " \t\n")) {
+ if (strncmp(s, "+unfact=", 8) == 0) {
+ s = s + 8;
+ info->meters = atof(s);
+ } else {
+ if (strncmp(s, "+", 1) == 0)
+ ++s;
+ if (nsize = strlen(s)) {
+ if (nopt >= MAX_PARGS) {
+ fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
+ G_fatal_error("Option input overflowed option table");
}
- }
- if (strncmp("zone=", s, 5) == 0) {
- sprintf(zonebuff, "%s", s + 5);
- sscanf(zonebuff, "%d", &(info->zone));
+
+ if (strncmp("zone=", s, 5) == 0) {
+ sprintf(zonebuff, "%s", s + 5);
+ sscanf(zonebuff, "%d", &(info->zone));
+ }
+
+ if (strncmp("proj=", s, 5) == 0) {
+ sprintf(info->proj, "%s", s + 5);
+ if (strncmp(info->proj, "ll", 2) == 0)
+ sprintf(buffa, "proj=latlong");
+ else
+ sprintf(buffa, s);
+ } else {
+ sprintf(buffa, s);
+ }
+ nsize = strlen(buffa);
+ if (!(opt_in[nopt] = (char *) malloc(nsize + 1)))
+ G_fatal_error("Option input memory failure");
+ sprintf(opt_in[nopt++], buffa);
}
}
+ s = 0;
}
- s = 0;
}
- if (strncmp(info->proj, "ll", 2) != 0) {
- if (!(pj = pj_init(nopt, opt_in))) {
- fprintf(stderr, "cannot initialize pj\ncause: ");
- fprintf(stderr, "%s\n", pj_strerrno(pj_errno));
- return -1;
- }
- info->pj = pj;
+ /* Set finder function for locating datum conversion tables PK */
+ pj_set_finder( FINDERFUNC );
+
+ if (!(pj = pj_init(nopt, opt_in))) {
+ fprintf(stderr, "cannot initialize pj\ncause: ");
+ fprintf(stderr, "%s\n", pj_strerrno(pj_errno));
+ return -1;
}
+ info->pj = pj;
+
return 1;
}
+/* This function is a bit useless now */
int pj_zero_proj(info)
struct pj_info *info;
@@ -207,3 +278,27 @@
return 0;
}
+
+/* set_proj_lib()
+ * 'finder function' for use with PROJ.4 pj_set_finder() function */
+
+const char * set_proj_lib (const char *name)
+{
+ const char *gisbase = G_gisbase();
+ static char *buf = NULL;
+ static int buf_len;
+ size_t len = strlen(gisbase) + sizeof(GRIDDIR)
+ + strlen(name) + 1;
+
+ if( buf_len < len )
+ {
+ if( buf != NULL )
+ G_free( buf );
+ buf_len = len + 20;
+ buf = G_malloc(buf_len);
+ }
+
+ sprintf (buf, "%s%s/%s", gisbase, GRIDDIR, name);
+
+ return buf;
+}
diff -ru grass50_exp_2003_01_10/src/libes/proj/pj_datums.c grass/src/libes/proj/pj_datums.c
--- grass50_exp_2003_01_10/src/libes/proj/pj_datums.c 2002-05-31 04:12:41.000000000 +0100
+++ grass/src/libes/proj/pj_datums.c 2003-01-15 12:34:16.549358000 +0000
@@ -56,7 +56,7 @@
/* id definition ellipse comments */
/* -- ---------- ------- -------- */
"ggrs87", "towgs84=-199.87,74.79,246.62", "grs80", "Greek Geodetic Reference System 1987",
-"nad27", "nadgrids=conus,ntv1 can.dat", "clark66", "North American Datum 1927",
+"nad27", "nadgrids=conus,ntv1_can.dat", "clark66", "North American Datum 1927",
"wgs84", "towgs84=0.0,0.0,0.0", "wgs84", "World Geodetic System 1984",
"wgs72", "towgs84=0.0,0.0,5.0", "wgs72", "World Geodetic System 1972",
"nad83", "towgs84=0.0,0.0,0.0", "grs80", "North American 1983",
@@ -79,5 +79,6 @@
"potsdam", "towgs84=606.0,23.0,413.0", "bessel", "Potsdam Rauenberg 1950 DHDN",
"carthage", "towgs84=-263.0,6.0,431.0", "clark80", "Carthage 1934 Tunisia",
"hermannskogel", "towgs84=653.0,-212.0,449.0", "bessel", "Hermannskogel",
+"ire65", "towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", "modif_airy", "Ireland 1965",
NULL, NULL, NULL, NULL ,
};
diff -ru grass50_exp_2003_01_10/src/libes/proj/pj_ellps.c grass/src/libes/proj/pj_ellps.c
--- grass50_exp_2003_01_10/src/libes/proj/pj_ellps.c 2002-04-26 04:12:53.000000000 +0100
+++ grass/src/libes/proj/pj_ellps.c 2003-01-15 12:35:14.217181000 +0000
@@ -9,7 +9,7 @@
"airy", "a=6377563.396", "rf=299.3249646", "Airy 1830",
"australian", "a=6378160.000", "rf=298.25", "Australian National",
"bessel", "a=6377397.155", "rf=299.1528128", "Bessel 1841",
-"bess nam", "a=6377483.865", "rf=299.1528128", "Bessel 1841, Namibia",
+"bess_nam", "a=6377483.865", "rf=299.1528128", "Bessel 1841, Namibia",
"clark66", "a=6378206.400", "rf=294.9786982", "Clarke 1866",
"clark80", "a=6378249.145", "rf=293.465", "Clarke 1880",
"everest", "a=6377276.345", "rf=300.8017", "Everest 1830, India",
@@ -24,7 +24,7 @@
"indonesian", "a=6378160.000", "rf=298.247", "Indonesian 1974",
"international", "a=6378388.000", "rf=297.0", "International 1924",
"krassovsky", "a=6378245.000", "rf=298.3", "Krassovsky 1940",
-"modif airy", "a=6377340.189", "rf=299.3249646", "Airy 1830 Modified",
+"modif_airy", "a=6377340.189", "rf=299.3249646", "Airy 1830 Modified",
"sam69", "a=6378160.000", "rf=298.25", "South American 1969",
"SAD-69", "a=6378160.000", "rf=298.25", "South American 1969",
"wgs72", "a=6378135.000", "rf=298.26", "World Geodetic System 1972",
@@ -32,7 +32,7 @@
"grs75", "a=6378140.000", "b=6356755.288", "GRS75",
"hayford", "a=6378388.000", "rf=297.000", "Hayford",
"apl4.9", "a=6378137.0", "rf=298.25", "Appl. Phys. 1965",
-"aust sa", "a=6378160.000", "rf=298.25", "Australian Natl & S. Amer. 1969",
+"aust_sa", "a=6378160.000", "rf=298.25", "Australian Natl & S. Amer. 1969",
"andrae", "a=6377104.43", "rf=300.0", "Andrae 1876",
"cpm", "a=6375738.7", "rf=334.29", "Comm. des Poids et Mesures 1799",
"delmbr", "a=6376428.0", "rf=311.5", "Delambre 1810",
@@ -47,7 +47,7 @@
"lerch", "a=6378139.0", "rf=298.257", "Lerch 1976",
"mprts", "a=6397300.0", "rf=191.0", "Maupertius 1738",
"merit", "a=6378137.000", "rf=298.257", "MERIT 1983",
-"new intl", "a=6378157.500", "b=6356772.2", "New International 1967",
+"new_intl", "a=6378157.500", "b=6356772.2", "New International 1967",
"nwl9d", "a=6378145.000", "rf=298.25", "Naval Weapons Lab., 1965",
"plessis", "a=6376523.0", "b=6355863.0", "Plessis 1817",
"seasia", "a=6378155.000", "b=6356773.3205", "South East Asia",
@@ -55,7 +55,7 @@
"walbeck", "a=6376896.000", "b=6355834.8467", "Walbeck",
"wgs60", "a=6378165.000", "rf=298.3", "WGS60",
"wgs66", "a=6378145.000", "rf=298.25", "WGS66",
-"everest m", "a=6377304.063", "rf=300.8017", "Modified Everest",
-"everest p", "a=6377309.613", "rf=300.8017", "Everest (Pakistan)",
+"everest_m", "a=6377304.063", "rf=300.8017", "Modified Everest",
+"everest_p", "a=6377309.613", "rf=300.8017", "Everest (Pakistan)",
0,0,0,0
};
diff -ru grass50_exp_2003_01_10/src/mapdev/v.in.gshhs/main.c grass/src/mapdev/v.in.gshhs/main.c
--- grass50_exp_2003_01_10/src/mapdev/v.in.gshhs/main.c 2002-02-22 04:12:37.000000000 +0000
+++ grass/src/mapdev/v.in.gshhs/main.c 2003-01-15 12:40:03.199002000 +0000
@@ -143,7 +143,6 @@
zone = region.zone;
/* In Info */
-pj_zero_proj(&info_in);
parms_in[0] = '\0';
pj_get_string(&info_in, parms_in);
diff -ru grass50_exp_2003_01_10/src/mapdev/v.proj/local_proto.h grass/src/mapdev/v.proj/local_proto.h
--- grass50_exp_2003_01_10/src/mapdev/v.proj/local_proto.h 2002-04-26 04:12:55.000000000 +0100
+++ grass/src/mapdev/v.proj/local_proto.h 2003-01-15 12:51:15.254350000 +0000
@@ -8,6 +8,3 @@
/* setenv.c */
int select_current_env(void);
int select_target_env(void);
-
-/* datum_shift.c */
-void set_datumshift(char *, char *, char *, char *);
diff -ru grass50_exp_2003_01_10/src/mapdev/v.proj/main.c grass/src/mapdev/v.proj/main.c
--- grass50_exp_2003_01_10/src/mapdev/v.proj/main.c 2002-12-06 04:12:33.000000000 +0000
+++ grass/src/mapdev/v.proj/main.c 2003-01-15 12:53:40.616753000 +0000
@@ -36,7 +36,6 @@
char *gbase;
static char *oform = (char *)0;
char opath[1024];
- char *hold;
char att_file[100], cat_file[100], date[40], mon[4];
FILE *wnd;
struct GModule *module;
@@ -54,14 +53,11 @@
} flag;
char buf[1024];
- /* added for datum conversion */
- char in_datum[64],out_datum[64],in_ellipse[64],out_ellipse[64];
-
G_gisinit (argv[0]);
module = G_define_module();
module->description =
- "Allows projection conversion of vector files (no datum transformation yet).";
+ "Allows projection conversion of vector files.";
/* set up the options and flags for the command line parser */
@@ -229,14 +225,6 @@
G_fatal_error(buffb) ;
}
- /*for datum conversion find input location datum and ellipse */
- *in_datum='\0';
- if((hold=G_database_datum_name()))
- strncpy(in_datum,hold,sizeof(in_datum));
- *in_ellipse='\0';
- if((hold=G_database_ellipse_name()))
- strncpy(in_ellipse,hold,sizeof(in_ellipse));
-
/*** Get projection info for input mapset ***/
in_proj_keys = G_get_projinfo();
if (in_proj_keys == NULL) {
@@ -281,14 +269,6 @@
select_current_env();
- /* for datum conversion find output location datum and ellipse */
- *out_datum='\0';
- if((hold=G_database_datum_name()))
- strncpy(out_datum,hold,sizeof(out_datum));
- *out_ellipse='\0';
- if((hold=G_database_ellipse_name()))
- strncpy(out_ellipse,hold,sizeof(out_ellipse));
-
/****** get the output projection parameters ******/
Out_proj = G_projection();
out_proj_keys = G_get_projinfo();
@@ -341,12 +321,9 @@
oform = "%.10f";
- /* determine which do_proj function to use */
- set_datumshift(in_datum,in_ellipse,out_datum,out_ellipse);
-
/*SE*/
- if(proj_f(&HE,&HS,&info_in,&info_out)<0) {
- fprintf(stderr,"Error in proj_f\n");
+ if(pj_do_proj(&HE,&HS,&info_in,&info_out)<0) {
+ fprintf(stderr,"Error in pj_do_proj\n");
exit(0);
}
E = HE;
@@ -355,8 +332,8 @@
HS = Map.head.S;
/*NE*/
- if(proj_f(&HE,&HN,&info_in,&info_out)<0) {
- fprintf(stderr,"Error in proj_f\n");
+ if(pj_do_proj(&HE,&HN,&info_in,&info_out)<0) {
+ fprintf(stderr,"Error in pj_do_proj\n");
exit(0);
}
N = HN;
@@ -365,8 +342,8 @@
else Out_Map.head.E = HE;
/*SW*/
- if(proj_f(&HW,&HS,&info_in,&info_out)<0) {
- fprintf(stderr,"Error in proj_f\n");
+ if(pj_do_proj(&HW,&HS,&info_in,&info_out)<0) {
+ fprintf(stderr,"Error in pj_do_proj\n");
exit(0);
}
W = HW;
@@ -376,8 +353,8 @@
HW = Map.head.W;
/*NW*/
- if(proj_f(&HW,&HN,&info_in,&info_out)<0) {
- fprintf(stderr,"Error in proj_f\n");
+ if(pj_do_proj(&HW,&HN,&info_in,&info_out)<0) {
+ fprintf(stderr,"Error in pj_do_proj\n");
exit(0);
}
if (HN < N) Out_Map.head.N = N;
@@ -460,8 +437,8 @@
{
X = Points->x[cnt];
Y = Points->y[cnt];
- if(proj_f(&X,&Y,&info_in,&info_out)<0) {
- fprintf(stderr,"Error in proj_f\n");
+ if(pj_do_proj(&X,&Y,&info_in,&info_out)<0) {
+ fprintf(stderr,"Error in pj_do_proj\n");
exit(0);
}
Points->x[cnt] = X;
@@ -495,7 +472,7 @@
strncmp(buffa,"P",1) == 0 )
{
sscanf(buffa,"%s %lf %lf %d", ctype, &X, &Y, &cat);
- proj_f(&X,&Y,&info_in,&info_out);
+ pj_do_proj(&X,&Y,&info_in,&info_out);
/* write line */
sprintf(buffb,"%1s %14.10f %14.10f %7d\n",ctype,X,Y,cat);
fputs(buffb,out1);
diff -ru grass50_exp_2003_01_10/src/misc/m.ll2db/main.c grass/src/misc/m.ll2db/main.c
--- grass50_exp_2003_01_10/src/misc/m.ll2db/main.c 2002-01-22 04:51:17.000000000 +0000
+++ grass/src/misc/m.ll2db/main.c 2003-01-15 12:40:34.506796000 +0000
@@ -149,7 +149,6 @@
}
/* In Info */
-pj_zero_proj(&info_in);
parms_in[0] = '\0';
pj_get_string(&info_in, parms_in);
diff -ru grass50_exp_2003_01_10/src/raster/r.proj/cmd/bordwalk.c grass/src/raster/r.proj/cmd/bordwalk.c
--- grass50_exp_2003_01_10/src/raster/r.proj/cmd/bordwalk.c 2002-06-28 04:14:47.000000000 +0100
+++ grass/src/raster/r.proj/cmd/bordwalk.c 2003-01-15 12:54:31.526154000 +0000
@@ -50,7 +50,7 @@
for (idx=from_hd->west+from_hd->ew_res/2; idx<from_hd->east; idx+=from_hd->ew_res) {
hx = idx;
hy = from_hd->north - from_hd->ns_res/2;
- if (proj_f(&hx, &hy, from_pj, to_pj) < 0)
+ if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
continue;
/* check if we are within the region, but allow for some 'almost inside' points */
/* (should probably be a factor based on input and output resolutions) */
@@ -74,7 +74,7 @@
for (idx=from_hd->north-from_hd->ns_res/2; idx>from_hd->south; idx-=from_hd->ns_res) {
hx = from_hd->east - from_hd->ew_res/2;
hy = idx;
- if (proj_f(&hx, &hy, from_pj, to_pj) < 0)
+ if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
continue;
if (!(hx<to_hd->west-to_hd->ew_res) && !(hx>to_hd->east+to_hd->ew_res) && !(hy<to_hd->south-to_hd->ns_res) && !(hy>to_hd->north+to_hd->ns_res)) {
xmin = !(hx > xmin) ? hx : xmin;
@@ -96,7 +96,7 @@
for (idx=from_hd->east-from_hd->ew_res/2; idx>from_hd->west; idx-=from_hd->ew_res) {
hx = idx;
hy = from_hd->south + from_hd->ns_res/2;
- if (proj_f(&hx, &hy, from_pj, to_pj) < 0)
+ if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
continue;
if (!(hx<to_hd->west-to_hd->ew_res) && !(hx>to_hd->east+to_hd->ew_res) && !(hy<to_hd->south-to_hd->ns_res) && !(hy>to_hd->north+to_hd->ns_res)) {
xmin = !(hx > xmin) ? hx : xmin;
@@ -118,7 +118,7 @@
for (idx=from_hd->south+from_hd->ns_res/2; idx<from_hd->north; idx+=from_hd->ns_res) {
hx = from_hd->west + from_hd->ew_res/2;
hy = idx;
- if (proj_f(&hx, &hy, from_pj, to_pj) < 0)
+ if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
continue;
if (!(hx<to_hd->west-to_hd->ew_res) && !(hx>to_hd->east+to_hd->ew_res) && !(hy<to_hd->south-to_hd->ns_res) && !(hy>to_hd->north+to_hd->ns_res)) {
xmin = !(hx > xmin) ? hx : xmin;
@@ -141,7 +141,7 @@
if (xmin > to_hd->west) {
hx = to_hd->west + to_hd->ew_res/2;
hy = to_hd->south + (to_hd->north - to_hd->south)/2;
- if (!(proj_f(&hx, &hy, to_pj, from_pj) < 0) &&
+ if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
!(hx<from_hd->west) && !(hx>from_hd->east) &&
!(hy<from_hd->south) && !(hy>from_hd->north))
xmin = to_hd->west + to_hd->ew_res/2;
@@ -150,7 +150,7 @@
if (xmax < to_hd->east) {
hx = to_hd->east - to_hd->ew_res/2;
hy = to_hd->south + (to_hd->north - to_hd->south)/2;
- if (!(proj_f(&hx, &hy, to_pj, from_pj) < 0) &&
+ if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
!(hx<from_hd->west) && !(hx>from_hd->east) &&
!(hy<from_hd->south) && !(hy>from_hd->north))
xmax = to_hd->east - to_hd->ew_res/2;
@@ -159,7 +159,7 @@
if (ymin > to_hd->south) {
hx = to_hd->west + (to_hd->east - to_hd->west)/2;
hy = to_hd->south + to_hd->ns_res/2;
- if (!(proj_f(&hx, &hy, to_pj, from_pj) < 0) &&
+ if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
!(hx<from_hd->west) && !(hx>from_hd->east) &&
!(hy<from_hd->south) && !(hy>from_hd->north))
ymin = to_hd->south + to_hd->ns_res/2;
@@ -168,7 +168,7 @@
if (ymax < to_hd->north) {
hx = to_hd->west + (to_hd->east - to_hd->west)/2;
hy = to_hd->north - to_hd->ns_res/2;
- if (!(proj_f(&hx, &hy, to_pj, from_pj) < 0) &&
+ if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
!(hx<from_hd->west) && !(hx>from_hd->east) &&
!(hy<from_hd->south) && !(hy>from_hd->north))
ymax = to_hd->north - to_hd->ns_res/2;
diff -ru grass50_exp_2003_01_10/src/raster/r.proj/cmd/main.c grass/src/raster/r.proj/cmd/main.c
--- grass50_exp_2003_01_10/src/raster/r.proj/cmd/main.c 2002-06-28 04:14:47.000000000 +0100
+++ grass/src/raster/r.proj/cmd/main.c 2003-01-15 12:57:10.297824000 +0000
@@ -62,8 +62,6 @@
#include "projects.h"
#include "r.proj.h"
-extern void set_datumshift(char *, char *, char *, char *);
-
/* modify this table to add new methods */
struct menu menu[] = {
{ p_nearest, "nearest", "nearest neighbor" },
@@ -79,11 +77,7 @@
char *mapname, /* ptr to name of output layer */
*setname, /* ptr to name of input mapset */
*ipolname, /* name of interpolation method */
- errbuf[256], /* buffer for error messages */
- *in_datum, /* data and ellipses for datum */
- *in_ellipse, /* conversion */
- *out_datum,
- *out_ellipse;
+ errbuf[256]; /* buffer for error messages */
int fdi, /* input map file descriptor */
fdo, /* output map file descriptor */
@@ -138,7 +132,7 @@
module = G_define_module();
module->description =
- "re-project a raster map from one location to the current location (no datum transformation yet)";
+ "re-project a raster map from one location to the current location.";
inmap = G_define_option();
@@ -230,12 +224,6 @@
if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
G_fatal_error("Can't get projection key values of output map");
- out_datum = G_database_datum_name();
- out_datum = out_datum ? G_store(out_datum) : "";
-
- out_ellipse = G_database_ellipse_name();
- out_ellipse = out_ellipse ? G_store(out_ellipse) : "";
-
/* Change the location */
G__create_alt_env();
G__setenv("GISDBASE", indbase->answer ? indbase->answer : G_gisdbase());
@@ -275,8 +263,6 @@
/* this call causes r.proj to read the entire map into memeory */
G_get_cellhd(inmap->answer, setname, &incellhd);
- /* this call causes r.proj to use the WIND file settings */
- /*G_get_window(&incellhd);*/
G_set_window(&incellhd);
cell_type = G_raster_map_type(inmap->answer, setname);
@@ -284,15 +270,6 @@
if (G_projection() == PROJECTION_XY)
G_fatal_error("Can't work with xy data");
- in_datum = G_database_datum_name();
- in_datum = in_datum ? G_store(in_datum) : "";
-
- in_ellipse = G_database_ellipse_name();
- in_ellipse = in_ellipse ? G_store(in_ellipse) : "";
-
- /* determine which do_proj function to use */
- set_datumshift(in_datum, in_ellipse, out_datum, out_ellipse);
-
/* Save default borders so we can show them later */
inorth = incellhd.north;
isouth = incellhd.south;
@@ -343,7 +320,7 @@
for(col=0;col<incellhd.cols;col++)
{
xcoord1=G_col_to_easting((double)(col+0.5),&incellhd);
- proj_f(&xcoord1,&ycoord1,&iproj,&oproj);
+ pj_do_proj(&xcoord1,&ycoord1,&iproj,&oproj);
if(xcoord1>outcellhd.east)outcellhd.east=xcoord1;
if(ycoord1>outcellhd.north)outcellhd.north=ycoord1;
if(xcoord1<outcellhd.west)outcellhd.west=xcoord1;
@@ -420,7 +397,7 @@
{
/* project coordinates in output matrix to */
/* coordinates in input matrix */
- if (proj_f(&xcoord1, &ycoord1, &oproj, &iproj) < 0)
+ if (pj_do_proj(&xcoord1, &ycoord1, &oproj, &iproj) < 0)
G_set_null_value(obufptr, 1, cell_type);
else
{
diff -ru grass50_exp_2003_01_10/src/raster/r.sun/main.c grass/src/raster/r.sun/main.c
--- grass50_exp_2003_01_10/src/raster/r.sun/main.c 2002-09-20 04:12:23.000000000 +0100
+++ grass/src/raster/r.sun/main.c 2003-01-15 12:43:41.591971000 +0000
@@ -1302,6 +1302,7 @@
int i, j, l;
/* double energy;*/
double lum, q1;
+ char buf[50];
fprintf(stderr,"\n\n");
@@ -1431,8 +1432,8 @@
if(pj_get_kv(&iproj,in_proj_info,in_unit_info) < 0)
G_fatal_error("Can't get projection key values of current location");
- pj_zero_proj(&oproj);
- sprintf(oproj.proj,"%s","ll");
+ buf[0] = '\0';
+ pj_get_string(&oproj, buf);
longitude = xp;
latitude = yp;
diff -ru grass50_exp_2003_01_10/src/raster/r.sunmask/g_solposition.c grass/src/raster/r.sunmask/g_solposition.c
--- grass50_exp_2003_01_10/src/raster/r.sunmask/g_solposition.c 2002-01-22 04:51:30.000000000 +0000
+++ grass/src/raster/r.sunmask/g_solposition.c 2003-01-15 12:46:08.349517000 +0000
@@ -60,6 +60,7 @@
struct pj_info oproj; /* output map proj parameters */
extern struct Cell_head window;
int inside;
+ char buf[50];
/* we don't like to run G_calc_solar_position in xy locations */
@@ -115,8 +116,8 @@
#endif
/* set output projection to lat/long for solpos*/
- pj_zero_proj(&oproj);
- sprintf(oproj.proj, "%s", "ll");
+ buf[0]='\0';
+ pj_get_string(&oproj, buf);
/* XX do the transform
* outx outy in_info out_info */
diff -ru grass50_exp_2003_01_10/src/sites/s.datum.shift/main.c grass/src/sites/s.datum.shift/main.c
--- grass50_exp_2003_01_10/src/sites/s.datum.shift/main.c 2001-01-23 21:25:57.000000000 +0000
+++ grass/src/sites/s.datum.shift/main.c 2003-01-15 12:47:35.759388000 +0000
@@ -166,8 +166,8 @@
G_fatal_error("Can't get projection key values of output map");
/* set up projection for lat/long reprojection */
- (void) pj_zero_proj(&pjll);
- sprintf(pjll.proj, "ll");
+ errbuf[0]='\0';
+ (void) pj_get_string(&pjll, errbuf);
}
if (molod->answer) {
diff -ru grass50_exp_2003_01_10/src/sites/s.proj/main.c grass/src/sites/s.proj/main.c
--- grass50_exp_2003_01_10/src/sites/s.proj/main.c 2002-05-31 04:13:04.000000000 +0100
+++ grass/src/sites/s.proj/main.c 2003-01-15 12:50:43.359089000 +0000
@@ -43,12 +43,7 @@
main(int argc, char **argv)
{
char *setname, /* ptr to name of input mapset */
- errbuf[256], /* buffer for error messages */
- in_datum[64], /* input datum name */
- in_ellipse[64], /* input ellipse name */
- out_datum[64], /* output datum name */
- out_ellipse[64], /* output datum name */
- *hold;
+ errbuf[256]; /* buffer for error messages */
int permissions, /* mapset permissions */
oldproj, newproj;
@@ -83,7 +78,7 @@
module = G_define_module();
module->description =
"Allows the user to re-project a sites file from one "
- "location to the current location (no datum transformation yet).";
+ "location to the current location.";
newproj = G_projection();
@@ -150,14 +145,6 @@
if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
G_fatal_error("Can't get projection key values of output map");
- /* for datum conversion, get output datum and ellipse */
- *out_datum='\0';
- if((hold=G_database_datum_name()))
- strncpy(out_datum,hold,sizeof(out_datum));
- *out_ellipse='\0';
- if((hold=G_database_ellipse_name()))
- strncpy(out_ellipse,hold,sizeof(out_ellipse));
-
/* Change the location */
G__create_alt_env();
G__setenv("GISDBASE", indbase->answer == NULL
@@ -166,14 +153,6 @@
G__setenv("LOCATION_NAME", inlocation->answer);
permissions = G__mapset_permissions(setname);
- /* for datum conversion, get output datum and ellipse */
- *in_datum='\0';
- if((hold=G_database_datum_name()))
- strncpy(in_datum,hold,sizeof(in_datum));
- *in_ellipse='\0';
- if((hold=G_database_ellipse_name()))
- strncpy(in_ellipse,hold,sizeof(in_ellipse));
-
if (permissions >= 0) {
/* if requested, list the raster files in source location - MN 5/2001*/
@@ -267,15 +246,12 @@
si = G_site_new_struct(rtype,ndim,nstr,ndec);
}
- /* find out which do_proj to use */
- set_datumshift(in_datum,in_ellipse,out_datum,out_ellipse);
-
while (G__site_get (infile, si, oldproj) >= 0) {
xcoord = si->east;
ycoord = si->north;
- if(proj_f(&xcoord, &ycoord, &iproj, &oproj) < 0 ) {
+ if(pj_do_proj(&xcoord, &ycoord, &iproj, &oproj) < 0 ) {
fprintf(stderr,"Error in pj_do_proj\n");
exit(0);
}
More information about the grass-dev
mailing list