[gdal-dev] Compound CRS in shapefile

[Javier Jimenez Shaw]

Sorry Even for the lack of detail. Let's see if now it is more clear.

I have done this small code example based on
https://gdal.org/tutorials/vector_api_tut.html#writing-to-ogr
but setting a SpatialReference, EPSG:7405 (no checks to make code simpler)

#include "gdal/ogr_spatialref.h"
#include "gdal/ogrsf_frmts.h"

int main()
{
    const char *pszDriverName = "ESRI Shapefile";
    GDALDriver *poDriver;

    GDALAllRegister();

    OGRSpatialReference sr;
    sr.SetFromUserInput("EPSG:7405");

    poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );

    GDALDataset *poDS;
    poDS = poDriver->Create( "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );

    OGRLayer *poLayer;
    poLayer = poDS->CreateLayer( "point_out", &sr, wkbPoint, NULL );

    OGRFeature *poFeature;
    poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
    OGRPoint pt;
    pt.setX( 0 );
    pt.setY( 1 );
    pt.setZ( 2 );
    poFeature->SetGeometry( &pt );
    OGRFeature::DestroyFeature( poFeature );

    GDALClose( poDS );
    return 0;
}

If I try to open the file point_out.shp in QGIS, it does not recognize the
CRS at all. There is no CRS set to that layer. Not even the horizontal one.
It also displays a question mark in the "Layer" pane.

The content of the .prj file is this:
COMPD_CS["OSGB 1936 / British National Grid + ODN
height",PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]],VERT_CS["ODN
height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]]]


Today, inspired by a shapefile created with ArcGIS, I removed the COMPD_CS
tag, leaving just a PROJCS and a VERT_CS separated by a comma. This is
recognized by QGIS as EPSG:27700 :
PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]],VERT_CS["ODN
height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]]

If I run
ogr2ogr out.shp point_out.shp -a_srs EPSG:7405
then it creates the out.prj just with the horizontal CRS, not a compound
one, neither adding the VERT_CS after the horizontal:
PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB
1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

I was expecting that the method "CreateLayer" was creating a compatible
format. However I do not know which one is the proper one. I would not like
to lose the vertical CRS information.

Thanks

.___ ._ ..._ .. . ._.  .___ .. __ . _. . __..  ... .... ._ .__
Entre dos pensamientos racionales
hay infinitos pensamientos irracionales.



On Tue, 20 Oct 2020 at 19:24, Even Rouault <even.rouault at spatialys.com>
wrote:

> Javier,
>
> > When I create a 3D shapefile with GDAL, and I provide as
> SpatialReference a
> > compound CRS, it generates the prj file with the COMPD_CS syntax. Ok so
> far.
>
> Always be specific: show the content of the .prj
>
> > However, when I open the shapefile in QGIS, it does not recognize the CRS
> > at all (neither the horizontal part)
>
> What do you call "recognize" exactly: it doesn't find any CRS information,
> or
> doesn't retrieve a EPSG code (retrieving a EPSG code from ESRI WKT is
> deeply
> dark magic so doesn't always work) ?
>
> > Is that intentional? QGIS is using GDAL to deal with shapefiles, afaik.
>
> Probably for QGIS mailing list or tracker then.
>
> I've tested the following:
>
> ogr2ogr out.shp poly.shp -a_srs EPSG:7405
>
> ogrinfo out.shp reports correctly a compoundCRS.
>
> Opening that in QGIS + GDAL + PROJ master shows EPSG:27700 (the horizontal
> part). Looking at QGIS QgsCoordinateReferenceSystem class, it only keeps
> the
> horizontal part of a CompoundCRS. I guess mostly for compatibility reasons
> with the rest of the code base.
>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201021/ad72837d/attachment-0001.html>