[gdal-dev] Compound CRS in shapefile

Javier Jimenez Shaw j1 at jimenezshaw.com
Thu Oct 22 06:41:53 PDT 2020


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 <even.rouault at spatialys.com>
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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201022/72f7cb94/attachment.html>


More information about the gdal-dev mailing list