[GRASS-dev] GRASS GIS + PROJ 6 + GDAL 2.5

GRASS GRASS developer
,
1

Hi all,

GRASS needs some adjustments in order to be compatible with PROJ 6 + GDAL 2.5

The plain text file “share/proj/epsg” no longer exists. This file is currently required by the GUI location wizard to retrieve a list of known EPSG codes. Now there is a sqlite db “proj.db”, and a new PROJ function proj_get_codes_from_database(). I suggest to enhance g.proj with a new option “list_codes_auth” which will print a list of codes for the given autority name (EPSG, IGNF, etc.). This option can be made backwards compatible to read “share/proj/epsg” with PROJ versions up to 5. The location wizard can then use this new output of g.proj to construct a list of known EPSG codes.

We need to take care of axis order if a CRS is created from EPSG because the axis order can now be northing, easting instead of traditional easting, northing. Maybe it is enough to enforce +axis=enu, or to use GDAL’s SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) when asking GDAL to construct a CRS definition from an EPSG code.

Markus M

2

Hi,

po 25. 2. 2019 v 15:10 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

GRASS needs some adjustments in order to be compatible with PROJ 6 + GDAL 2.5

The plain text file "share/proj/epsg" no longer exists. This file is currently required by the GUI location wizard to retrieve a list of known EPSG codes. Now there is a sqlite db "proj.db", and a new PROJ function proj_get_codes_from_database(). I suggest to enhance g.proj with a new option "list_codes_auth" which will print a list of codes for the given autority name (EPSG, IGNF, etc.). This option can be made backwards compatible to read "share/proj/epsg" with PROJ versions up to 5. The location wizard can then use this new output of g.proj to construct a list of known EPSG codes.

sounds like a good plan :slight_smile: Worth to create new enhancement ticket? Martin

--
Martin Landa
http://geo.fsv.cvut.cz/gwiki/Landa
http://gismentors.cz/mentors/landa

3

On lundi 25 février 2019 15:10:10 CET Markus Metz wrote:

Hi all,

GRASS needs some adjustments in order to be compatible with PROJ 6 + GDAL
2.5

The plain text file "share/proj/epsg" no longer exists. This file is
currently required by the GUI location wizard to retrieve a list of known
EPSG codes. Now there is a sqlite db "proj.db", and a new PROJ function
proj_get_codes_from_database(). I suggest to enhance g.proj with a new
option "list_codes_auth" which will print a list of codes for the given
autority name (EPSG, IGNF, etc.). This option can be made backwards
compatible to read "share/proj/epsg" with PROJ versions up to 5. The
location wizard can then use this new output of g.proj to construct a list
of known EPSG codes.

If you want to do a nice GUI, proj_get_crs_info_list_from_database() is an
enhanced version of proj_get_codes_from_database(), that will retrieve
more information besides the code: name, area of use, type of CRS (geographic,
projected, compound, etc...)
in a very fast way (you can do the same with proj_get_codes_from_database(),
and constructing a CRS for each code, but this is slower)

We need to take care of axis order if a CRS is created from EPSG because
the axis order can now be northing, easting instead of traditional easting,
northing. Maybe it is enough to enforce +axis=enu,

You can only use this if you use a PROJ string when invokating
proj_create_crs_to_crs(), and this is the default axis order for PROJ string anyway.

or to use GDAL's
SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)

Yes, if you use GDAL OGRSpatialReference this is the way to go.

If you don't use GDAL SRS and just PROJ API, you might also just reuse
https://github.com/pramsey/postgis/blob/svn-trunk-proj/liblwgeom/lwgeom_transform.c#L298

which is called for the source and target CRS of the PJ object created
with proj_create_crs_to_crs()
https://github.com/pramsey/postgis/blob/svn-trunk-proj/liblwgeom/lwgeom_transform.c#L357

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

4

On Mon, Feb 25, 2019 at 9:42 PM Even Rouault <even.rouault@spatialys.com> wrote:

