[GRASS5] PROJ_INFO and pj_get_kv() and r.in.gdal

Markus_Neteler1 · October 6, 2001, 6:54am

Hi all,

while working with/on r.in.gdal I have realized that three
entries are not written to PROJ_INFO when creating a new
location from a geocoded dataset. These entries are
"name", "a" and "es". Example:

cat PROJ_INFO
name: Transverse Mercator
datum: rome40
dx: -225.000000
dy: -65.000000
dz: 9.000000
proj: tmerc
ellps: international
a: 6378388.0000000000
es: 0.0067226700
f: 297.0000000000
lat_0: 0.0000000000
lon_0: 9.0
k_0: 0.9996000000
x_0: 1500000.0000000000

If "name", "a" and "es" the r.proj, v.proj and s.proj fail because
the function pj_get_kv() fails. This function is in
src/libes/proj/get_proj.c

But...
In my opinion there is no need to change r.in.gdal ( see wkt_to_grass()
function there which defines the proj_info and proj_units struct,
which is later used by G__make_location() ) because
"name", "a" and "es"
are redundant. The needed values are "proj" and "ellps".

Question: Can we simplify pj_get_kv() or do we have to update
wkt_to_grass() in src/raster/r.in.gdal/main.c

If you look into src/raster/r.in.gdal/main.c you see my additions
to optionally re-project GCPs in a SAR image from lat/long to
a user defined target projection. But since pj_get_kv() fails
here due to the missing values, I had to hard-code lat/long as
input "projection" (which might not be true for all SAR data).

Comments would be welcome. I have spend some time to update wkt_to_grass()
but didn't manage (a bit too complex for me...)

Thanks for hints,

Markus

Andreas_Lange1 · October 9, 2001, 11:35am

Hi Markus, hi all,

just a comment:

Markus Neteler wrote:

Hi all,

while working with/on r.in.gdal I have realized that three
entries are not written to PROJ_INFO when creating a new
location from a geocoded dataset. These entries are
"name", "a" and "es". Example:

cat PROJ_INFO
name: Transverse Mercator
datum: rome40
dx: -225.000000
dy: -65.000000
dz: 9.000000
proj: tmerc
ellps: international
a: 6378388.0000000000
es: 0.0067226700
f: 297.0000000000
lat_0: 0.0000000000
lon_0: 9.0
k_0: 0.9996000000
x_0: 1500000.0000000000

If "name", "a" and "es" the r.proj, v.proj and s.proj fail because
the function pj_get_kv() fails. This function is in
src/libes/proj/get_proj.c

But...
In my opinion there is no need to change r.in.gdal ( see wkt_to_grass()
function there which defines the proj_info and proj_units struct,
which is later used by G__make_location() ) because
"name", "a" and "es"
are redundant. The needed values are "proj" and "ellps".

I think that each location must be self-containted in the way that all
parameters should be numerically defined in the PROJ_INFO file. Else you
can get very confusing results if you transfer a location from one
installation to another installation with different etc/ellipse.table
and/or etc/datum.table files. This is possible because e. g. sometimes
there are different interpretations of ellipsoids and/or datums.
IMHO the "name" and "ellps" strings are only for infomational purposes
and all functions should read out the values for calculation. But i
can't remember if the functions really behave that way.
For the ellipsoid you need a and one of es or 1/f or b, but not all. But
it does not harm to have both es and 1/f in the file.

Question: Can we simplify pj_get_kv() or do we have to update
wkt_to_grass() in src/raster/r.in.gdal/main.c

I would ask to fix r.in.gdal instead of messing up the working library
code. Or is the problem in the library code? Or in *.proj modules?

If you look into src/raster/r.in.gdal/main.c you see my additions
to optionally re-project GCPs in a SAR image from lat/long to
a user defined target projection. But since pj_get_kv() fails
here due to the missing values, I had to hard-code lat/long as
input "projection" (which might not be true for all SAR data).

Comments would be welcome. I have spend some time to update wkt_to_grass()
but didn't manage (a bit too complex for me...)

I have no time currently, otherwise i would try to update this myself.

cu,

Andreas

--
Andreas Lange, 65187 Wiesbaden, Germany, Tel. +49 611 807850
url: http://mitglied.tripod.de/AndreasLange
mail: Andreas.Lange_at_Rhein-Main.de - A.C.Lange_at_GMX.net

Markus_Neteler1 · October 9, 2001, 12:06pm

Hi Andreas,

even I prefer to have the wkt_to_grass() function fixed in
r.in.gdal (to also provide "name", "a" and "es" in the proj_info
struct).

I have tried to fix, but didn't manage. Any PROJ experts out
there?

This is one of the few things missing for the stable release.

Ciao

Markus

