Cannot create a new GRASS Location using a custom proj4 string

[voncasec]

I have tried this using both grass6.4.3svn release branch (revision 54578) and grass7.0svn trunk (revision 54587).

Using the wxGUI: Location Wizard -> [input location name] -> Specify projection and datum terms using custom PROJ.4 parameters -> [input proj.4 string (eg. +proj=utm +zone=10 +datum=nad83)] -> NEXT

Error: No output format specified, define one of flags -p,-g,-j or -w

[annakrat]

This patch seems to solve this problem, however I don't really understand what's going on here. There are also other places when g.proj is called without specifying flags.

--- location_wizard/wizard.py   (revision 54601)
+++ location_wizard/wizard.py   (working copy)
@@ -1535,7 +1535,8 @@
             ret, out, err = RunCommand('g.proj',
                                        read = True, getErrorMsg = True,
                                        proj4 = self.customstring, 
-                                       datum_trans = '-1')
+                                       datum_trans = '-1',
+                                       flags = 'jf')

[hamish]

Hi, testing with the latest relbr64 and devbr6,

in the wxPy loc'n wizard, defining by proj4 terms, using the string:

+proj=utm +zone=10 +datum=nad83

it asks me to select from transform options 0-6, then when I press Ok to that, I get a slightly different error printed in the terminal:

Traceback (most recent call last):
  File "/home/hbowman/dev/grass/svn/relbr_6_4/dist.x86_64-unknown-linux-gnu/
  etc/wxpython/location_wizard/wizard.py", line 1517, in OnPageChanging

    if len(dtrans) == 0:
TypeError: object of type 'int' has no len()

the error that it's checking for is 'Datum transform is required.', -- but I don't think it's true, grass is happy enough to let you leave the datum transform undefined in the PROJ_INFO file. (then it gives you a message about it, "default transform will be used") ??

If I change the +datum= to wgs84 I get the same "Error: No output format specified, define one of flags -p,-g,-j or -w" error box pops up. (that message is from g.proj)

from the grass6 command line, this shows the 0-6 transform list:

g.proj proj4='+proj=utm +zone=10 +datum=nad83' datumtrans=-1

but this shows the 'please use a flag' error message:

g.proj proj4='+proj=utm +zone=10 +datum=wgs84' datumtrans=-1

Adding the g.proj '-t' flag to "Force override of datum transformation information in input co-ordinate system" makes it print the transform opts correctly (then exit) for both. I assume the wgs84 continues on as it only has the 1 case to choose from, so it gets back to main() and then wonders what to do next since no flags were given.

So a quick fix is to add -t, but why is that needed?

The problem seems to have come in with r53453; adding '-jf' flags is not the correct thing, that just makes it fail differently.

Hamish

[hamish]

some debug:

# ok datum name, mult options:
G:g.proj > g.proj proj4='+proj=utm +zone=10 +datum=nad83' datumtrans=-1
..force=0
D0/0: set_datumtrans(): GPJ__get_datum_params() status=1
---
1
Used in whole nad83 region
towgs84=0.000,0.000,0.000
Default 3-Parameter Transformation (May not be optimum for older datums; use this only if no more appropriate options are available.)
---
2
Used in Florida
[...]

(ok)


# bogus datum name:
G:g.proj > g.proj proj4='+proj=utm +zone=10 +datum=nad839999' datumtrans=-1
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
..force=0
D0/0: set_datumtrans(): GPJ__get_datum_params() status=-1
D0/0: set_datumtrans(): Datum name either invalid or not supplied.
ERROR: No output format specified, define one of flags -p, -g, -j, or -w