On lundi 25 février 2019 15:10:10 CET Markus Metz wrote:

Hi all,

GRASS needs some adjustments in order to be compatible with PROJ 6 + GDAL
2.5

The plain text file “share/proj/epsg” no longer exists. This file is
currently required by the GUI location wizard to retrieve a list of known
EPSG codes. Now there is a sqlite db “proj.db”, and a new PROJ function
proj_get_codes_from_database(). I suggest to enhance g.proj with a new
option “list_codes_auth” which will print a list of codes for the given
autority name (EPSG, IGNF, etc.). This option can be made backwards
compatible to read “share/proj/epsg” with PROJ versions up to 5. The
location wizard can then use this new output of g.proj to construct a list
of known EPSG codes.

If you want to do a nice GUI, proj_get_crs_info_list_from_database() is an
enhanced version of proj_get_codes_from_database(), that will retrieve
more information besides the code: name, area of use, type of CRS (geographic,
projected, compound, etc…)
in a very fast way (you can do the same with proj_get_codes_from_database(),
and constructing a CRS for each code, but this is slower)

Thanks for the hint!

I have added a new option list_codes to g.proj that will print codes, names, and proj definitions for the given authority in trunk r74174. This new option works with PROJ 4, 5, and 6.

The GUI is now (r74175) no longer reading the file share/proj/epsg which no longer exists in PROJ 6, instead it uses g.proj list_codes=EPSG to get the list of EPSG codes. This is important e.g. for the location wizard to provide a list of known EPSG codes.

The GUI in GRASS 7.6 is still trying to read share/proj/epsg and thus not compatible with PROJ 6.

There are other authorities besides EPSG, these are supported with g.proj list_codes= if compiled against PROJ 6. Maybe we should support them too.

I will look at the issues mentioned below in the next few days, should not be too much work I hope.

Markus M

We need to take care of axis order if a CRS is created from EPSG because
the axis order can now be northing, easting instead of traditional easting,
northing. Maybe it is enough to enforce +axis=enu,

You can only use this if you use a PROJ string when invokating
proj_create_crs_to_crs(), and this is the default axis order for PROJ string anyway.

or to use GDAL’s
SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)

Yes, if you use GDAL OGRSpatialReference this is the way to go.

If you don’t use GDAL SRS and just PROJ API, you might also just reuse
https://github.com/pramsey/postgis/blob/svn-trunk-proj/liblwgeom/lwgeom_transform.c#L298

which is called for the source and target CRS of the PJ object created
with proj_create_crs_to_crs()
https://github.com/pramsey/postgis/blob/svn-trunk-proj/liblwgeom/lwgeom_transform.c#L357

Even

–
Spatialys - Geospatial professional services
http://www.spatialys.com

5

On Thu, Mar 7, 2019 at 5:17 PM Markus Metz <markus.metz.giswork@gmail.com> wrote:

On Mon, Feb 25, 2019 at 9:42 PM Even Rouault <even.rouault@spatialys.com> wrote:

On lundi 25 février 2019 15:10:10 CET Markus Metz wrote:

Hi all,

GRASS needs some adjustments in order to be compatible with PROJ 6 + GDAL
2.5

The plain text file “share/proj/epsg” no longer exists. This file is
currently required by the GUI location wizard to retrieve a list of known
EPSG codes. Now there is a sqlite db “proj.db”, and a new PROJ function
proj_get_codes_from_database(). I suggest to enhance g.proj with a new
option “list_codes_auth” which will print a list of codes for the given
autority name (EPSG, IGNF, etc.). This option can be made backwards
compatible to read “share/proj/epsg” with PROJ versions up to 5. The
location wizard can then use this new output of g.proj to construct a list
of known EPSG codes.

If you want to do a nice GUI, proj_get_crs_info_list_from_database() is an
enhanced version of proj_get_codes_from_database(), that will retrieve
more information besides the code: name, area of use, type of CRS (geographic,
projected, compound, etc…)
in a very fast way (you can do the same with proj_get_codes_from_database(),
and constructing a CRS for each code, but this is slower)