On Tue, Oct 09, 2001 at 01:35:52PM +0200, Andreas Lange wrote:

Hi Markus, hi all,

just a comment:

Markus Neteler wrote:
>
> Hi all,
>
> while working with/on r.in.gdal I have realized that three
> entries are not written to PROJ_INFO when creating a new
> location from a geocoded dataset. These entries are
> "name", "a" and "es". Example:
>
> cat PROJ_INFO
> name: Transverse Mercator
> datum: rome40
> dx: -225.000000
> dy: -65.000000
> dz: 9.000000
> proj: tmerc
> ellps: international
> a: 6378388.0000000000
> es: 0.0067226700
> f: 297.0000000000
> lat_0: 0.0000000000
> lon_0: 9.0
> k_0: 0.9996000000
> x_0: 1500000.0000000000
>
> If "name", "a" and "es" the r.proj, v.proj and s.proj fail because
> the function pj_get_kv() fails. This function is in
> src/libes/proj/get_proj.c
>
> But...
> In my opinion there is no need to change r.in.gdal ( see wkt_to_grass()
> function there which defines the proj_info and proj_units struct,
> which is later used by G__make_location() ) because
> "name", "a" and "es"
> are redundant. The needed values are "proj" and "ellps".

I think that each location must be self-containted in the way that all
parameters should be numerically defined in the PROJ_INFO file. Else you
can get very confusing results if you transfer a location from one
installation to another installation with different etc/ellipse.table
and/or etc/datum.table files. This is possible because e. g. sometimes
there are different interpretations of ellipsoids and/or datums.
IMHO the "name" and "ellps" strings are only for infomational purposes
and all functions should read out the values for calculation. But i
can't remember if the functions really behave that way.
For the ellipsoid you need a and one of es or 1/f or b, but not all. But
it does not harm to have both es and 1/f in the file.

>
> Question: Can we simplify pj_get_kv() or do we have to update
> wkt_to_grass() in src/raster/r.in.gdal/main.c

I would ask to fix r.in.gdal instead of messing up the working library
code. Or is the problem in the library code? Or in *.proj modules?

>
> If you look into src/raster/r.in.gdal/main.c you see my additions
> to optionally re-project GCPs in a SAR image from lat/long to
> a user defined target projection. But since pj_get_kv() fails
> here due to the missing values, I had to hard-code lat/long as
> input "projection" (which might not be true for all SAR data).
>
> Comments would be welcome. I have spend some time to update wkt_to_grass()
> but didn't manage (a bit too complex for me...)

I have no time currently, otherwise i would try to update this myself.

cu,

Andreas

--
Andreas Lange, 65187 Wiesbaden, Germany, Tel. +49 611 807850
url: http://mitglied.tripod.de/AndreasLange
mail: Andreas.Lange_at_Rhein-Main.de - A.C.Lange_at_GMX.net
_______________________________________________
grass5 mailing list
grass5@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass5

Frank_Warmerdam · October 11, 2001, 2:03am

Markus Neteler wrote:

Hi Andreas,

even I prefer to have the wkt_to_grass() function fixed in
r.in.gdal (to also provide "name", "a" and "es" in the proj_info
struct).

I have tried to fix, but didn't manage. Any PROJ experts out
there?

This is one of the few things missing for the stable release.

Markus,

I appologise for the lengthy delay. I have finally clarified the
definition of es, and where to lookup name and added these along with
a to the wkt_to_grass() function in r.in.gdal. I have done a few tests
when creating new locations and it seems to work fine. Could you try it
out in the situation you were running into the problem?

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

Markus_Neteler1 · October 11, 2001, 9:27am

On Wed, Oct 10, 2001 at 10:03:44PM -0400, Frank Warmerdam wrote:

Markus Neteler wrote:

> Hi Andreas,
>
> even I prefer to have the wkt_to_grass() function fixed in
> r.in.gdal (to also provide "name", "a" and "es" in the proj_info
> struct).
>
> I have tried to fix, but didn't manage. Any PROJ experts out
> there?
>
> This is one of the few things missing for the stable release.

Markus,

I appologise for the lengthy delay. I have finally clarified the
definition of es, and where to lookup name and added these along with
a to the wkt_to_grass() function in r.in.gdal. I have done a few tests
when creating new locations and it seems to work fine. Could you try it
out in the situation you were running into the problem?

Frank,

thanks for the update. I have submitted a change to r.in.gdal to
use pj_get_kv() instead of the hard-coded lat/long.

     /* set input projection*/
     if (pj_get_kv(&iproj, &proj_info, &proj_units) < 0)
         G_fatal_error("Can't get projection key values of input data...");

However, this function still fails:
r.in.gdal ...
[...]
GCPs have the following OpenGIS WKT Coordinate System:
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6326]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG",4326]]
cannot initialize pj
cause: no arguments in initialization list
ERROR: Can't get projection key values of input data (no GCPs present?)

