Oh, I forgot to mention the versions I am using: to compile and run the C++ code: PROJ 7.0.1, GDAL 3.1.0 (sorry, that is the environment I had configured).
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. 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) The prj that a friend gave me was this: 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]] He combined NAD83 with Canadian to not be "the normal one" on purpose. This shapefile is also understood by other software, I guess Civil3D. Opening the shp file created with the C++ code (that is with COMPD_CS tag), in ArcGIS ArcMap 10.2.2, gives this info: Data Type: Shapefile Feature Class Shapefile: C:\Users\fherl\Desktop\point_out\point_out.shp Geometry Type: Point Coordinates have Z values: Yes Coordinates have measures: Yes Projected Coordinate System: British_National_Grid Projection: Transverse_Mercator False_Easting: 400000,00000000 False_Northing: -100000,00000000 Central_Meridian: -2,00000000 Scale_Factor: 0,99960127 Latitude_Of_Origin: 49,00000000 Linear Unit: Meter Geographic Coordinate System: GCS_OSGB_1936 Datum: D_OSGB_1936 Prime Meridian: Greenwich Angular Unit: Degree 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 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. Cheers Javier .___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__ Entre dos pensamientos racionales hay infinitos pensamientos irracionales. On Wed, 21 Oct 2020 at 18:18, Even Rouault <[email protected]> wrote: > > 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"],PARAME > > > TER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETE > > > R["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]]] > > > > Works for me with QGIS, GDAL and PROJ master. Before a yesterday fix in > PROJ, > this was identified as EPSG:27700 (only horizontal part), and now this is > identified as EPSG:7405 > > > > > 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",SP > > > HEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["D > > > egree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["Fal > > > se_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Centr > > > al_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],PARAMETER["Latitud > > e_Of_Origin",49.0],UNIT["Meter",1.0]],VERT_CS["ODN > > height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]] > > Ah, this rings a bell. PROJ doesn't recognized a PROJCS[],VERT_CS[] as a > compound CRS. It will stop parsing at the end of the PROJCS[]. Perhaps it > should be improvded to recognize the ArcGIS flavor... Could you paste an > exact > output of a .prj from ArcGIS (since you say the above one is 'inspired > from'), > and the ArcGIS version that outputs it ? And on the reverse, does it like > a > .proj with a COMPD_CS like: > > 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]],VERTCS["Newlyn",VDATUM["Ordnance_Datum_Newlyn"],PARAMETER["Vertical_Shift", > 0.0],PARAMETER["Direction",1.0],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: > > I get a COMPD_CS with GDAL and PROJ master > > Even > > -- > Spatialys - Geospatial professional services > http://www.spatialys.com >
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
