[gdal-dev] Compound CRS in shapefile
Javier Jimenez Shaw
j1 at jimenezshaw.com
Wed Oct 21 07:21:16 PDT 2020
I forgot a few lines in the C++ code, but they are not needed to reproduce
the problem
if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
{
exit( 1 );
}
.___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__
Entre dos pensamientos racionales
hay infinitos pensamientos irracionales.
On Wed, 21 Oct 2020 at 16:16, Javier Jimenez Shaw <j1 at jimenezshaw.com>
wrote:
> Sorry Even for the lack of detail. Let's see if now it is more clear.
>
> I have done this small code example based on
> https://gdal.org/tutorials/vector_api_tut.html#writing-to-ogr
> but setting a SpatialReference, EPSG:7405 (no checks to make code simpler)
>
> #include "gdal/ogr_spatialref.h"
> #include "gdal/ogrsf_frmts.h"
>
> int main()
> {
> const char *pszDriverName = "ESRI Shapefile";
> GDALDriver *poDriver;
>
> GDALAllRegister();
>
> OGRSpatialReference sr;
> sr.SetFromUserInput("EPSG:7405");
>
> poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );
>
> GDALDataset *poDS;
> poDS = poDriver->Create( "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );
>
> OGRLayer *poLayer;
> poLayer = poDS->CreateLayer( "point_out", &sr, wkbPoint, NULL );
>
> OGRFeature *poFeature;
> poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
> OGRPoint pt;
> pt.setX( 0 );
> pt.setY( 1 );
> pt.setZ( 2 );
> poFeature->SetGeometry( &pt );
> OGRFeature::DestroyFeature( poFeature );
>
> GDALClose( poDS );
> return 0;
> }
>
> 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"],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]]]
>
>
> 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",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]]
>
> 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:
> 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]]
>
> 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.
>
> Thanks
>
> .___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__
> Entre dos pensamientos racionales
> hay infinitos pensamientos irracionales.
>
>
>
> On Tue, 20 Oct 2020 at 19:24, Even Rouault <even.rouault at spatialys.com>
> wrote:
>
>> Javier,
>>
>> > When I create a 3D shapefile with GDAL, and I provide as
>> SpatialReference a
>> > compound CRS, it generates the prj file with the COMPD_CS syntax. Ok so
>> far.
>>
>> Always be specific: show the content of the .prj
>>
>> > However, when I open the shapefile in QGIS, it does not recognize the
>> CRS
>> > at all (neither the horizontal part)
>>
>> What do you call "recognize" exactly: it doesn't find any CRS
>> information, or
>> doesn't retrieve a EPSG code (retrieving a EPSG code from ESRI WKT is
>> deeply
>> dark magic so doesn't always work) ?
>>
>> > Is that intentional? QGIS is using GDAL to deal with shapefiles, afaik.
>>
>> Probably for QGIS mailing list or tracker then.
>>
>> I've tested the following:
>>
>> ogr2ogr out.shp poly.shp -a_srs EPSG:7405
>>
>> ogrinfo out.shp reports correctly a compoundCRS.
>>
>> Opening that in QGIS + GDAL + PROJ master shows EPSG:27700 (the
>> horizontal
>> part). Looking at QGIS QgsCoordinateReferenceSystem class, it only keeps
>> the
>> horizontal part of a CompoundCRS. I guess mostly for compatibility
>> reasons
>> with the rest of the code base.
>>
>> Even
>>
>> --
>> Spatialys - Geospatial professional services
>> http://www.spatialys.com
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201021/9b3cb551/attachment.html>
More information about the gdal-dev
mailing list