Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for NADCON5 grids and transformation method #3510

Merged
merged 12 commits into from
Dec 17, 2022
Merged

Add support for NADCON5 grids and transformation method #3510

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
fa97a9c
Add a +proj=gridshift (general grid shift) method to be used for NADC…
rouault Dec 14, 2022
7888f3b
GTiff grid reading: speed-up retrieving data in last accessed tile
rouault Dec 14, 2022
b15909e
gridshift: add a +no_z_transform flag
rouault Dec 14, 2022
b55fd3c
Add NADCON5 transformation method support and database entries
rouault Dec 14, 2022
526178d
metadata.sql: bump PROJ_DATA.VERSION to 1.13
rouault Dec 14, 2022
ce95d2e
Database: add concatenated operations for NADCON5 transformations
rouault Dec 14, 2022
69cfda9
Add a GenericShiftGrid::valuesAt() method for faster value retrieval
rouault Dec 14, 2022
f29145c
gridshift: use GenericShiftGrid::valuesAt()
rouault Dec 14, 2022
9c3f30f
pj_get_suggested_operation(): make sure to use NADCON5 grids rather t…
rouault Dec 14, 2022
c159d42
bench_proj_trans: add options to add random noise
rouault Dec 14, 2022
4ca9950
gridshift: cache grid values when interpolating repeatedly in the sam…
rouault Dec 14, 2022
e732237
pj_get_suggested_operation(): tune so that NADCON5 grid is used for N…
rouault Dec 15, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added data/tests/us_noaa_nadcon5_nad83_2007_nad83_2011_alaska_extract.tif
View file
Open in desktop
Binary file not shown.
Binary file added data/tests/us_noaa_nadcon5_nad83_2007_nad83_2011_conus_extract.tif
View file
Open in desktop
Binary file not shown.
70 changes: 70 additions & 0 deletions docs/source/operations/transformations/gridshift.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
.. _gridshift:

================================================================================
General grid shift
================================================================================

.. versionadded:: 9.2.0

Translation of geodetic coordinates using a grid shift.

+-----------------+-------------------------------------------------------------------+
| **Alias** | gridshift |
+-----------------+-------------------------------------------------------------------+
| **Domain** | 2D and 3D |
+-----------------+-------------------------------------------------------------------+
| **Input type** | Geodetic coordinates (horizontal), meters (vertical) |
+-----------------+-------------------------------------------------------------------+
| **Output type** | Geodetic coordinates (horizontal), meters (vertical) |
+-----------------+-------------------------------------------------------------------+

The transformation may apply horizontal geodetic offsetting and/or vertical
(ellipsoidal or orthometric height) offsetting, depending on the type of the
grid(s).

This is a generalization of the :ref:`hgridshift` and :ref:`vgridshift` methods,
that may be used in particular for US NADCON5 grids that contain both horizontal
geodetic and ellipsoidal height offsets.


Example
-------------------------------------------------------------------------------

Transformation of a point from NAD83(NSRS2007) to NAD83(2011) on conterminous USA
using a NADCON5 grid that contains horizontal and ellipsoidal height offsets.

::

+proj=gridshift +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif

Parameters
################################################################################

Required
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

.. option:: +grids=<list>

Comma-separated list of grids to load. If a grid is prefixed by an ``@`` the
grid is considered optional and PROJ will the not complain if the grid is
not available.

Grids must be in GeoTIFF format (:ref:`geodetictiffgrids`) and have an
explicit TYPE metadata item whose value is ``HORIZONTAL_OFFSET``,
``GEOGRAPHIC_3D_OFFSET``, ``ELLIPSOIDAL_HEIGHT_OFFSET``
``VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL`` or ``VERTICAL_OFFSET_VERTICAL_TO_VERTICAL``.

Optional
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

.. option:: +interpolation=bilinear/biquadratic

Default is bilinear, unless the grid contains a ``interpolation_method``
metadata item specifying the method.
Biquadratic is typically used for NADCON5 grids, and is defined in
`NOAA Technical Memorandum NOS NGS 84 - Biquadratic Interpolation
<https://geodesy.noaa.gov/library/pdfs/NOAA_TM_NOS_NGS_0084.pdf>`__

.. option:: +no_z_transform

If specified, vertical coordinate transformation will be skipped.
Comment on lines +68 to +70
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What use-case is behind this addition? A similar thing could be done with pop and push

Copy link
Member Author

@rouault rouault Dec 16, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use-case is for example NAD83 to NAD83(2011) for which elliposoidal height changes are not significant given that NAD83 is 2D-only officially and the grid for the transformation NAD83 to NAD83(HARN) has no ellipsoidal height offset. Hence it is useless to pay the price for the ellipsoidal height offset interpolation for the 3 next grids of the concatenated transformation.

Below performance with the pipeline output by "projinfo -s NAD83 -t NAD83(2011)" that will insert +no_z_transform where it makes sense

$ PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/bench_proj_trans \
    --pipeline "+proj=pipeline +step +proj=axisswap +order=2,1 \
                     +step +proj=unitconvert +xy_in=deg +xy_out=rad \
                     +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_1986_nad83_harn_conus.tif \
                     +step +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif \
                     +step +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif \
                     +step +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif \
                     +step +proj=unitconvert +xy_in=rad +xy_out=deg \
                     +step +proj=axisswap +order=2,1" 40 -100 --noise-x 1 --noise-y 1
40 -100 -> 39.9999988651031 -99.9999991253051
Duration: 3293 ms
Throughput: 1.52 million coordinates/s

compared to using push/pop and letting z interpolation:

$ PROJ_DATA=/home/even/proj/PROJ-data/us_noaa:data bin/bench_proj_trans \
       --pipeline "+proj=pipeline +step +proj=axisswap +order=2,1 \
                        +step +proj=unitconvert +xy_in=deg +xy_out=rad \
                        +step +proj=push +v_3 \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_1986_nad83_harn_conus.tif \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif \
                        +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif \
                        +step +proj=pop +v_3 \
                        +step +proj=unitconvert +xy_in=rad +xy_out=deg \
                        +step +proj=axisswap +order=2,1" 40 -100 --noise-x 1 --noise-y 1
40 -100 -> 39.9999988651031 -99.9999991253051
Duration: 3604 ms
Throughput: 1.39 million coordinates/s

Juts seeing that the difference is no longer that big since I added various optimizations in grid caching & interpolation, but it was more significant at an earlier stage of developpement. For some grids like Alaska where there's a dedicated subgrid for the ellipsoidal height offset, skipping the Z interpolation saves I/O, which can speed up things in a remote grid use case.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use-case is for example NAD83 to NAD83(2011) for which elliposoidal height changes are not significant given that NAD83 is 2D-only officially and the grid for the transformation NAD83 to NAD83(HARN) has no ellipsoidal height offset. Hence it is useless to pay the price for the ellipsoidal height offset interpolation for the 3 next grids of the concatenated transformation.

I would expect that NAD83(1986) does not have ellipsoidal or XYZ forms. I guess you are saying that once you have started with a 2D datum the point is to skip additional work and that the user has to treat the result as 2D only even though it is a value in a 3D CRS. Did I follow that right?

Copy link
Member Author

@rouault rouault Dec 16, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did I follow that right?

yes, taking into account the ellipsoidal height offsets of the HARN -> FBN -> NSRS2007 -> 2011 transformations when the starting point is 86 would give a false sense of accuracy given that the offset between 86 and HARN is probably much more bigger but nobody knows its value.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the explanation. It makes a lot of sense to optimize performance!

2 changes: 2 additions & 0 deletions docs/source/operations/transformations/hgridshift.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@ about all three formats can be found in the GDAL documentation and/or driver sou
code. GDAL reads and writes all three formats. Using GDAL for construction of
new grids is recommended.

To apply as well ellipsoidal height differences sometimes present in some grids
(such as US NADCON5 grids), use the :ref:`gridshift` method.

Temporal gridshifting
################################################################################
Expand Down
1 change: 1 addition & 0 deletions docs/source/operations/transformations/index.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ systems are based on different datums.
defmodel
deformation
geogoffset
gridshift
helmert
horner
molodensky
Expand Down
32 changes: 30 additions & 2 deletions docs/source/specifications/geodetictiffgrids.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,21 @@ is an easy way to inspect such grid files:
as a speed per year for temporal gridshifting.
Corresponds to PROJ :ref:`hgridshift` method.

- ``GEOGRAPHIC_3D_OFFSET``: implies the presence of at least 3 samples.
The first sample must contain the latitude offset, the second
sample must contain the longitude offset and the third one the ellipsoidal
height difference. Typically used for NADCON5 grids.
Added in PROJ 9.2
Corresponds to PROJ :ref:`gridshift` method.

- ``ELLIPSOIDAL_HEIGHT_OFFSET``: implies the presence of one sample with
the ellipsoidal height difference. Generally used in combination with
another grid of type ``HORIZONTAL_OFFSET`` to perform Geographic 3D
offseting when the horizontal and vertical grids do not have the same
resolution, as found in some NADCON5 grids.
Added in PROJ 9.2
Corresponds to PROJ :ref:`gridshift` method.

- ``VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL``: implies the presence of at least one sample.
The first sample must contain the vertical adjustment. Must be used when
the source/interpolation CRS is a Geographic CRS and the target CRS a Vertical CRS.
Expand Down Expand Up @@ -253,14 +268,21 @@ is an easy way to inspect such grid files:

Values recognized by PROJ for this Item are currently:

+ ``latitude_offset``: valid for TYPE=HORIZONTAL_OFFSET. Sample values should be
+ ``latitude_offset``: valid for TYPE=HORIZONTAL_OFFSET or
GEOGRAPHIC_3D_OFFSET . Sample values should be
the value to add a latitude expressed in the CRS encoded in the GeoKeys
to obtain a latitude value expressed in the target CRS.

+ ``longitude_offset``: valid for TYPE=HORIZONTAL_OFFSET. Sample values should be
+ ``longitude_offset``: valid for TYPE=HORIZONTAL_OFFSET or
GEOGRAPHIC_3D_OFFSET . Sample values should be
the value to add a longitude expressed in the CRS encoded in the GeoKeys
to obtain a longitude value expressed in the target CRS.

+ ``ellipsoidal_height_offset``: valid for TYPE=ELLIPSOIDAL_HEIGHT_OFFSET or
GEOGRAPHIC_3D_OFFSET . Sample values should be the value to add to the
ellipsoidal height of the source CRS to obtain the ellipsoidal height of
the target CRS.

+ ``geoid_undulation``: valid for TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL.
For a source CRS being a geographic CRS and a target CRS being a vertical CRS,
sample values should be the value to add to a geoid-related height (that
Expand Down Expand Up @@ -415,6 +437,10 @@ is an easy way to inspect such grid files:
of this grid).
Will be ignored by PROJ (this information can be inferred by the grids extent)

* The ``interpolation_method`` metadata item may be present to indicate the
interpolation method to apply. ``bilinear`` or ``biquadratic`` are supported.
If not specified, defaults to ``bilinear``.

Example
+++++++

Expand Down Expand Up @@ -648,5 +674,7 @@ georeferencing is contained within the first 40 KB of the file.
Revisions
+++++++++

* v0.3: addition of TYPE=GEOGRAPHIC_3D_OFFSET, ELLIPSOIDAL_HEIGHT_OFFSET and
interpolation_method (PROJ 9.2)
* v0.2: addition of "arc-seconds per year" as a valid unit (PROJ 9.1.1)
* v0.1: initial version for PROJ 7.0
1 change: 1 addition & 0 deletions scripts/reference_exported_symbols.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,7 @@ osgeo::proj::GenericShiftGrid::gridAt(double, double) const
osgeo::proj::GenericShiftGridSet::~GenericShiftGridSet()
osgeo::proj::GenericShiftGridSet::GenericShiftGridSet()
osgeo::proj::GenericShiftGridSet::gridAt(double, double) const
osgeo::proj::GenericShiftGridSet::gridAt(std::string const&, double, double) const
osgeo::proj::GenericShiftGridSet::open(pj_ctx*, std::string const&)
osgeo::proj::GenericShiftGridSet::reassign_context(pj_ctx*)
osgeo::proj::GenericShiftGridSet::reopen(pj_ctx*)
Expand Down
1 change: 1 addition & 0 deletions scripts/reformat_cpp.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp \
"$TOPDIR"/src/transformations/tinshift_exceptions.hpp \
"$TOPDIR"/src/transformations/tinshift_impl.hpp \
"$TOPDIR"/src/transformations/tinshift.cpp \
"$TOPDIR"/src/transformations/gridshift.cpp \
"$TOPDIR"/src/quadtree.hpp \
; do
if ! echo "$i" | grep -q "lru_cache.hpp"; then
Expand Down
Loading