[PROJ] More confusion regarding proj 6

[Nyall Dawson]

On Tue, 30 Apr 2019 at 20:44, Even Rouault <even.rouault at spatialys.com> 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