[gdal-dev] Compound CRS in shapefile

[Even Rouault]

> 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