[PROJ] More confusion regarding proj 6

[Even Rouault]

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