cs2cs (5.0.0-rc2) 100% CPU running rdnap tests

[sebastic]

The proj-rdnap Debian package contains a test script which runs cs2cs to test the conversion using the rdtrans2008.gsb & naptrans2008.gtx grids provided by the Dutch kadaster.

Running the test script in verbose mode shows that it hangs on the first test using cs2cs from proj-5.0.0-rc2:

PROJ_LIB=. perl ./testrdtrans2008.pl -v
Direction: From ETRS89 to RD/NAP

Test:   01 Texel
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  53.160753042 4.824761912 42.8614

Debug output of the cs2cs command:

echo "53.160753042 4.824761912 42.8614" | PROJ_LIB=. PROJ_DEBUG=3 cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
get_init: searching cache for key: [epsg:4258]

get_init: searching on in init files for [epsg:4258]

get_init_string: searching for section [4258] in init file [epsg]

pj_open_lib(epsg): call fopen(./epsg) - succeeded

key=epsg:4258, value: [proj=longlat ellps=GRS80 towgs84=0,0,0,0,0,0,0 no_defs]

get_init: got [proj=longlat ellps=GRS80 towgs84=0,0,0,0,0,0,0 no_defs], paralist[0,1]: [proj=longlat,ellps=GRS80]

pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
get_init: searching cache for key: [rdnap:rdnap]

get_init: searching on in init files for [rdnap:rdnap]

get_init_string: searching for section [rdnap] in init file [rdnap]

pj_open_lib(rdnap): call fopen(./rdnap) - succeeded

key=rdnap:rdnap, value: [proj=sterea lat_0=52.15616055555555 lon_0=5.38763888888889 k=0.9999079 x_0=155000 y_0=463000 ellps=bessel nadgrids=rdtrans2008.gsb geoidgrids=naptrans2008.gtx units=m no_defs]

get_init: got [proj=sterea lat_0=52.15616055555555 lon_0=5.38763888888889 k=0.9999079 x_0=155000 y_0=463000 ellps=bessel nadgrids=rdtrans2008.gsb geoidgrids=naptrans2008.gtx units=m no_defs], paralist[0,1]: [proj=sterea,lat_0=52.15616055555555]

pj_ellipsoid - final: a=6377397.155 f=1/299.153, errno=0
pj_ellipsoid - final:    ellps=bessel
pj_open_lib(rdtrans2008.gsb): call fopen(./rdtrans2008.gsb) - succeeded

NTv2 NL_ALL   63x65: LL=(2.5,50.5) UR=(7.66666667,55.8333333)

NTv2 NL_LAND  561x421: LL=(3,50.5) UR=(7.66666667,54)

pj_open_lib(rdtrans2008.gsb): call fopen(./rdtrans2008.gsb) - succeeded

NTv2 NL_ALL   63x65: LL=(2.5,50.5) UR=(7.66666667,55.8333333)

NTv2 NL_LAND  561x421: LL=(3,50.5) UR=(7.66666667,54)

[kbevers]

This must be related to the NTv2 file having sub-grids. The grid loading and reading code has been refactored quite a lot and the sub-grid features were not tested i think. Where can I get a copy of this grid?

[sebastic]

Where can I get a copy of this grid?

From the Debian package git repository, the files are linked in the OP as is the patch that creates the test script.

[kbevers]

Sorry, totally missed that. Thanks. I'll have a look at it.

[kbevers]

I've looked into this now. I got a segfault when trying to reproduce the above. That is fixed in #792. Hopefully it also fixes this issue.

[sebastic]

I've included 332bf64 as a patch in the Debian package, and confirmed that fixes the issue with the proj-rdnap test script.

The expected output for coordinates outside the grid has also been updated to support inf for height.

[busstoptaktik]

The expected output for coordinates outside the grid has also been updated to support inf for height

@sebastic: I'm missing context here, so just responding "just in case": Please make sure you check against HUGE_VAL and not whatever goes for infinity on whatever platform: HUGE_VAL equals DBL_MAX on platforms not supporting proper infinities.

[sebastic]

There are a couple of test that have coordinates outside of the grid.

From Use of RDTRANS2008 and NAPTRANS2008.pdf:

Points 07 - 10 are outside the region where interpolation between either the NLGEO2004 geoid or the RD correction grid points is possible. RD is defined only within the region enclosed by the following points (in RD), outside this region RD coordinates can be computed, but the output should be handled with care.

Corners of the validity region for RD:

┌────────┬────────┐
│  x (m) │  y (m) │
├────────┼────────┤
│ 140000 │ 630000 │
│ 100000 │ 600000 │
│  80000 │ 500000 │
│  -8000 │ 390000 │
│  -8000 │ 335000 │
│ 100000 │ 335000 │
│ 160000 │ 288000 │
│ 220000 │ 288000 │
│ 301000 │ 450000 │
│ 301000 │ 615000 │
│ 260000 │ 630000 │
└────────┴────────┘

To not fail the build due to these tests failing to report coordinates there is special handling for the output of cs2cs for these tests, a regular expression is applied to the height which if it matches lets the test pass.

With proj 4.9.x the output was '*' for lat/lon and a (negative) floating point value for the height, now inf is reported for the height instead.

Verbose output of the test script for some context:

# PROJ_LIB=. perl testrdtrans2008.pl -v                          
Direction: From ETRS89 to RD/NAP

Test:   01 Texel
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  53.160753042 4.824761912 42.8614
Output: 117380.1203     575040.3398 1.0000
Expect: 117380.1200     575040.3400 1.0000
Test OK: From ETRS89 to RD/NAP - 01 Texel (Not identical, but within margin)