Perhaps it is a pointers task since gcc -Wall reports:

main.c:373: warning: passing arg 2 of 'pj_get_kv' from incompatible pointer
type
main.c:373: warning: passing arg 3 of 'pj_get_kv' from incompatible
pointer type

I have the feeling that &proj_info, &proj_units are not initialized
properly (how to check that?).
Please note: in g.region and r.sunmask I got that working, also
to set the pj_get_kv(&oproj, out_proj_info, out_unit_info)

Pointers are obviously still a secret for me.

Best regards

Markus

Glynn_Clements · October 11, 2001, 10:03am

Markus Neteler wrote:

Perhaps it is a pointers task since gcc -Wall reports:

main.c:373: warning: passing arg 2 of 'pj_get_kv' from incompatible pointer type
main.c:373: warning: passing arg 3 of 'pj_get_kv' from incompatible pointer type

Yep; You're passing "struct Key_Value **" for the second and third
arguments, but it should be "struct Key_Value *".

Try removing the ampersands.

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

Markus_Neteler1 · October 11, 2001, 10:08am

On Thu, Oct 11, 2001 at 11:03:22AM +0100, Glynn Clements wrote:

Markus Neteler wrote:

> Perhaps it is a pointers task since gcc -Wall reports:
>
> main.c:373: warning: passing arg 2 of 'pj_get_kv' from incompatible pointer type
> main.c:373: warning: passing arg 3 of 'pj_get_kv' from incompatible pointer type

Yep; You're passing "struct Key_Value **" for the second and third
arguments, but it should be "struct Key_Value *".

Try removing the ampersands.

Thanks for the hint!
However, I tried already -> core dump...

Best regards

Markus

Markus_Neteler1 · October 11, 2001, 10:16am

On Thu, Oct 11, 2001 at 12:08:54PM +0200, Markus Neteler wrote:

On Thu, Oct 11, 2001 at 11:03:22AM +0100, Glynn Clements wrote:
>
> Markus Neteler wrote:
>
> > Perhaps it is a pointers task since gcc -Wall reports:
> >
> > main.c:373: warning: passing arg 2 of 'pj_get_kv' from incompatible pointer type
> > main.c:373: warning: passing arg 3 of 'pj_get_kv' from incompatible pointer type
>
> Yep; You're passing "struct Key_Value **" for the second and third
> arguments, but it should be "struct Key_Value *".
>
> Try removing the ampersands.

Thanks for the hint!
However, I tried already -> core dump...

Here a gdb output:

Program received signal SIGSEGV, Segmentation fault.
0x08053786 in G_find_key_value (key=0x8094c85 "meters", kv=0x0)
at key_value1.c:106
106 for (n = 0; n < kv->nitems; n++)

key_value1.c is in /src/libes/gis I think.

Shall I test more with gdb?

Markus

Frank_Warmerdam · October 11, 2001, 3:02pm

Markus Neteler wrote:

Here a gdb output:

Program received signal SIGSEGV, Segmentation fault.
0x08053786 in G_find_key_value (key=0x8094c85 "meters", kv=0x0)
at key_value1.c:106
106 for (n = 0; n < kv->nitems; n++)

key_value1.c is in /src/libes/gis I think.

Shall I test more with gdb?

Markus,

The code is fundamentally flawed (proj_info is usually uninitialized).

I will restructure this morning.

Later,

--
---------------------------------------+--------------------------------------
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

Frank_Warmerdam · October 11, 2001, 4:05pm

Folks,

o I have restructured r.in.gdal so that the GCP reprojection should work
properly.

  o I have modified pj_get_kv() so that it works even if name is not in the
    key/value list, or if it is out of order. The old code seemed to depend
    on name: being the first key in the list. r.in.gdal now generates the
    name parameter in the PROJ_INFO file, but doesn't generally put it first.

Let me know if any problems remain.

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

Markus_Neteler1 · October 12, 2001, 3:50pm

Frank,

thank you a lot for finishing the GCPs support in r.in.gdal.

On Thu, Oct 11, 2001 at 12:05:04PM -0400, Frank Warmerdam wrote:

Folks,

o I have restructured r.in.gdal so that the GCP reprojection should work
properly.

yes, it does. Today was the final run to import the SAR, rectification
also went properly.

  o I have modified pj_get_kv() so that it works even if name is not in the
    key/value list, or if it is out of order. The old code seemed to depend
    on name: being the first key in the list. r.in.gdal now generates the
    name parameter in the PROJ_INFO file, but doesn't generally put it first.

Thanks for this fix as well.

My apologize to take so much of your time here. However, the module
is again more user friendly now.

Best regards

Markus