<div dir="ltr"><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature">I forgot a few lines in the C++ code, but they are not needed to reproduce the problem<br><br><span style="font-family:monospace"> if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )<br> {<br> exit( 1 );<br> }</span><br><br>.___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__<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 16:16, Javier Jimenez Shaw <<a href="mailto:j1@jimenezshaw.com">j1@jimenezshaw.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"><div dir="ltr"><div>Sorry Even for the lack of detail. Let's see if now it is more clear.</div><div><br></div><div>I have done this small code example based on <a href="https://gdal.org/tutorials/vector_api_tut.html#writing-to-ogr" target="_blank">https://gdal.org/tutorials/vector_api_tut.html#writing-to-ogr</a></div><div>but setting a SpatialReference, EPSG:7405 (no checks to make code simpler)<br></div><div><br></div><div><span style="font-family:monospace">#include "gdal/ogr_spatialref.h"<br>#include "gdal/ogrsf_frmts.h"<br></span></div><div><span style="font-family:monospace"><br></span></div><div><span style="font-family:monospace">int main()<br>{<br> const char *pszDriverName = "ESRI Shapefile";<br> GDALDriver *poDriver;<br><br> GDALAllRegister();<br><br> OGRSpatialReference sr;<br> sr.SetFromUserInput("EPSG:7405");<br><br> poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );<br><br> GDALDataset *poDS;<br> poDS = poDriver->Create( "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );<br><br> OGRLayer *poLayer;<br> poLayer = poDS->CreateLayer( "point_out", &sr, wkbPoint, NULL );<br><br> OGRFeature *poFeature;<br> poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );<br> OGRPoint pt;<br> pt.setX( 0 );<br> pt.setY( 1 );</span></div><div><span style="font-family:monospace"> pt.setZ( 2 );<br> poFeature->SetGeometry( &pt );<br> OGRFeature::DestroyFeature( poFeature );<br><br> GDALClose( poDS );<br> return 0;<br>}</span></div><div><br></div><div>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.<br></div><div><br></div><div>The content of the .prj file is this:<br></div><div><span style="font-family:monospace">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]],VERT_CS["ODN height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]]]</span></div><div><br></div><div><br></div><div>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 :<br></div><div></div><div><span style="font-family:monospace">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]],VERT_CS["ODN height",VERT_DATUM["Ordnance Datum Newlyn",2005],UNIT["Meter",1.0]]</span></div><div><br></div><div>If I run</div><div><span style="font-family:monospace">ogr2ogr out.shp point_out.shp -a_srs EPSG:7405</span></div><div>then it creates the out.prj just with the horizontal CRS, not a compound one, neither adding the VERT_CS after the horizontal: <br></div><div><span style="font-family:monospace"></span></div><div><span style="font-family:monospace">PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB 1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]</span></div><div><br></div><div>I was expecting that the method "CreateLayer" was creating a compatible format. However I do not know which one is the proper one. I would not like to lose the vertical CRS information.</div><div><br></div><div>Thanks<br></div><div><br></div><div><div><div dir="ltr">.___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__<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 Tue, 20 Oct 2020 at 19:24, 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">Javier,<br>
<br>
> When I create a 3D shapefile with GDAL, and I provide as SpatialReference a<br>
> compound CRS, it generates the prj file with the COMPD_CS syntax. Ok so far.<br>
<br>
Always be specific: show the content of the .prj<br>
<br>
> However, when I open the shapefile in QGIS, it does not recognize the CRS<br>
> at all (neither the horizontal part)<br>
<br>
What do you call "recognize" exactly: it doesn't find any CRS information, or <br>
doesn't retrieve a EPSG code (retrieving a EPSG code from ESRI WKT is deeply <br>
dark magic so doesn't always work) ?<br>
<br>
> Is that intentional? QGIS is using GDAL to deal with shapefiles, afaik.<br>
<br>
Probably for QGIS mailing list or tracker then.<br>
<br>
I've tested the following:<br>
<br>
ogr2ogr out.shp poly.shp -a_srs EPSG:7405<br>
<br>
ogrinfo out.shp reports correctly a compoundCRS.<br>
<br>
Opening that in QGIS + GDAL + PROJ master shows EPSG:27700 (the horizontal <br>
part). Looking at QGIS QgsCoordinateReferenceSystem class, it only keeps the <br>
horizontal part of a CompoundCRS. I guess mostly for compatibility reasons <br>
with the rest of the code base.<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>
_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a></blockquote></div>
</blockquote></div>