<div dir="ltr"><div>Oh, I forgot to mention the versions I am using:</div><div>to compile and run the C++ code: PROJ 7.0.1, GDAL 3.1.0 (sorry, that is the environment I had configured).</div><div><br></div><div>To run ogr2ogr unfortunately I forgot to change the path to the build version, and was using the installed one, GDAL 2.2.3. Now repeating with GDAL 3.1.0, ogr2ogr dumps COMPD_CS..., as yours.<br></div><div><br></div><div>About QGIS, I have not compiled it. So I am testing with the installed versions in Ubuntu 18.04 (GDAL 2.2.3) and Ubuntu 20.04 (GDAL 3.0.4, PROJ 6.3.1)</div><div><br></div><div>The prj that a friend gave me was this:</div><div><span style="font-family:monospace">PROJCS["NAD_1983_2011_StatePlane_Colorado_Central_FIPS_0502_Ft_US",GEOGCS["GCS_NAD_1983_2011",DATUM["D_NAD_1983_2011",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",3000000.000316083],PARAMETER["False_Northing",999999.999996],PARAMETER["Central_Meridian",-105.5],PARAMETER["Standard_Parallel_1",38.45],PARAMETER["Standard_Parallel_2",39.75],PARAMETER["Latitude_Of_Origin",37.83333333333334],UNIT["Foot_US",0.3048006096012192]],VERTCS["CGVD2013_height",VDATUM["Canadian_Geodetic_Vertical_Datum_of_2013"],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]</span></div><div><br></div><div>He combined NAD83 with Canadian to not be "the normal one" on purpose. This shapefile is also understood by other software, I guess Civil3D.</div><div><br></div><div>Opening the shp file created with the C++ code (that is with COMPD_CS tag), in ArcGIS ArcMap 10.2.2, gives this info:</div><div><span style="font-family:monospace">Data Type: Shapefile Feature Class <br>Shapefile: C:\Users\fherl\Desktop\point_out\point_out.shp<br>Geometry Type: Point<br>Coordinates have Z values: Yes <br>Coordinates have measures: Yes <br><br>Projected Coordinate System: British_National_Grid<br>Projection: Transverse_Mercator<br>False_Easting: 400000,00000000<br>False_Northing: -100000,00000000<br>Central_Meridian: -2,00000000<br>Scale_Factor: 0,99960127<br>Latitude_Of_Origin: 49,00000000<br>Linear Unit: Meter<br><br>Geographic Coordinate System: GCS_OSGB_1936<br>Datum: D_OSGB_1936<br>Prime Meridian: Greenwich<br>Angular Unit: Degree</span></div><div><br></div><div>Another friend tried in ArcGIS Pro 2.4.3 reading the shapefile created by the C++ code (that is with COMPD_CS tag) and it recognizes it, with both horizontal EPSG:27700 and vertical Newlyn EPSG:5701 CRSs<br></div><div><br></div><div>What I am afraid now is... what is the "proper" format for a compound CRS in the prj file of a shapefile? If GDAL is now exporting it as COMPD_CS, but not everybody understands it, neither the horizontal part, it becomes a compatibility problem.</div><div><br></div><div></div><div>Cheers</div><div>Javier<br></div><div dir="ltr"><div><div><div dir="ltr" data-smartmail="gmail_signature">.___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__<br>Entre dos pensamientos racionales <br>hay infinitos pensamientos irracionales.<br><br></div></div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, 21 Oct 2020 at 18:18, Even Rouault <<a href="mailto:even.rouault@spatialys.com" target="_blank">even.rouault@spatialys.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">> If I try to open the file point_out.shp in QGIS, it does not recognize the<br>
> CRS at all. There is no CRS set to that layer. Not even the horizontal one.<br>
> It also displays a question mark in the "Layer" pane.<br>
> <br>
> The content of the .prj file is this:<br>
> COMPD_CS["OSGB 1936 / British National Grid + ODN<br>
> height",PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_<br>
> 1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0]<br>
> ,UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAME<br>
> TER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETE<br>
> R["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],PARAMETER[<br>
> "Latitude_Of_Origin",49.0],UNIT["Meter",1.0]],VERT_CS["ODN<br>
> height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]]]<br>
> <br>
<br>
Works for me with QGIS, GDAL and PROJ master. Before a yesterday fix in PROJ, <br>
this was identified as EPSG:27700 (only horizontal part), and now this is <br>
identified as EPSG:7405<br>
<br>
> <br>
> Today, inspired by a shapefile created with ArcGIS, I removed the COMPD_CS<br>
> tag, leaving just a PROJCS and a VERT_CS separated by a comma. This is<br>
> recognized by QGIS as EPSG:27700 :<br>
> PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SP<br>
> HEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["D<br>
> egree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["Fal<br>
> se_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Centr<br>
> al_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],PARAMETER["Latitud<br>
> e_Of_Origin",49.0],UNIT["Meter",1.0]],VERT_CS["ODN<br>
> height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]]<br>
<br>
Ah, this rings a bell. PROJ doesn't recognized a PROJCS[],VERT_CS[] as a <br>
compound CRS. It will stop parsing at the end of the PROJCS[]. Perhaps it <br>
should be improvded to recognize the ArcGIS flavor... Could you paste an exact <br>
output of a .prj from ArcGIS (since you say the above one is 'inspired from'), <br>
and the ArcGIS version that outputs it ? And on the reverse, does it like a <br>
.proj with a COMPD_CS like:<br>
<br>
COMPD_CS["OSGB 1936 / British National Grid + ODN <br>
height",PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",<br>
6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",<br>
0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",<br>
400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",<br>
0.9996012717],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",<br>
1.0]],VERTCS["Newlyn",VDATUM["Ordnance_Datum_Newlyn"],PARAMETER["Vertical_Shift",<br>
0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]]<br>
<br>
?<br>
<br>
> <br>
> If I run<br>
> ogr2ogr out.shp point_out.shp -a_srs EPSG:7405<br>
> then it creates the out.prj just with the horizontal CRS, not a compound<br>
> one, neither adding the VERT_CS after the horizontal:<br>
<br>
I get a COMPD_CS with GDAL and PROJ master<br>
<br>
Even<br>
<br>
-- <br>
Spatialys - Geospatial professional services<br>
<a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
</blockquote></div></div>