> 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
