compound cs with vertical datum fails to validate and transform

[alon]

Hi,

I am using GDAL 1.11.0 compiled with proj 4.8.0 under windows. I have PROJ_LIB defined. and egm86_15.gtx file inside the directory.

this could be a proj library problem. if this is the case, let me know.

consider the following scenario:

OGRSpatialReference srs1; srs1.importFromEPSG(4326); srs1.SetVertCS("EGM96 geoid height", "EGM96 geoid"); srs1.SetExtension("VERT_DATUM", "PROJ4_GRIDS", "egm96_15.gtx");

OGRErr srs1v = srs1.Validate();

srs1v returns 5 (OGRErr srs1v = srs1.Validate())

when I export this to well known text, I get this.

COMPD_CS[

GEOGCS["WGS 84",

DATUM["WGS_1984",

SPHEROID["WGS 84",6378137,298.257223563,

AUTHORITY["EPSG","7030"]],

AUTHORITY["EPSG","6326"]],

PRIMEM["Greenwich",0,

AUTHORITY["EPSG","8901"]],

UNIT["degree",0.0174532925199433,

AUTHORITY["EPSG","9122"]],

AUTHORITY["EPSG","4326"]],

VERT_CS["EGM96 geoid height",

VERT_DATUM["EGM96 geoid",2005,

EXTENSION["PROJ4_GRIDS","egm96_15.gtx"]],

AXIS["Up",UP]]]

however:

OGRSpatialReference srs3; srs3.SetFromUserInput("epsg:4326+5773");

OGRErr srs2v = srs2.Validate();

is validated OK (0)

COMPD_CS["WGS 84 + EGM96 geoid height",

GEOGCS["WGS 84",

DATUM["WGS_1984",

SPHEROID["WGS 84",6378137,298.257223563,

AUTHORITY["EPSG","7030"]],

AUTHORITY["EPSG","6326"]],

PRIMEM["Greenwich",0,

AUTHORITY["EPSG","8901"]],

UNIT["degree",0.0174532925199433,

AUTHORITY["EPSG","9122"]],

AUTHORITY["EPSG","4326"]],

VERT_CS["EGM96 geoid height",

VERT_DATUM["EGM96 geoid",2005,

AUTHORITY["EPSG","5171"], EXTENSION["PROJ4_GRIDS","egm96_15.gtx"]],

UNIT["metre",1,

AUTHORITY["EPSG","9001"]],

AXIS["Up",UP], AUTHORITY["EPSG","5773"]]]

the coordiante system is actually identical. the only difference is EPSG authority. why the first one fails?

now I want to "read" the vetical datum ( in lat long, WGS84) to do that I transfrom from 4326 to this compound CS.

OGRSpatialReference srs; srs.importFromEPSG(4326); OGRCoordinateTransformation *oct = OGRCreateCoordinateTransformation(&srs, &srs2); it can be srs1 or srs2, same results.

when I transform all the coordinates from -180 to 180, -90 to 90 in quarter of a degree step, I get 2159 calls succeed and 1034641 fails. I get the z value even when call fails. the return x and y values are converted to radians. why is that? why the call sometime fails while returning results that seems accurate?

[Even Rouault]

When enabling debug traces, the reason for the failure appears clearly : "OGRSpatialReference::Validate: No UNIT child in VERT_CS."

Regarding transforms, I get satisfactory results. See below. Failures can happen if you hit the borders (longitude=-180 / 180, latitude = -90 / 90), but otherwise works fine.

$ PROJ_LIB=$PWD gdaltransform --debug on -s_srs EPSG:4326 -t_srs EPSG:4326+5773
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs 
OGRCT: Target: +proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +vunits=m +no_defs 
OGRCT: Source: +proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +vunits=m +no_defs 
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs 

2 49
2 49 -44.6434936523438

-179 89
-179 89 -12.8665313720703

179 89
179 89 -12.8683443069458

-100 50
-100 50 23.7892208099366

[alon]

Hi,

Thanks for the quick reply. as for the missing units, my mistake.

my second mistake is that I assumed Transform retuns OGRErr and not BOOL.

I noticed that I get errors also when the long value is 179.75. is there a way to fix this?

Regards, Alon.

[Even Rouault]

Regarding erros related to long value 179.75, running gdalinfo on the gtx files shows :

Origin = (-180.125000000000000,90.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Corner Coordinates:
Upper Left  (-180.1250000,  90.1250000) (180d 7'30.00"W, 90d 7'30.00"N)
Lower Left  (-180.1250000, -90.1250000) (180d 7'30.00"W, 90d 7'30.00"S)
Upper Right ( 179.8750000,  90.1250000) (179d52'30.00"E, 90d 7'30.00"N)
Lower Right ( 179.8750000, -90.1250000) (179d52'30.00"E, 90d 7'30.00"S)
Center      (  -0.1250000,   0.0000000) (  0d 7'30.00"W,  0d 0' 0.01"N)

So 179.75 would seem to be in validity area of the grid. But proj doesn't use GDAL to read georeferencing, so I'm not sure if takes the same conventions regarding pixel-is-area vs pixel-is-center. In any case, this would be a proj.4 issue to report in Proj bug tracker ( ​http://trac.osgeo.org/proj/wiki ). Closing this as a non GDAL bug.

Note: similar proj.4 ticket loosely related ​http://trac.osgeo.org/proj/ticket/126