Test:   02 Noord-Groningen
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  53.419482050 6.776726674 42.3586
Output: 247380.5600     604580.7798 1.9999
Expect: 247380.5600     604580.7800 2.0000
Test OK: From ETRS89 to RD/NAP - 02 Noord-Groningen (Not identical, but within margin)

Test:   03 Amersfoort
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  52.155172897 5.387203657 43.2551
Output: 155000.0000     462999.9998 -0.0000
Expect: 155000.0000     463000.0000 0.0000
Test OK: From ETRS89 to RD/NAP - 03 Amersfoort (Not identical, but within margin)

Test:   04 Zeeuws-Vlaanderen
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  51.368607152 3.397588595 47.4024
Output: 16460.9098      377380.2298 3.0000
Expect: 16460.9100      377380.2300 3.0000
Test OK: From ETRS89 to RD/NAP - 04 Zeeuws-Vlaanderen (Not identical, but within margin)

Test:   05 Zuid-Limburg
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  50.792584908 5.773795547 174.9478
Output: 182260.4500     311480.6701 129.0000
Expect: 182260.4500     311480.6700 129.000
Test OK: From ETRS89 to RD/NAP - 05 Zuid-Limburg (Not identical, but within margin)

Test:   06 Maasvlakte
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  51.947393898 4.072887101 47.5968
Output: 64640.8899      440700.0099 3.9998
Expect: 64640.8900      440700.0101 4.0000
Test OK: From ETRS89 to RD/NAP - 06 Maasvlakte (Not identical, but within margin)

Test:   07 No_rd&geoid
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  50.687420405 4.608971812
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * -inf
Expect: 100000.6700     300000.8900
Test OK: From ETRS89 to RD/NAP - 07 No_rd&geoid (Expected output overriden: *   * ^-?(\d+\.\d+|inf)$)

Test:   08 No_geoid
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  51.136825197 4.601375361
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * -inf
Expect: 100000.6700     350000.8900
Test OK: From ETRS89 to RD/NAP - 08 No_geoid (Expected output overriden: *      * ^-?(\d+\.\d+|inf)$)

Test:   09 No_rd
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  52.482440839 4.268403889
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * -inf
Expect: 79000.0100      500000.2300
Test OK: From ETRS89 to RD/NAP - 09 No_rd (Expected output overriden: * * ^-?(\d+\.\d+|inf)$)

Test:   10 edge_rd
Exec:   cs2cs -r +init=epsg:4258 +to +init=rdnap:rdnap -f %.4f
Input:  51.003976532 3.891247830
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * -inf
Expect: 50000.4500      335999.6700
Test OK: From ETRS89 to RD/NAP - 10 edge_rd (Expected output overriden: *       * ^-?(\d+\.\d+|inf)$)

Direction: From RD/NAP to ETRS89

Test:   01 Texel
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  117380.1200 575040.3400 1.0000
Output: 53.160753044    4.824761908 42.861397591
Expect: 53.160753042    4.824761912 42.8614
Test OK: From RD/NAP to ETRS89 - 01 Texel (Not identical, but within margin)

Test:   02 Noord-Groningen
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  247380.5600 604580.7800 2.0000
Output: 53.419482052    6.776726674 42.358705049
Expect: 53.419482050    6.776726674 42.3586
Test OK: From RD/NAP to ETRS89 - 02 Noord-Groningen (Not identical, but within margin)

Test:   03 Amersfoort
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  155000.0000 463000.0000 0.0000
Output: 52.155172899    5.387203658 43.255106924
Expect: 52.155172897    5.387203657 43.2551
Test OK: From RD/NAP to ETRS89 - 03 Amersfoort (Not identical, but within margin)

Test:   04 Zeeuws-Vlaanderen
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  16460.9100 377380.2300 3.0000
Output: 51.368607153    3.397588597 47.402397195
Expect: 51.368607152    3.397588595 47.4024
Test OK: From RD/NAP to ETRS89 - 04 Zeeuws-Vlaanderen (Not identical, but within margin)

Test:   05 Zuid-Limburg
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  182260.4500 311480.6700 129.000
Output: 50.792584907    5.773795547 174.947778837
Expect: 50.792584908    5.773795547 174.9478
Test OK: From RD/NAP to ETRS89 - 05 Zuid-Limburg (Not identical, but within margin)

Test:   06 Maasvlakte
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  64640.8900 440700.0101 4.0000
Output: 51.947393899    4.072887103 47.596976393
Expect: 51.947393898    4.072887101 47.5968
Test OK: From RD/NAP to ETRS89 - 06 Maasvlakte (Not identical, but within margin)

Test:   07 No_rd&geoid
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  100000.6700 300000.8900
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * inf
Expect: 50.687420405    4.608971812
Test OK: From RD/NAP to ETRS89 - 07 No_rd&geoid (Expected output overriden: *   * ^-?(\d+\.\d+|inf)$)

Test:   08 No_geoid
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  100000.6700 350000.8900
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * inf
Expect: 51.136825197    4.601375361
Test OK: From RD/NAP to ETRS89 - 08 No_geoid (Expected output overriden: *      * ^-?(\d+\.\d+|inf)$)

Test:   09 No_rd
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  79000.0100 500000.2300
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * inf
Expect: 52.482440839    4.268403889
Test OK: From RD/NAP to ETRS89 - 09 No_rd (Expected output overriden: * * ^-?(\d+\.\d+|inf)$)

Test:   10 edge_rd
Exec:   cs2cs -s +init=rdnap:rdnap +to +init=epsg:4258 -f %.9f
Input:  50000.4500 335999.6700
STDERR: Rel. 5.0.0, 15 february 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): point not within available datum shift grids
Output: *       * inf
Expect: 51.003976532    3.891247830
Test OK: From RD/NAP to ETRS89 - 10 edge_rd (Expected output overriden: *       * ^-?(\d+\.\d+|inf)$)