Thanks for the hint!

I have added a new option list_codes to g.proj that will print codes, names, and proj definitions for the given authority in trunk r74174. This new option works with PROJ 4, 5, and 6.

Unfortunately it does not work yet with PROJ 4 because new fns have been added to the old proj_api.h…

Markus M

The GUI is now (r74175) no longer reading the file share/proj/epsg which no longer exists in PROJ 6, instead it uses g.proj list_codes=EPSG to get the list of EPSG codes. This is important e.g. for the location wizard to provide a list of known EPSG codes.

The GUI in GRASS 7.6 is still trying to read share/proj/epsg and thus not compatible with PROJ 6.

There are other authorities besides EPSG, these are supported with g.proj list_codes= if compiled against PROJ 6. Maybe we should support them too.

I will look at the issues mentioned below in the next few days, should not be too much work I hope.

Markus M

We need to take care of axis order if a CRS is created from EPSG because
the axis order can now be northing, easting instead of traditional easting,
northing. Maybe it is enough to enforce +axis=enu,

You can only use this if you use a PROJ string when invokating
proj_create_crs_to_crs(), and this is the default axis order for PROJ string anyway.

or to use GDAL’s
SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)

Yes, if you use GDAL OGRSpatialReference this is the way to go.

If you don’t use GDAL SRS and just PROJ API, you might also just reuse
https://github.com/pramsey/postgis/blob/svn-trunk-proj/liblwgeom/lwgeom_transform.c#L298

which is called for the source and target CRS of the PJ object created
with proj_create_crs_to_crs()
https://github.com/pramsey/postgis/blob/svn-trunk-proj/liblwgeom/lwgeom_transform.c#L357

Even

–
Spatialys - Geospatial professional services
http://www.spatialys.com

6

Markus Metz:

The GUI is now (r74175) no longer reading the file share/proj/epsg which
no longer exists in PROJ 6, instead it uses g.proj list_codes=EPSG to get
the list of EPSG codes. This is important e.g. for the location wizard to
provide a list of known EPSG codes.

Obviously, this works also in the TUI (aka the command line) and looks
great!

→ g.proj list_codes=EPSG |grep 2100
2100|GGRS87 / Greek Grid|+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +datum=GGRS87 +units=m +no_defs
3035|ETRS89 / LAEA Europe|+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
5633|PTRA08 / LAEA Europe|+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
5635|REGCAN95 / LAEA Europe|+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
5636|TUREF / LAEA Europe|+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
5638|ISN2004 / LAEA Europe|+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
32100|NAD83 / Montana|+proj=lcc +lat_1=49 +lat_2=45 +lat_0=44.25 +lon_0=-109.5 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs

Thanks

7

Since rgdal::make_EPSG() is facing the same problems of listing tabulated
EPSG fields as g.proj -l, I was very happy to see Markus' code in
g.proj/main.c mentioned in this thread, and have used this approach in
https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=markup&root=rgdal

However, there are plenty of messages such as: "proj_as_proj_string:
Unsupported conversion method: Lambert Conic Conformal (West Orientated)". I
haven't installed GRASS trunk with PROJ6, so I can't see whether g.proj -l
also sees the same messages. If it does, maybe we could ask on the proj list
how they might be captured for summary reporting. I think they are coming
from line 5758 in src/iso19111/coordinateoperation.cpp or maybe line 906 in
same file. Maybe PROJ now has an error handler that

Another question concerns the issue of whether one needs to free objects
created, in particular proj_crs_info and pj. Not so important for g.proj,
which exists when done, but important for rgdal whose functions don't exit.

Anyway, very helpful to see that Markus is looking at the same issues as we
are!

Roger

-----
Roger Bivand
NHH Norwegian School of Economics, Bergen, Norway
--
Sent from: http://osgeo-org.1560.x6.nabble.com/Grass-Dev-f3991897.html

8

On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no> wrote:

