[gdal-dev] BAG CRS

Even Rouault even.rouault at spatialys.com
Fri Apr 22 09:47:17 PDT 2022 

André,

hum ok. it might perhaps be possible that WKT2 could come when extract 
the vertical part but I still can't reproduce that trying:

$ gdal_translate autotest/gcore/data/byte.tif byte.bag -a_srs 
'COMPOUNDCRS["foo",PROJCRS["NAD27 / UTM zone 
11N",BASEGEOGCRS["NAD27",DATUM["North American Datum 
1927",ELLIPSOID["Clarke 
1866",6378206.4,294.978698213898,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4267]],CONVERSION["UTM 
zone 11N",METHOD["Transverse 
Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural 
origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude 
of natural 
origin",-117,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale 
factor at natural 
origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False 
easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False 
northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering 
survey, topographic mapping."],AREA["North America - between 120°W and 
114°W - onshore. Canada - Alberta; British Columbia; Northwest 
Territories; Nunavut. Mexico. United States (USA) - California; Idaho; 
Nevada; Oregon; 
Washington."],BBOX[26.93,-120,78.13,-114]],ID["EPSG",26711]],VERTCRS["ellipse",VDATUM["NAD83(2011) 
/ UTM zone 10N + ellipse"],CS[vertical,1],AXIS["ellipsoid height 
(h)",up,LENGTHUNIT["metre",1]]]]'

$ gdalinfo byte.bag -mdd xml:BAG | grep VERT_CS
<gco:CharacterString>VERT_CS["ellipse",VERT_DATUM["NAD83(2011) / UTM 
zone 10N + ellipse",2005],UNIT["metre",1],AXIS["Ellipsoid 
height",UP]]</gco:CharacterString>

(the ,2005 here isn't really appropriate. it should rather be ,2002,  
but such formulations of CompoundCRS with a VerticalCRS that expressed 
an ellipsoidal height are not ISO-19111 compliant)

Anyway if you can find a way to reproduce the issue with a CRS provided 
through the GetProjection()/SetProjection() API of GDAL, a fix would be 
welcome. If it is the user who provides WKT2 with a manual VAR_VERT_WKT 
option, then I'd say it is their responsibility to provide the right value.

Even

Le 22/04/2022 à 18:35, Vautour, André (INT) a écrit :
>
> Sorry Even,
>
> I should have done more research on my end up front before sending 
> this. A client told us his WKT 2.0 BAG was created with GDAL. Looking 
> at it some more, it was the vertical CRS that was WKT 2.0.
>
> Their vertical CRS looked as follows:
>
> <gco:CharacterString>VERTCRS["ellipse",VDATUM["NAD83(2011) / UTM zone 
> 10N + ellipse"],CS[vertical,1],AXIS["ellipsoid height 
> (h)",up,LENGTHUNIT["metre",1]]]</gco:CharacterString>
>
> Cheers,
>
> André
>
> *From:* Even Rouault <even.rouault at spatialys.com>
> *Sent:* April 22, 2022 13:22
> *To:* Vautour, André (INT) <Andre.Vautour at Teledyne.com>; 
> gdal-dev at lists.osgeo.org
> *Subject:* Re: [gdal-dev] BAG CRS
>
> External Email
>
> André,
>
> I don't confirm this trying:
>
> $ gdal_translate autotest/gcore/data/byte.tif byte.bag
>
> $ gdalinfo byte.bag -mdd xml:BAG | grep PROJCS
>             <gco:CharacterString>PROJCS["NAD27 / UTM zone 
> 11N",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 
> 1866",6378206.4,294.978698213898,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26711"]]</gco:CharacterString>
>
> This is WKT1
>
> I'm wondering if the confusion might come from looking at the CRS 
> reported by gdalinfo, which will always be WKT2 in recent GDAL 
> versions whatever the encoding in the source dataset itself ?
>
> Or perhaps you try to write a Geographic 3D CRS ? in which case WKT2 
> will be used since there's no WKT1 representation for this
>
> $ gdalwarp -overwrite autotest/gcore/data/byte.tif byte.bag -t_srs 
> EPSG:4979 -ot float32
>
> $ gdalinfo byte.bag -mdd xml:BAG | grep GEOGCRS
>             <gco:CharacterString>GEOGCRS["WGS 84",ENSEMBLE["World 
> Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 
> (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World 
> Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 
> (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World 
> Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 
> (G2139)"],ELLIPSOID["WGS 
> 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["geodetic 
> latitude 
> (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic 
> longitude 
> (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal 
> height (h)",up,ORDER[3],LENGTHUNIT["metre",1]], [... snip ....] 
> ,ID["EPSG",4979]]</gco:CharacterString>
>
> (I see that gdalinfo also issues a "ERROR 1: PROJ: 
> proj_create_compound_crs: components of the compound CRS do not belong 
> to one of the allowed combinations of 
> http://docs.opengeospatial.org/as/18-005r4/18-005r4.html#34" when 
> reading this, which I suspect comes from the fact that it tries to 
> create a CompoundCRS from this 3D GeogCRS and the "VERT_CS["unknown", 
> VERT_DATUM["unknown", 2000]]", so writing 3D GeogCRS isn't really 
> appropriate for BAG)
>
> Even
>
> Le 22/04/2022 à 18:01, Vautour, André (INT) a écrit :
>
>     Hi all,
>
>     It has just come to my attention that the BAG driver in GDAL is
>     writing the CRS with a WKT codeSpace and a WKT 2.0 string as the
>     code. While I am fairly sure we started that incorrect practice
>     here at CARIS with a WKT 1.0 string, this is the first time I’ve
>     seen a WKT 2.0 string written in there, and I am therefore
>     concerned about interoperability.
>
>     I think we can all agree that what makes the most sense is to
>     write an EPSG codebase with the EPSG code if we have it, and I
>     think falling back to a WKT 1.0 string if it is not EPSG is a
>     reasonable fallback if it is not EPSG, given how it has been used
>     in the past, but I am worried that older software won’t recognize
>     the WKT 2.0 string.
>
>     Would you be willing to take a pull request that does that?
>     Worst-case, would you be fine with changing that WKT 2.0 string to
>     a WKT 1.0 string?
>
>     Regards,
>
>     André
>
>
>
>     _______________________________________________
>
>     gdal-dev mailing list
>
>     gdal-dev at lists.osgeo.org
>
>     https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- 
> http://www.spatialys.com
> My software is free, but my time generally not.

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20220422/6a8b6f7b/attachment.html>

More information about the gdal-dev mailing list