[PROJ] More confusion regarding proj 6

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[PROJ] More confusion regarding proj 6

Nyall Dawson
Hi list,

It's time for another set of questions from a mere projection mortal
trying to make way with this:

Q1:

I have the WKT string:

  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"]]

I call proj_create_from_wkt using this string, and get an object back.
But calling methods like proj_get_id_auth_name or proj_get_id_code
gives no result here. I'd expect EPSG/4326. What am I missing?

Q2:

If I call proj_create using "+proj=longlat +ellps=WGS84
+towgs84=0,0,0,0,0,0,0 +no_defs +type=crs", I get something back. What
is this something? Trying to call proj_crs_get_coordinate_system on
the returned value throws the error "Object is not a SingleCRS". How
can I get a coordinate system from this result?

Q3:

Sometimes the objects returned by database creation are CompoundCRS. I
gather I should be using proj_crs_get_sub_crs to get the horizontal
crs from these, but the docs state: "Index of the CRS component
(typically 0 = horizontal, 1 = vertical)". Can I just blindly call use
an index of 0 and hope for the best? I can't see any other methods to
work with CompoundCRS objects, so I can't iterate through the indexes
as I don't know how many there are. Or are there always 2?

That's it for now... there's sure to be more to follow ;)

Nyall
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: More confusion regarding proj 6

Even Rouault-2
Hi Nyall,

>
> Q1:
>
> I have the WKT string:
>
>   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"]]
>
> I call proj_create_from_wkt using this string, and get an object back.
> But calling methods like proj_get_id_auth_name or proj_get_id_code
> gives no result here. I'd expect EPSG/4326. What am I missing?

OK, I admit this is subtle. The issue here is that there is a TOWGS84[]
clause, hence the returned object is not a GeographicCRS, but a BoundCRS of a
GeographicCRS (a rather useless one, because the DATUM is already WGS84, so
the TOWGS84[] clause could have been omitted). The BoundCRS has not the ID
attached to it, only its base CRS.

Let me advertize again the projinfo utility as a investigation tool, because
projinfo 'GEOGCS["WGS...' will return the WKT2_2018 representation, which
shows that it is a BOUNDCRS)

The following will give you what you expect:

#include "proj.h"
#include <assert.h>
#include <stdio.h>

int main()
{
    PJ_CONTEXT* ctx = NULL;
    PJ* p = proj_create_from_wkt(ctx,
        "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\"]]", NULL, NULL, NULL);
   
    assert(proj_get_type(p) == PJ_TYPE_BOUND_CRS);
    PJ* base = proj_get_source_crs(ctx, p);
   
    printf("%s %s\n", proj_get_id_auth_name(base, 0), proj_get_id_code(base,
0));
    proj_destroy(base);
    proj_destroy(p);
    return 0;
}

>
> Q2:
>
> If I call proj_create using "+proj=longlat +ellps=WGS84
> +towgs84=0,0,0,0,0,0,0 +no_defs +type=crs", I get something back. What
> is this something? Trying to call proj_crs_get_coordinate_system on
> the returned value throws the error "Object is not a SingleCRS". How
> can I get a coordinate system from this result?

Same as above. A boundCRS (and projinfo is still of help here :-)) You can get
the coordinate system of its sourceCRS by first fetching it with
proj_get_source_crs().

>
> Q3:
>
> Sometimes the objects returned by database creation are CompoundCRS. I
> gather I should be using proj_crs_get_sub_crs to get the horizontal
> crs from these, but the docs state: "Index of the CRS component
> (typically 0 = horizontal, 1 = vertical)". Can I just blindly call use
> an index of 0 and hope for the best? I can't see any other methods to
> work with CompoundCRS objects, so I can't iterate through the indexes
> as I don't know how many there are. Or are there always 2?

There should be (normally) at least 2 components in a CompoundCRS. In theory
there might be more, but that would be completely exotic objects (spatio-
temporal CRS, spatio-temporal-parametric CRS) and you shouldn't encounter them
in practice.

proj_crs_get_sub_crs(crs, idx) will cleanly return NULL if you try to access a
component that does not exist. Yes, there is not function to return the number
of components, so checking for NULL is the way to go.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: More confusion regarding proj 6

Nyall Dawson
On Tue, 30 Apr 2019 at 20:44, Even Rouault <[hidden email]> wrote:

>
> Hi Nyall,
> >
> > Q1:
> >
> > I have the WKT string:
> >
> >   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"]]
> >
> > I call proj_create_from_wkt using this string, and get an object back.
> > But calling methods like proj_get_id_auth_name or proj_get_id_code
> > gives no result here. I'd expect EPSG/4326. What am I missing?
>
> OK, I admit this is subtle. The issue here is that there is a TOWGS84[]
> clause, hence the returned object is not a GeographicCRS, but a BoundCRS of a
> GeographicCRS (a rather useless one, because the DATUM is already WGS84, so
> the TOWGS84[] clause could have been omitted). The BoundCRS has not the ID
> attached to it, only its base CRS.
>
> Let me advertize again the projinfo utility as a investigation tool, because
> projinfo 'GEOGCS["WGS...' will return the WKT2_2018 representation, which
> shows that it is a BOUNDCRS)
>
> The following will give you what you expect:
>
> #include "proj.h"
> #include <assert.h>
> #include <stdio.h>
>
> int main()
> {
>     PJ_CONTEXT* ctx = NULL;
>     PJ* p = proj_create_from_wkt(ctx,
>         "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\"]]", NULL, NULL, NULL);
>
>     assert(proj_get_type(p) == PJ_TYPE_BOUND_CRS);
>     PJ* base = proj_get_source_crs(ctx, p);
>
>     printf("%s %s\n", proj_get_id_auth_name(base, 0), proj_get_id_code(base,
> 0));
>     proj_destroy(base);
>     proj_destroy(p);
>     return 0;
> }
>
> >
> > Q2:
> >
> > If I call proj_create using "+proj=longlat +ellps=WGS84
> > +towgs84=0,0,0,0,0,0,0 +no_defs +type=crs", I get something back. What
> > is this something? Trying to call proj_crs_get_coordinate_system on
> > the returned value throws the error "Object is not a SingleCRS". How
> > can I get a coordinate system from this result?
>
> Same as above. A boundCRS (and projinfo is still of help here :-)) You can get
> the coordinate system of its sourceCRS by first fetching it with
> proj_get_source_crs().
>
> >
> > Q3:
> >
> > Sometimes the objects returned by database creation are CompoundCRS. I
> > gather I should be using proj_crs_get_sub_crs to get the horizontal
> > crs from these, but the docs state: "Index of the CRS component
> > (typically 0 = horizontal, 1 = vertical)". Can I just blindly call use
> > an index of 0 and hope for the best? I can't see any other methods to
> > work with CompoundCRS objects, so I can't iterate through the indexes
> > as I don't know how many there are. Or are there always 2?
>
> There should be (normally) at least 2 components in a CompoundCRS. In theory
> there might be more, but that would be completely exotic objects (spatio-
> temporal CRS, spatio-temporal-parametric CRS) and you shouldn't encounter them
> in practice.
>
> proj_crs_get_sub_crs(crs, idx) will cleanly return NULL if you try to access a
> component that does not exist. Yes, there is not function to return the number
> of components, so checking for NULL is the way to go.

Perfect -- many thanks Even! Sounds like most of my issues are
stemming from similar underlying causes, and I'll need to wrap up a
generic "getHorizontalSingleCRSfromPJ" function which handles these
different situations.

Nyall

>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj