Skip to content

Commit

Permalink
Include US VERTCON
Browse files Browse the repository at this point in the history
  • Loading branch information
dvdoug committed Oct 7, 2021
1 parent 7cffb98 commit b7cee81
Show file tree
Hide file tree
Showing 10 changed files with 100 additions and 6 deletions.
Binary file added resources/papers/NOAA_TR_NOS_NGS_0068.pdf
Binary file not shown.
24 changes: 24 additions & 0 deletions src/CoordinateOperation/CRSTransformations.php
Original file line number Diff line number Diff line change
Expand Up @@ -55143,6 +55143,30 @@ class CRSTransformations
'accuracy' => 0.1,
'reversible' => true,
],
[
'operation' => 'urn:ogc:def:coordinateOperation:EPSG::7969',
'name' => 'NGVD29 height (m) to NAVD88 height (1)',
'source_crs' => 'urn:ogc:def:crs:EPSG::7968',
'target_crs' => 'urn:ogc:def:crs:EPSG::5703',
'accuracy' => 0.02,
'reversible' => true,
],
[
'operation' => 'urn:ogc:def:coordinateOperation:EPSG::7970',
'name' => 'NGVD29 height (m) to NAVD88 height (2)',
'source_crs' => 'urn:ogc:def:crs:EPSG::7968',
'target_crs' => 'urn:ogc:def:crs:EPSG::5703',
'accuracy' => 0.02,
'reversible' => true,
],
[
'operation' => 'urn:ogc:def:coordinateOperation:EPSG::7971',
'name' => 'NGVD29 height (m) to NAVD88 height (3)',
'source_crs' => 'urn:ogc:def:crs:EPSG::7968',
'target_crs' => 'urn:ogc:def:crs:EPSG::5703',
'accuracy' => 0.02,
'reversible' => true,
],
[
'operation' => 'urn:ogc:def:coordinateOperation:EPSG::7977',
'name' => 'HKPD depth to HKCD depth (1)',
Expand Down
7 changes: 7 additions & 0 deletions src/CoordinateOperation/CoordinateOperationMethods.php
Original file line number Diff line number Diff line change
Expand Up @@ -753,6 +753,13 @@ class CoordinateOperationMethods
*/
public const EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NZLVD = 'urn:ogc:def:method:EPSG::1071';

/**
* Vertical Offset by Grid Interpolation (VERTCON)
* Any NAD realization may be used as the Interpolation CRS; bi-linear interpolation is used. Input expects
* longitudes to be positive west.
*/
public const EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_VERTCON = 'urn:ogc:def:method:EPSG::9658';

/**
* Vertical Offset by Grid Interpolation (gtx).
*/
Expand Down
18 changes: 18 additions & 0 deletions src/CoordinateOperation/CoordinateOperationParams.php
Original file line number Diff line number Diff line change
Expand Up @@ -47834,6 +47834,24 @@ class CoordinateOperationParams
'reverses' => true,
],
],
'urn:ogc:def:coordinateOperation:EPSG::7969' => [
'verticalOffsetFile' => [
'reverses' => true,
'fileProvider' => 'PHPCoord\\CoordinateOperation\\GTXNGVD29NAVD88CONUSWestProvider',
],
],
'urn:ogc:def:coordinateOperation:EPSG::7970' => [
'verticalOffsetFile' => [
'reverses' => true,
'fileProvider' => 'PHPCoord\\CoordinateOperation\\GTXNGVD29NAVD88CONUSCentralProvider',
],
],
'urn:ogc:def:coordinateOperation:EPSG::7971' => [
'verticalOffsetFile' => [
'reverses' => true,
'fileProvider' => 'PHPCoord\\CoordinateOperation\\GTXNGVD29NAVD88CONUSEastProvider',
],
],
'urn:ogc:def:coordinateOperation:EPSG::7977' => [
'verticalOffset' => [
'value' => -0.146,
Expand Down
15 changes: 15 additions & 0 deletions src/CoordinateOperation/CoordinateOperations.php
Original file line number Diff line number Diff line change
Expand Up @@ -8840,6 +8840,21 @@ class CoordinateOperations
'method' => 'urn:ogc:def:method:EPSG::9616',
'extent_code' => ['2530'],
],
'urn:ogc:def:coordinateOperation:EPSG::7969' => [
'name' => 'NGVD29 height (m) to NAVD88 height (1)',
'method' => 'urn:ogc:def:method:EPSG::1084',
'extent_code' => ['2950'],
],
'urn:ogc:def:coordinateOperation:EPSG::7970' => [
'name' => 'NGVD29 height (m) to NAVD88 height (2)',
'method' => 'urn:ogc:def:method:EPSG::1084',
'extent_code' => ['2949'],
],
'urn:ogc:def:coordinateOperation:EPSG::7971' => [
'name' => 'NGVD29 height (m) to NAVD88 height (3)',
'method' => 'urn:ogc:def:method:EPSG::1084',
'extent_code' => ['2948'],
],
'urn:ogc:def:coordinateOperation:EPSG::7977' => [
'name' => 'HKPD depth to HKCD depth (1)',
'method' => 'urn:ogc:def:method:EPSG::9616',
Expand Down
5 changes: 5 additions & 0 deletions src/CoordinateOperation/GTXGrid.php
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,11 @@ public function getAdjustment(GeographicPoint $point): Metre
$longitude = $point->getLongitude()->getValue();
$offset = $this->interpolateBilinear($longitude, $latitude)[0];

// These are in millimeters for some reason... :/
if (in_array($this->getBasename(), ['vertconc.gtx', 'vertcone.gtx', 'vertconw.gtx'], true)) {
$offset /= 1000;
}

return new Metre($offset);
}

Expand Down
1 change: 0 additions & 1 deletion src/EPSG/Import/EPSGCodegenFromDataImport.php
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,6 @@ class EPSGCodegenFromDataImport
1106, // Geographic3D to GravityRelatedHeight (ITAL2005)
9620, // Norway Offshore Interpolation
9634, // Maritime Provinces polynomial interpolation
9658, // Vertical Offset by Grid Interpolation (VERTCON)

// only distributed as .dll, can't use
1036, // Cartesian Grid Offsets from Form Function
Expand Down
19 changes: 14 additions & 5 deletions src/EPSG/Import/EPSGImporter.php
Original file line number Diff line number Diff line change
Expand Up @@ -70,12 +70,21 @@ public function createSQLiteDB(): void
* AusGeoidv2 abuses the NTv2 file format to have very large files, and a confusing implementation (latitude shifts are actually height offsets),
* so use a GTX conversion of those files instead...
*/
$sqlite->exec('UPDATE epsg_coordoperation SET coord_op_method_code = 1088 WHERE coord_op_method_code = 1083');
$sqlite->exec('UPDATE epsg_coordoperation SET coord_op_method_code = 9665 WHERE coord_op_method_code = 1048');
$sqlite->exec('UPDATE epsg_coordoperation SET coord_op_method_code = 1088 WHERE coord_op_code IN (9465, 9466, 9467, 9693)');
$sqlite->exec('UPDATE epsg_coordoperation SET coord_op_method_code = 9665 WHERE coord_op_code IN (5656, 5657, 8451, 9461, 9692)');

$sqlite->exec("UPDATE epsg_coordoperationparamvalue SET param_value_file_ref = REPLACE(param_value_file_ref, '.gsb', '.gtx') WHERE coord_op_method_code IN (1083, 1048)");
$sqlite->exec('UPDATE epsg_coordoperationparamvalue SET coord_op_method_code = 1088 WHERE coord_op_method_code = 1083');
$sqlite->exec('UPDATE epsg_coordoperationparamvalue SET coord_op_method_code = 9665 WHERE coord_op_method_code = 1048');
$sqlite->exec("UPDATE epsg_coordoperationparamvalue SET param_value_file_ref = REPLACE(param_value_file_ref, '.gsb', '.gtx') WHERE coord_op_code IN (9465, 9466, 9467, 9693) OR coord_op_code IN (5656, 5657, 8451, 9692)");
$sqlite->exec('UPDATE epsg_coordoperationparamvalue SET coord_op_method_code = 1088 WHERE coord_op_code IN (9465, 9466, 9467, 9693)');
$sqlite->exec('UPDATE epsg_coordoperationparamvalue SET coord_op_method_code = 9665 WHERE coord_op_code IN (5656, 5657, 8451, 9692)');

/*
* VERTCON files/extents described in EPSG are the old VERTCON1/2 .94 file format, not the VERTCON3 .b files
* which cover geographically different extents. Therefore we use the official NOAA VDatum GTX conversions of the
* older files instead to avoid yet-another-grid-format implementation.
*/
$sqlite->exec('UPDATE epsg_coordoperation SET coord_op_method_code = 1084 WHERE coord_op_code IN (7969, 7970, 7971)');
$sqlite->exec("UPDATE epsg_coordoperationparamvalue SET param_value_file_ref = REPLACE(param_value_file_ref, '.94', '.gtx') WHERE coord_op_code IN (7969, 7970, 7971)");
$sqlite->exec('UPDATE epsg_coordoperationparamvalue SET coord_op_method_code = 1084 WHERE coord_op_code IN (7969, 7970, 7971)');

$sqlite->exec('VACUUM');
$sqlite->close();
Expand Down
3 changes: 3 additions & 0 deletions src/EPSG/Import/FilenameToProviderMap.php
Original file line number Diff line number Diff line change
Expand Up @@ -191,4 +191,7 @@
'Ranc08_Circe.mnt' => CoordinateOperation\IGNFHeightRGNC9193NGNC08NewCaledoniaProvider::class,
'EGM08_REDNAP.txt' => CoordinateOperation\IGNESHeightETRS89REDNAPSpainProvider::class,
'EGM08_REDNAP_Canarias.txt' => CoordinateOperation\IGNESHeightETRS89REDNAPCanaryIslandsProvider::class,
'vertconc.gtx' => CoordinateOperation\GTXNGVD29NAVD88CONUSCentralProvider::class,
'vertcone.gtx' => CoordinateOperation\GTXNGVD29NAVD88CONUSEastProvider::class,
'vertconw.gtx' => CoordinateOperation\GTXNGVD29NAVD88CONUSWestProvider::class,
];
14 changes: 14 additions & 0 deletions tests/VerticalPointTest.php
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
use PHPCoord\CoordinateOperation\CoordinateOperations;
use PHPCoord\CoordinateOperation\CRSTransformations;
use PHPCoord\CoordinateOperation\GTXDunedin1958NZVD2016Provider;
use PHPCoord\CoordinateOperation\GTXNGVD29NAVD88CONUSCentralProvider;
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
use PHPCoord\CoordinateReferenceSystem\Vertical;
Expand Down Expand Up @@ -153,6 +154,19 @@ public function testVerticalOffsetGTXReverse(): void
self::assertEqualsWithDelta(50.000, $to->getHeight()->getValue(), 0.001);
}

public function testVerticalOffsetGTXUSVERTCON(): void
{
if (!class_exists(GTXNGVD29NAVD88CONUSCentralProvider::class)) {
self::markTestSkipped('Requires phpcoord/datapack-northamerica');
}
$horizontalPoint = GeographicPoint::create(new Degree(29.4667897), new Degree(-98.4803739), null, Geographic2D::fromSRID(Geographic2D::EPSG_NAD83));
$from = VerticalPoint::create(new Metre(247.47), Vertical::fromSRID(Vertical::EPSG_NGVD29_HEIGHT_M));
$toCRS = Vertical::fromSRID(Vertical::EPSG_NAVD88_HEIGHT);
$to = $from->verticalOffsetGTX($toCRS, (new GTXNGVD29NAVD88CONUSCentralProvider())->provideGrid(), '', false, $horizontalPoint);

self::assertEqualsWithDelta(247.599, $to->getHeight()->getValue(), 0.001);
}

/**
* @group integration
* @dataProvider supportedOperations
Expand Down

0 comments on commit b7cee81

Please sign in to comment.