(doesn't quit when it should)



# ok datum name, single option:
G:g.proj > g.proj proj4='+proj=utm +zone=10 +datum=wgs84' datumtrans=-1
..force=0
D0/0: set_datumtrans(): GPJ__get_datum_params() status=1
ERROR: No output format specified, define one of flags -p, -g, -j, or -w



# same as above w/ -t
G:g.proj > g.proj proj4='+proj=utm +zone=10 +datum=wgs84' datumtrans=-1 -t
..force=1
D0/0: set_datumtrans(): GPJ__get_datum_params() status=1
---
1
Used in whole wgs84 region
towgs84=0.000,0.000,0.000
Default 3-Parameter Transformation (May not be optimum for older datums; use this only if no more appropriate options are available.)

[hamish]

suggestion: for 6.4.3 throw the '-t' flag in there as a quick hack-around to get it working for release, then continue to work on the real problem in dev branches.

Hamish

[hamish]

Added the -t flag to g.proj in devbr6 for testing, now for the utm 10 example above I still get the python error "object of type 'int' has no len()" from comment:2, but at least the proj4 terms are populated and the location creates itself ok. I notice that after selecting the default 3-term datum transform that +towgs84=0,0,0,0,0,0,0. No problem since it's the same as 0,0,0, but perhaps a sign of trouble elsewhere. -- it is, trying again with the Florida transform grid ends up with the same result, +towgs84=0,0,0,0,0,0,0 and no grid file selected in the PROJ_INFO file.

different than in comment:2 (which was run on linux), on wingrass 6.5svn I get the "define one of flags -p,-g,-j or -w" error box for nad83 too.

Hamish

[hamish]

besides the one on line 1500, maybe also needed in wizard.py line 901, 1346?

tried manually adding to latest wingrass 6.4svn build, selected the grid file by epsg code 27200, but it jumped back to the 7-term default proj4 terms anyway. argh.

specifying from the search list worked (yay), specifying by proj4 terms fails with:

...
    if len(dtrans) == 0:
TypeError: object of type 'int' has no len()

as in comment:2.

Hamish

[hamish]

Replying to hamish:

besides the one on line 1500, maybe also needed in wizard.py line 901, 1346?

committed to relbr6 in r56143.

specifying by proj4 terms fails with:

...
    if len(dtrans) == 0:
TypeError: object of type 'int' has no len()

as in comment:2.

no idea about the above Traceback error, but we're getting closer to having something working for 6.4.3.

Hamish

[hamish]

current situation of loc wiz + datum transform options in devbr6 nightly wingrass build: (which uses proj4 4.8.0)

test 1: datum transform selection

proj: nzmg, datum: nzgd49, (epsg: 27200) datum transform: "2" NTv2 grid file

• setup by terms: works!
• setup by epsg code: grid file selected but +towgs84= 7-term in Summary page & loc'n.
• setup by proj terms: "+proj=nzmg +datum=nzgd49"; grid file selected but no sign of it in the Summary page and loc'n created with the +towgs84 7-term transform.

test 2: no datum

• setup by terms: proj -> ll, ellipsoid only -> sphere -> error. (see trac #1967)

these two things are the last majors blocking release AFAICT.

Hamish

[hamish]

traceback when specifying by proj4 terms hopefully fixed in devbr6 by r56160.

[hamish]

(r56160 fixed by r56161, the "test 2: no datum" case is hopefully now fixed too in devbr6 and trunk)

datum transform selection by epsg code and proj4 terms still broken.

Hamish

[hamish]

setup by custom proj terms hopefully fixed & working ok in devbr6 as of r56176 (please test with your local!), similar fixes for epsg remain todo.

an existing bug in the custom proj terms (maybe elsewhere too) is if you select datum transform opt 0 from the list (i.e. don't specify transform opts) it picks the grass default one anyway.

Hamish

[hamish]

Replying to hamish:

setup by custom proj terms hopefully fixed & working ok in devbr6 as of r56176 (please test with your local!), similar fixes for epsg remain todo.

epsg, custom proj terms still not right in wingrass+proj 4.8.0, but looking ok on linux+proj 4.7.0.

It seems that the g.proj -t flag is needed with the 'g.proj -jw' calls to get the datum transform on the summary page on wingrass+4.8.0, and the -t flag is also needed in grass.create_location() (in core.py, at least for epsg-mode) for the g.proj call there.

On wingrass 6.5 I tried 'g.proj -c epsg=27200 datumtrans=2 loc=test_without_t' right on the command line, and without using the -t flag meant that the datum transform was not added to the resulting PROJ_INFO file, instead the default proj4-heuristics option was used. With the -t flag it respected the datumtrans= option.

an existing bug in the custom proj terms (maybe elsewhere too) is if you select datum transform opt 0 from the list (i.e. don't specify transform opts) it picks the grass default one anyway.

this is a bug in g.proj ignoring datumtrans=0, or rather falling back to the grass datum table default (this time, proj4 heuristics choose differently). As I'm not entirely sure if that should be classed as a bug or intended feature, I'll file it as another ticket. (seems ok in wingrass 6.5 + proj 4.8.0 but not linux 6.5 + proj 4.7.0)

Hamish

[hamish]

Hi,

'g.proj -t' work-around in the wxGUI location wizard now applied in all branches, which seems to make things better when used with proj 4.8.0.

please everyone test all ways of creating all types of locations, both with proj 4.7.0 and 4.8.0. Don't trust the summary screen, you need to confirm the result was as expected with 'g.proj -p' from within the new session.

still haven't done anything about datumtrans=0, my feeling is that it should remove all datum transform terms to match its description "Do not apply any datum transforms", so leave it undefined even if proj 4.8.0 has decided to choose some default terms itself (i.e. strip away the trying-to-be-helpfulness).

thanks, Hamish

[hamish]

Helmut wrote:

some teste here ​http://lists.osgeo.org/pipermail/grass-user/2013-May/068240.html ​http://lists.osgeo.org/pipermail/grass-user/2013-May/068241.html ​http://lists.osgeo.org/pipermail/grass-user/2013-May/068242.html

...

EPSG 3416 ETRS89 / Austria Lambert

can you comment on if the results are as you expect? the PROJ_INFO file (g.proj -p) is the main result to be concerned with.

Moritz added:

​http://lists.osgeo.org/pipermail/grass-user/2013-May/068239.html

Helmut noticed the same:

it seems "0 Do not apply any datum transformation" isn't applied?

the PROJ.4 function called by 'g.proj -c' sees there is missing datum transform opts and decides that it should add one. Which one it decides to add depends on your version of PROJ.4.

filed as a new bug for 6.4.4 (#1995).

The old text based g.setproj location creation tool doesn't have this problem since it creates the PROJ_INFO from terms directly. The wxGUI location wizard during location creation goes something like grass terms -> +proj terms (summary page) -> WKT -> grass terms. ... much more opportunity for something to happen, but fortunately the coefficients are pretty well defined. Also it is different if you define by terms (parsed via g.proj), or set +proj terms directly (not parsed until the last step, so no defaults added). e.g. if you define the location by +proj4 terms, and enter "+proj=nzmg +datum=nzgd49" and select "0: do not apply transform" the summary page (parsed by g.proj) shows added default transform terms, but the actual PROJ_INFO in the location (correctly) doesn't have them.

without abandoning the 'g.proj -c' method, I'm not sure of the best way to solve the problem of avoiding PROJ.4 being helpful when you want it to be dumb. maybe the best thing to do is to tell people to look at 'g.proj -p' just to verify it's what you wanted in a final pop up after the location is created (so we can just display the PROJ_INFO file as-is), but before it asks about setting default bounds and creating a new mapset.

(actually it may be OGR channeling PROJ 4.8.0 via GPJ_wkt_to_grass() and GPJ_osr_to_grass(), and the 4.8.0 thing mostly when starting with EPSG codes w/o specified transform terms)

Select WKT or PRJ file

...

no asking which transformation should be used.

ok, filed that as a new bug for 6.4.4 (#1994).

"Used in whole hermannskogel region: towwgs84=653.000,-212.000,449.000" is preselected, but not the best e.g. if location used for Austria.

that gets back to the whole mess of trying to decide which is the best to choose. proj 4.8.0 likes 7-terms (may be more "correct" in smaller areas) over 3-terms (may be less accuracy at focus but covers a much wider area). And sometimes the 7-term for a place is great and the 3-term was poorly derived. It's hard to say. Generally if using the old-ways GRASS's internal will suggest the 3-term as the default. See text comments at the start of $GISBASE/etc/datum.table. There is no real right or wrong answer, the user will have to do their homework and make an informed decision for themselves.

(it's all a bit complicated, the above is my best guess and probably not 100% accurate)

The original +proj=utm +zone=10 +datum=nad83 problem in this ticket is now fixed in all versions, and besides for the two new tickets most things seem to be working as expected, so closing the ticket.

Perhaps it would still be good to pop up a window as soon as the location is created displaying what is in the final PROJ_INFO file, just to be sure. it wouln't have to say that we somewhat doubt our code, just be informative saying this is what it ended up with, "ok, continue".

Hamish