Since rgdal::make_EPSG() is facing the same problems of listing tabulated
EPSG fields as g.proj -l, I was very happy to see Markus’ code in
g.proj/main.c mentioned in this thread, and have used this approach in
https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=markup&root=rgdal

However, there are plenty of messages such as: “proj_as_proj_string:
Unsupported conversion method: Lambert Conic Conformal (West Orientated)”. I
haven’t installed GRASS trunk with PROJ6, so I can’t see whether g.proj -l

also sees the same messages.

I get the same messages with g.proj, but no error code from PROJ error handlers.

If it does, maybe we could ask on the proj list

yes, we should ask. For now g.proj is ignoring those codes for which no proj string can be retrieved.

Markus M

how they might be captured for summary reporting. I think they are coming
from line 5758 in src/iso19111/coordinateoperation.cpp or maybe line 906 in
same file. Maybe PROJ now has an error handler that

Another question concerns the issue of whether one needs to free objects
created, in particular proj_crs_info and pj. Not so important for g.proj,
which exists when done, but important for rgdal whose functions don’t exit.

Anyway, very helpful to see that Markus is looking at the same issues as we
are!

Roger

Roger Bivand
NHH Norwegian School of Economics, Bergen, Norway

Sent from: http://osgeo-org.1560.x6.nabble.com/Grass-Dev-f3991897.html

grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

9

On vendredi 8 mars 2019 12:27:37 CET Markus Metz wrote:

On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no> wrote:
> Since rgdal::make_EPSG() is facing the same problems of listing tabulated
> EPSG fields as g.proj -l, I was very happy to see Markus' code in
> g.proj/main.c mentioned in this thread, and have used this approach in

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=mar
kup&root=rgdal
> However, there are plenty of messages such as: "proj_as_proj_string:
> Unsupported conversion method: Lambert Conic Conformal (West

Orientated)". I

> haven't installed GRASS trunk with PROJ6, so I can't see whether g.proj -l
> also sees the same messages.

I get the same messages with g.proj, but no error code from PROJ error
handlers.

> If it does, maybe we could ask on the proj list

yes, we should ask. For now g.proj is ignoring those codes for which no
proj string can be retrieved.

I'm not sure what the question is exactly :slight_smile: Is that proj_context_errno()
doesn't return an error in that case ?
That's related to the error being emitted from the new code emitting C++
exceptions, which are caught by the C wrapper. It doesn't set errno style
errors currently. I'm not a big fan of the errno approach, and find the proper
error string currently emitted to be more expressive, but that can be
discussed if some harmonization should be done.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

10

On Fri, 8 Mar 2019, Even Rouault wrote:

On vendredi 8 mars 2019 12:27:37 CET Markus Metz wrote:

On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no> wrote:

Since rgdal::make_EPSG() is facing the same problems of listing tabulated
EPSG fields as g.proj -l, I was very happy to see Markus' code in
g.proj/main.c mentioned in this thread, and have used this approach in

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=mar
kup&root=rgdal

However, there are plenty of messages such as: "proj_as_proj_string:
Unsupported conversion method: Lambert Conic Conformal (West

Orientated)". I

haven't installed GRASS trunk with PROJ6, so I can't see whether g.proj -l
also sees the same messages.

I get the same messages with g.proj, but no error code from PROJ error
handlers.

If it does, maybe we could ask on the proj list

yes, we should ask. For now g.proj is ignoring those codes for which no
proj string can be retrieved.

I'm not sure what the question is exactly :slight_smile: Is that proj_context_errno()
doesn't return an error in that case ?

Yes, for example:

...
     for (i = 0; i < crs_cnt; i++) {
         const char *proj_definition;

         pj = proj_create_from_database(ctx, proj_crs_info[i]->auth_name,
             proj_crs_info[i]->code, PJ_CATEGORY_CRS, 0, NULL);
Rprintf("proj_create_from_database %s: %d\n", proj_crs_info[i]->code, proj_context_errno(ctx));
         proj_definition = proj_as_proj_string(ctx, pj, PJ_PROJ_5, NULL);
Rprintf("proj_as_proj_string %s: %d\n", proj_crs_info[i]->code, proj_context_errno(ctx));
...

showing in one case:

proj_create_from_database 7082: 0
proj_as_proj_string: Unsupported conversion method: Polar Stereographic (variant C)
proj_as_proj_string 7082: 0

so the message is being delivered to the console (stderr?), no errno (or other notification seems to be set. Is there a way of checking first before calling proj_as_proj_string() to avoid the message?

I did find the *_destroy() functions instead of legacy pj_free() etc., thanks for prompting me to re-read proj.h.

I'm assuming there is no error handler to hook onto, to divert errors to R's error handler. In this case, it doesn't matter for replacing legacy functionality.

Roger

That's related to the error being emitted from the new code emitting C++
exceptions, which are caught by the C wrapper. It doesn't set errno style
errors currently. I'm not a big fan of the errno approach, and find the proper
error string currently emitted to be more expressive, but that can be
discussed if some harmonization should be done.

Even

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

11

On Fri, Mar 8, 2019 at 12:41 PM Even Rouault <even.rouault@spatialys.com> wrote:

On vendredi 8 mars 2019 12:27:37 CET Markus Metz wrote:

On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no> wrote:

Since rgdal::make_EPSG() is facing the same problems of listing tabulated
EPSG fields as g.proj -l, I was very happy to see Markus’ code in
g.proj/main.c mentioned in this thread, and have used this approach in

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=mar
kup&root=rgdal

However, there are plenty of messages such as: "proj_as_proj_string:
Unsupported conversion method: Lambert Conic Conformal (West

Orientated)". I

haven’t installed GRASS trunk with PROJ6, so I can’t see whether g.proj -l
also sees the same messages.

I get the same messages with g.proj, but no error code from PROJ error
handlers.

If it does, maybe we could ask on the proj list

yes, we should ask. For now g.proj is ignoring those codes for which no
proj string can be retrieved.

I’m not sure what the question is exactly :slight_smile:

My question is, what should GRASS do with a SRS without proj string? The best answer is probably: use WKT instead of a proj string. But most of the GRASS code uses (for historical reasons) some form of proj syntax, also to store projection info for a GRASS location. If there is no proj string, GRASS thinks there is no SRS and falls back to a generic xy system which is not always correct.

Markus M

12

On vendredi 8 mars 2019 17:01:57 CET Markus Metz wrote:

On Fri, Mar 8, 2019 at 12:41 PM Even Rouault <even.rouault@spatialys.com>

wrote:
> On vendredi 8 mars 2019 12:27:37 CET Markus Metz wrote:
> > On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no>

wrote:
> > > Since rgdal::make_EPSG() is facing the same problems of listing

tabulated

> > > EPSG fields as g.proj -l, I was very happy to see Markus' code in
> > > g.proj/main.c mentioned in this thread, and have used this approach in

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=mar
> > kup&root=rgdal
> >
> > > However, there are plenty of messages such as: "proj_as_proj_string:
> > > Unsupported conversion method: Lambert Conic Conformal (West
> >
> > Orientated)". I
> >
> > > haven't installed GRASS trunk with PROJ6, so I can't see whether

g.proj -l

> > > also sees the same messages.
> >
> > I get the same messages with g.proj, but no error code from PROJ error
> > handlers.
> >
> > > If it does, maybe we could ask on the proj list
> >
> > yes, we should ask. For now g.proj is ignoring those codes for which no
> > proj string can be retrieved.
>
> I'm not sure what the question is exactly :slight_smile:

My question is, what should GRASS do with a SRS without proj string? The
best answer is probably: use WKT instead of a proj string. But most of the
GRASS code uses (for historical reasons) some form of proj syntax, also to
store projection info for a GRASS location. If there is no proj string,
GRASS thinks there is no SRS and falls back to a generic xy system which is
not always correct.

The WKT will not help here. You will be able to get here, but you won't be
able to do coordinate transformation with it. This is just that this EPSG
method is not implemented or mapped to a PROJ projection method

--
Spatialys - Geospatial professional services
http://www.spatialys.com

13

I'm assuming there is no error handler to hook onto, to divert errors to
R's error handler.

Actually you can install your own error handler with proj_log_func()

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

14

On Fri, Mar 8, 2019 at 5:56 PM Even Rouault <even.rouault@spatialys.com> wrote:

On vendredi 8 mars 2019 17:01:57 CET Markus Metz wrote:

My question is, what should GRASS do with a SRS without proj string? The
best answer is probably: use WKT instead of a proj string. But most of the
GRASS code uses (for historical reasons) some form of proj syntax, also to
store projection info for a GRASS location. If there is no proj string,
GRASS thinks there is no SRS and falls back to a generic xy system which is
not always correct.

The WKT will not help here. You will be able to get here, but you won’t be
able to do coordinate transformation with it. This is just that this EPSG
method is not implemented or mapped to a PROJ projection method

I understand that coordinate transformation with PROJ will not be possible. But having a valid SRS definition would at least allow to compare two spatial references, i.e. OSRIsSame() should work. Having a valid SRS instead of dumping it would also help if the corresponding PROJ method becomes implemented in the future. That’s why I am wondering about an alternative to the current proj-like projection information in GRASS.

Markus M

15

On vendredi 8 mars 2019 18:22:42 CET Markus Metz wrote:

On Fri, Mar 8, 2019 at 5:56 PM Even Rouault <even.rouault@spatialys.com>

wrote:
> On vendredi 8 mars 2019 17:01:57 CET Markus Metz wrote:
> > My question is, what should GRASS do with a SRS without proj string? The
> > best answer is probably: use WKT instead of a proj string. But most of

the

> > GRASS code uses (for historical reasons) some form of proj syntax, also

to

> > store projection info for a GRASS location. If there is no proj string,
> > GRASS thinks there is no SRS and falls back to a generic xy system

which is

> > not always correct.
>
> The WKT will not help here. You will be able to get here, but you won't be
> able to do coordinate transformation with it. This is just that this EPSG
> method is not implemented or mapped to a PROJ projection method

I understand that coordinate transformation with PROJ will not be possible.
But having a valid SRS definition would at least allow to compare two
spatial references, i.e. OSRIsSame() should work. Having a valid SRS
instead of dumping it would also help if the corresponding PROJ method
becomes implemented in the future. That's why I am wondering about an
alternative to the current proj-like projection information in GRASS.

Yes, for having a reference SRS definition, you can use the WKT export. By
"useful", I meant being able to do coordinate transformation.
To compare SRS for equality, rather than using full string comparison, I'd
suggest using proj_is_equivalent_to() which has a few options for different
levels of equality.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

16

On Fri, 8 Mar 2019, Even Rouault wrote:

On vendredi 8 mars 2019 18:22:42 CET Markus Metz wrote:

On Fri, Mar 8, 2019 at 5:56 PM Even Rouault <even.rouault@spatialys.com>

wrote:

On vendredi 8 mars 2019 17:01:57 CET Markus Metz wrote:

My question is, what should GRASS do with a SRS without proj string? The
best answer is probably: use WKT instead of a proj string. But most of

the

GRASS code uses (for historical reasons) some form of proj syntax, also

to

store projection info for a GRASS location. If there is no proj string,
GRASS thinks there is no SRS and falls back to a generic xy system

which is

not always correct.

The WKT will not help here. You will be able to get here, but you won't be
able to do coordinate transformation with it. This is just that this EPSG
method is not implemented or mapped to a PROJ projection method

I understand that coordinate transformation with PROJ will not be possible.
But having a valid SRS definition would at least allow to compare two
spatial references, i.e. OSRIsSame() should work. Having a valid SRS
instead of dumping it would also help if the corresponding PROJ method
becomes implemented in the future. That's why I am wondering about an
alternative to the current proj-like projection information in GRASS.

Might it be possible for proj_get_crs_info_list_from_database() to return a field indicating whether the record in proj_crs_info can be exported to proj_string?

I found example code for proj_log_func() in GDAL in ogr/ogr_proj_p.cpp, and am using that to do nothing on apparent PJ_LOG_ERROR messages reported in proj_as_proj_string(). The messages are suppressed, and I destroy the context after looping through the projections found. I'm a bit concerned that if something goes wrong, turning off error handling isn't safe.

In rgdal::make_EPSG(), the proj_string is returned as (null), and can be grepped, so it does provide information.

Roger

Yes, for having a reference SRS definition, you can use the WKT export. By
"useful", I meant being able to do coordinate transformation.
To compare SRS for equality, rather than using full string comparison, I'd
suggest using proj_is_equivalent_to() which has a few options for different
levels of equality.

Even

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

17

Dear Roger,

On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no> wrote:

Since rgdal::make_EPSG() is facing the same problems of listing tabulated
EPSG fields as g.proj -l, I was very happy to see Markus’ code in
g.proj/main.c mentioned in this thread, and have used this approach in
https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=markup&root=rgdal

However, there are plenty of messages such as: “proj_as_proj_string:
Unsupported conversion method: Lambert Conic Conformal (West Orientated)”. I
haven’t installed GRASS trunk with PROJ6, so I can’t see whether g.proj -l
also sees the same messages. If it does, maybe we could ask on the proj list
how they might be captured for summary reporting. I think they are coming
from line 5758 in src/iso19111/coordinateoperation.cpp or maybe line 906 in
same file. Maybe PROJ now has an error handler that

Another question concerns the issue of whether one needs to free objects
created, in particular proj_crs_info and pj. Not so important for g.proj,
which exists when done, but important for rgdal whose functions don’t exit.

Anyway, very helpful to see that Markus is looking at the same issues as we
are!

Regarding parsing the traditional PROJ init files, I looked at pj_init.c in PROJ 5 and earlier and attempted to improve parsing those traditional init files in GRASS trunk r74224. I guess that R as well as GRASS is trying to support PROJ versions from 4(.8) onwards.

Regards,

Markus M

18

On Tue, 12 Mar 2019, Markus Metz wrote:

Dear Roger,

On Fri, Mar 8, 2019 at 10:49 AM Roger Bivand <Roger.Bivand@nhh.no> wrote:

Since rgdal::make_EPSG() is facing the same problems of listing tabulated
EPSG fields as g.proj -l, I was very happy to see Markus' code in
g.proj/main.c mentioned in this thread, and have used this approach in

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_info6.cpp?view=markup&root=rgdal

However, there are plenty of messages such as: "proj_as_proj_string:
Unsupported conversion method: Lambert Conic Conformal (West

Orientated)". I

haven't installed GRASS trunk with PROJ6, so I can't see whether g.proj -l
also sees the same messages. If it does, maybe we could ask on the proj

list

how they might be captured for summary reporting. I think they are coming
from line 5758 in src/iso19111/coordinateoperation.cpp or maybe line 906

in

same file. Maybe PROJ now has an error handler that

Another question concerns the issue of whether one needs to free objects
created, in particular proj_crs_info and pj. Not so important for g.proj,
which exists when done, but important for rgdal whose functions don't

exit.

Anyway, very helpful to see that Markus is looking at the same issues as

we

are!

Regarding parsing the traditional PROJ init files, I looked at pj_init.c in
PROJ 5 and earlier and attempted to improve parsing those traditional init
files in GRASS trunk r74224. I guess that R as well as GRASS is trying to
support PROJ versions from 4(.8) onwards.

Thanks for letting me know where to look, I'll add this to our list of hints (Edzer is shepherding the sf package). And yes, we keep getting bitten by our bugs upsetting users of PROJ 4.8.0; Edzer was trying to set up check instances, but there is always something (recently stale proj.pc files reporting the wrong PROJ version).

Best wishes,

Roger

Regards,

Markus M

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en