From d747be2506bab01f37fc0229b83e4beb4aa4b0e2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 9 Nov 2025 21:44:13 +0100 Subject: [PATCH 1/2] Doc: move the description of the TIN JSON format into its own page --- .../operations/transformations/tinshift.rst | 101 +---------------- docs/source/specifications/index.rst | 2 + docs/source/specifications/tin_json.rst | 104 ++++++++++++++++++ 3 files changed, 110 insertions(+), 97 deletions(-) create mode 100644 docs/source/specifications/tin_json.rst diff --git a/docs/source/operations/transformations/tinshift.rst b/docs/source/operations/transformations/tinshift.rst index 30d191d8ef..07700dc740 100644 --- a/docs/source/operations/transformations/tinshift.rst +++ b/docs/source/operations/transformations/tinshift.rst @@ -116,100 +116,7 @@ the inverse transformation), otherwise which triangle will be selected is unspecified. Besides that, the triangulation does not need to have particular properties (like being a Delaunay triangulation) -File format -+++++++++++ - -The triangulation is stored in a text-based format, using JSON as a serialization. - -Below a minimal example, from the KKJ to ETRS89 transformation, with just a -single triangle: - -.. literalinclude:: ../../../../data/tests/tinshift_crs_implicit.json - :language: json - - -So after the generic metadata, we define the input and output CRS (informative -only), and that the transformation affects horizontal components of -coordinates. We name the columns of the ``vertices`` and ``triangles`` arrays. -We defined the source and target coordinates of each vertex, and define a -triangle by referring to the index of its vertices in the ``vertices`` array. - -More formally, the specific items for the triangulation file are: - -input_crs - String identifying the CRS of source coordinates - in the vertices. Typically ``EPSG:XXXX``. If the transformation is for vertical - component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY - where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). - For example, for the KKJ->ETRS89 transformation, this is EPSG:2393 - (``KKJ / Finland Uniform Coordinate System``). The input coordinates are assumed - to be passed in the "normalized for visualisation" / "GIS friendly" order, - that is longitude, latitude for geographic coordinates and - easting, northing for projected coordinates. - - -output_crs - String identifying the CRS of target coordinates in the vertices. - Typically ``EPSG:XXXX``. If the transformation is for vertical component, - this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where - XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). - For example, for the KKJ->ETRS89 transformation, this is EPSG:3067 - (\"ETRS89 / TM35FIN(E,N)\"). The output coordinates will be returned in - the "normalized for visualisation" / "GIS friendly" order, - that is longitude, latitude for geographic coordinates and - easting, northing for projected coordinates. - - -transformed_components - Array which may contain one or two strings: "horizontal" when horizontal - components of the coordinates are transformed and/or "vertical" when the - vertical component is transformed. - - -fallback_strategy - String identifying how to treat points that do not fall into any of the - specified triangles. This item is available for ``format_version`` >= 1.1. - Possible values are ``none``, ``nearest_side`` and ``nearest_centroid``. The - default is ``none`` and signifies, that points that fall outside the - specified triangles are not transformed. This is also the behavior for - ``format_version`` before 1.1. If ``fallback_strategy`` is set to - ``nearest_side``, then points that do not fall into any triangle are - transformed according to the triangle closest to them by euclidean distance. - If ``fallback_strategy`` is set to ``nearest_centroid``, then points that do - not fall into any triangle are transformed according to the triangle with the - closest centroid (intersection of its medians). - - -vertices_columns - Specify the name of the columns of the rows in the ``vertices`` - array. There must be exactly as many elements in ``vertices_columns`` as in a - row of ``vertices``. The following names have a special meaning: ``source_x``, - ``source_y``, ``target_x``, ``target_y``, ``source_z``, ``target_z`` and - ``offset_z``. ``source_x`` and ``source_y`` are compulsory. - ``source_x`` is for the source longitude (in degree) or easting. - ``source_y`` is for the source latitude (in degree) or northing. - ``target_x`` and ``target_y`` are compulsory when ``horizontal`` is specified - in ``transformed_components``. (``source_z`` and ``target_z``) or - ``offset_z`` are compulsory when ``vertical`` is specified in ``transformed_components`` - - -triangles_columns - Specify the name of the columns of the rows in the - ``triangles`` array. There must be exactly as many elements in ``triangles_columns`` - as in a row of ``triangles``. The following names have a special meaning: - ``idx_vertex1``, ``idx_vertex2``, ``idx_vertex3``. They are compulsory. - - -vertices - An array whose items are themselves arrays with as many columns as - described in ``vertices_columns``. - - -triangles - An array whose items are themselves arrays with as many columns as - described in ``triangles_columns``. - The value of the ``idx_vertexN`` columns must be indices - (between 0 and len(``vertices``-1) of items of the ``vertices`` array. - -A `JSON schema `_ is available -for this file format. +JSON File format +++++++++++++++++ + +See :ref:`tin_json`. diff --git a/docs/source/specifications/index.rst b/docs/source/specifications/index.rst index 1853ca53a5..2dfff9f7b1 100644 --- a/docs/source/specifications/index.rst +++ b/docs/source/specifications/index.rst @@ -12,3 +12,5 @@ for the sake of broader interoperability. projjson geodetictiffgrids + tin_json + diff --git a/docs/source/specifications/tin_json.rst b/docs/source/specifications/tin_json.rst new file mode 100644 index 0000000000..e898126d39 --- /dev/null +++ b/docs/source/specifications/tin_json.rst @@ -0,0 +1,104 @@ +.. _tin_json: + +================================================================================ +Triangulated irregular network (TIN) JSON format +================================================================================ + +The TIN JSON format defines the format for a triangulation model stored in a JSON +file. + +That file format may be used by the :ref:`tinshift ` operation since +PROJ 7.2 + +Below a minimal example, from the KKJ to ETRS89 transformation, with just a +single triangle: + +.. literalinclude:: ../../../data/tests/tinshift_crs_implicit.json + :language: json + + +So after the generic metadata, we define the input and output CRS (informative +only), and that the transformation affects horizontal components of +coordinates. We name the columns of the ``vertices`` and ``triangles`` arrays. +We defined the source and target coordinates of each vertex, and define a +triangle by referring to the index of its vertices in the ``vertices`` array. + +More formally, the specific items for the triangulation file are: + +input_crs + String identifying the CRS of source coordinates + in the vertices. Typically ``EPSG:XXXX``. If the transformation is for vertical + component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY + where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). + For example, for the KKJ->ETRS89 transformation, this is EPSG:2393 + (``KKJ / Finland Uniform Coordinate System``). The input coordinates are assumed + to be passed in the "normalized for visualisation" / "GIS friendly" order, + that is longitude, latitude for geographic coordinates and + easting, northing for projected coordinates. + + +output_crs + String identifying the CRS of target coordinates in the vertices. + Typically ``EPSG:XXXX``. If the transformation is for vertical component, + this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where + XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). + For example, for the KKJ->ETRS89 transformation, this is EPSG:3067 + (\"ETRS89 / TM35FIN(E,N)\"). The output coordinates will be returned in + the "normalized for visualisation" / "GIS friendly" order, + that is longitude, latitude for geographic coordinates and + easting, northing for projected coordinates. + + +transformed_components + Array which may contain one or two strings: "horizontal" when horizontal + components of the coordinates are transformed and/or "vertical" when the + vertical component is transformed. + + +fallback_strategy + String identifying how to treat points that do not fall into any of the + specified triangles. This item is available for ``format_version`` >= 1.1. + Possible values are ``none``, ``nearest_side`` and ``nearest_centroid``. The + default is ``none`` and signifies, that points that fall outside the + specified triangles are not transformed. This is also the behavior for + ``format_version`` before 1.1. If ``fallback_strategy`` is set to + ``nearest_side``, then points that do not fall into any triangle are + transformed according to the triangle closest to them by euclidean distance. + If ``fallback_strategy`` is set to ``nearest_centroid``, then points that do + not fall into any triangle are transformed according to the triangle with the + closest centroid (intersection of its medians). + + +vertices_columns + Specify the name of the columns of the rows in the ``vertices`` + array. There must be exactly as many elements in ``vertices_columns`` as in a + row of ``vertices``. The following names have a special meaning: ``source_x``, + ``source_y``, ``target_x``, ``target_y``, ``source_z``, ``target_z`` and + ``offset_z``. ``source_x`` and ``source_y`` are compulsory. + ``source_x`` is for the source longitude (in degree) or easting. + ``source_y`` is for the source latitude (in degree) or northing. + ``target_x`` and ``target_y`` are compulsory when ``horizontal`` is specified + in ``transformed_components``. (``source_z`` and ``target_z``) or + ``offset_z`` are compulsory when ``vertical`` is specified in ``transformed_components`` + + +triangles_columns + Specify the name of the columns of the rows in the + ``triangles`` array. There must be exactly as many elements in ``triangles_columns`` + as in a row of ``triangles``. The following names have a special meaning: + ``idx_vertex1``, ``idx_vertex2``, ``idx_vertex3``. They are compulsory. + + +vertices + An array whose items are themselves arrays with as many columns as + described in ``vertices_columns``. + + +triangles + An array whose items are themselves arrays with as many columns as + described in ``triangles_columns``. + The value of the ``idx_vertexN`` columns must be indices + (between 0 and len(``vertices``-1) of items of the ``vertices`` array. + +A `JSON schema `_ is available +for this file format. From 863f356ae3ba546e9b4a1a5a7539abf220e3ce3b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 14 Nov 2025 15:37:03 +0100 Subject: [PATCH 2/2] tinshift: support TIN GeoPackage files Such files are equivalent to already supported TIN JSON files, but scale better for arbitrarily large triangulations, in particular for network-based access. The format is defined in source/specifications/tin_gpkg.rst and the https://github.com/OSGeo/PROJ-data/blob/master/grid_tools/tin_json_to_gpkg_tin.py Python script may be used o convert an existing TIN JSON into a TIN GeoPackage. --- .github/workflows/fedora_rawhide/start.sh | 4 + .github/workflows/windows.yml | 2 +- data/CMakeLists.txt | 2 +- data/sql/grid_alternatives.sql | 3 + data/sql/proj_db_table_defs.sql | 6 +- data/tests/tinshift_crs_implicit.gpkg | Bin 0 -> 22528 bytes data/tests/tinshift_empty_file.gpkg | 0 .../tinshift_fallback_nearest_centroid.gpkg | Bin 0 -> 21504 bytes .../tests/tinshift_fallback_nearest_side.gpkg | Bin 0 -> 21504 bytes data/tests/tinshift_simplified_kkj_etrs.gpkg | Bin 0 -> 24576 bytes data/tests/tinshift_simplified_n60_n2000.gpkg | Bin 0 -> 24576 bytes docs/source/install.rst | 7 +- .../operations/transformations/tinshift.rst | 18 +- docs/source/specifications/index.rst | 1 + docs/source/specifications/tin_gpkg.rst | 241 +++++ docs/source/specifications/tin_json.rst | 3 + scripts/reference_exported_symbols.txt | 2 +- src/filemanager.cpp | 8 +- src/grids.cpp | 6 +- src/init.cpp | 6 +- src/iso19111/factory.cpp | 6 +- src/iso19111/io.cpp | 3 +- src/lib_proj.cmake | 11 +- src/proj_internal.h | 2 +- src/sqlite3_utils.cpp | 196 ++++ src/sqlite3_utils.hpp | 2 + src/transformations/tinshift.cpp | 70 +- src/transformations/tinshift_gpkg.cpp | 945 ++++++++++++++++++ src/transformations/tinshift_gpkg.hpp | 135 +++ src/transformations/tinshift_iface.hpp | 33 + .../{tinshift.hpp => tinshift_json.hpp} | 37 +- ...tions.hpp => tinshift_json_exceptions.hpp} | 6 +- ...nshift_impl.hpp => tinshift_json_impl.hpp} | 30 +- test/CMakeLists.txt | 5 + test/gie/tinshift_gpkg.gie | 58 ++ test/gie/tinshift_gpkg_network.gie | 27 + test/unit/test_network.cpp | 6 +- test/unit/test_tinshift.cpp | 46 +- 38 files changed, 1833 insertions(+), 94 deletions(-) create mode 100644 data/tests/tinshift_crs_implicit.gpkg create mode 100644 data/tests/tinshift_empty_file.gpkg create mode 100644 data/tests/tinshift_fallback_nearest_centroid.gpkg create mode 100644 data/tests/tinshift_fallback_nearest_side.gpkg create mode 100644 data/tests/tinshift_simplified_kkj_etrs.gpkg create mode 100644 data/tests/tinshift_simplified_n60_n2000.gpkg create mode 100644 docs/source/specifications/tin_gpkg.rst create mode 100644 src/transformations/tinshift_gpkg.cpp create mode 100644 src/transformations/tinshift_gpkg.hpp create mode 100644 src/transformations/tinshift_iface.hpp rename src/transformations/{tinshift.hpp => tinshift_json.hpp} (91%) rename src/transformations/{tinshift_exceptions.hpp => tinshift_json_exceptions.hpp} (95%) rename src/transformations/{tinshift_impl.hpp => tinshift_json_impl.hpp} (96%) create mode 100644 test/gie/tinshift_gpkg.gie create mode 100644 test/gie/tinshift_gpkg_network.gie diff --git a/.github/workflows/fedora_rawhide/start.sh b/.github/workflows/fedora_rawhide/start.sh index dfc0a220b3..f776ff2236 100755 --- a/.github/workflows/fedora_rawhide/start.sh +++ b/.github/workflows/fedora_rawhide/start.sh @@ -26,8 +26,10 @@ ctest -j$(nproc) # Try EMBED_RESOURCE_DIRECTORY option wget https://raw.githubusercontent.com/OSGeo/PROJ-data/refs/heads/master/us_nga/us_nga_egm96_15.tif +wget https://cdn.proj.org/fi_nls_ykj_etrs35fin.gpkg mkdir grids mv us_nga_egm96_15.tif grids +mv fi_nls_ykj_etrs35fin.gpkg grids echo "Build with -DEMBED_RESOURCE_FILES=ON -DEMBED_RESOURCE_DIRECTORY=$PWD/grids" CC=clang CXX=clang++ cmake .. -DEMBED_RESOURCE_DIRECTORY=$PWD/grids make -j$(nproc) @@ -36,6 +38,8 @@ echo 49 2 0 | bin/cs2cs "WGS84 + EGM96 height" EPSG:4979 echo 49 2 0 | bin/cs2cs "WGS84 + EGM96 height" EPSG:4979 | grep 44.643 >/dev/null || (echo "Expected 49dN 2dE 44.643 as a result" && /bin/false) echo 0 0 0 | bin/cct +init=ITRF2000:ITRF96 echo 0 0 0 | bin/cct +init=ITRF2000:ITRF96 | grep 0.0067 >/dev/null || (echo "Expected 0.0067 0.0061 -0.0185 as a result" && /bin/false) +echo 3432087 6995748 0 | bin/cct +proj=tinshift +file=fi_nls_ykj_etrs35fin.gpkg +echo 3432087 6995748 0 | bin/cct +proj=tinshift +file=fi_nls_ykj_etrs35fin.gpkg | grep 431943.0905 >/dev/null || (echo "Expected 431943.0905 6992816.7826 0 as a result" && /bin/false) echo "Build with -DEMBED_RESOURCE_FILES=ON -DEMBED_RESOURCE_DIRECTORY=$PWD/grids -DUSE_ONLY_EMBEDDED_RESOURCE_FILES=ON" CC=clang CXX=clang++ cmake .. -DEMBED_RESOURCE_DIRECTORY=$PWD/grids -DUSE_ONLY_EMBEDDED_RESOURCE_FILES=ON diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 9f86d768bd..a5ca709c70 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -55,7 +55,7 @@ jobs: shell: cmd if: steps.cache.outputs.cache-hit != 'true' run: | - vcpkg install sqlite3[core,tool] tiff curl --triplet=${{ env.ARCH }}-windows + vcpkg install sqlite3[core,tool,rtree] tiff curl --triplet=${{ env.ARCH }}-windows - name: Build shell: cmd diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt index 8df4644152..487c85b24a 100644 --- a/data/CMakeLists.txt +++ b/data/CMakeLists.txt @@ -31,7 +31,7 @@ set(ALL_SQL_IN "${CMAKE_CURRENT_BINARY_DIR}/all.sql.in") set(PROJ_DB "${CMAKE_CURRENT_BINARY_DIR}/proj.db") include(sql_filelist.cmake) -set(PROJ_DB_SQL_EXPECTED_MD5 "cca51e9e250f9c979e8f18bb483433e5") +set(PROJ_DB_SQL_EXPECTED_MD5 "a71b63a199da368414be93d819bf2e5b") add_custom_command( OUTPUT ${PROJ_DB} diff --git a/data/sql/grid_alternatives.sql b/data/sql/grid_alternatives.sql index e82c1bd019..89e465996c 100644 --- a/data/sql/grid_alternatives.sql +++ b/data/sql/grid_alternatives.sql @@ -183,6 +183,9 @@ VALUES ('fi_nls_n43_n60.json','fi_nls_n43_n60.json',NULL,'JSON','tinshift',0,NULL,'https://cdn.proj.org/fi_nls_n43_n60.json',1,1,NULL), ('fi_nls_n60_n2000.json','fi_nls_n60_n2000.json',NULL,'JSON','tinshift',0,NULL,'https://cdn.proj.org/fi_nls_n60_n2000.json',1,1,NULL), ('fi_nls_ykj_etrs35fin.json','fi_nls_ykj_etrs35fin.json',NULL,'JSON','tinshift',0,NULL,'https://cdn.proj.org/fi_nls_ykj_etrs35fin.json',1,1,NULL), +('NOT-YET-IN-GRID-TRANSFORMATION-fi_nls_n43_n60.gpkg','fi_nls_n43_n60.gpkg',NULL,'GPKG','tinshift',0,NULL,'https://cdn.proj.org/fi_nls_n43_n60.gpkg',1,1,NULL), +('NOT-YET-IN-GRID-TRANSFORMATION-fi_nls_n60_n2000.gpkg','fi_nls_n60_n2000.gpkg',NULL,'GPKG','tinshift',0,NULL,'https://cdn.proj.org/fi_nls_n60_n2000.gpkg',1,1,NULL), +('NOT-YET-IN-GRID-TRANSFORMATION-fi_nls_ykj_etrs35fin.gpkg','fi_nls_ykj_etrs35fin.gpkg',NULL,'GPKG','tinshift',0,NULL,'https://cdn.proj.org/fi_nls_ykj_etrs35fin.gpkg',1,1,NULL), -- fr_ign - IGN France ('rgf93_ntf.gsb','fr_ign_ntf_r93.tif','ntf_r93.gsb','GTiff','hgridshift',1,NULL,'https://cdn.proj.org/fr_ign_ntf_r93.tif',1,1,NULL), diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql index 073245a315..0b9f17b9bd 100644 --- a/data/sql/proj_db_table_defs.sql +++ b/data/sql/proj_db_table_defs.sql @@ -896,7 +896,7 @@ CREATE TABLE grid_alternatives( original_grid_name TEXT NOT NULL PRIMARY KEY, -- original grid name (e.g. Und_min2.5x2.5_egm2008_isw=82_WGS84_TideFree.gz). For LOS/LAS format, the .las files proj_grid_name TEXT NOT NULL, -- PROJ >= 7 grid name (e.g us_nga_egm08_25.tif) old_proj_grid_name TEXT, -- PROJ < 7 grid name (e.g egm08_25.gtx) - proj_grid_format TEXT NOT NULL, -- 'GTiff', 'GTX', 'NTv2', JSON + proj_grid_format TEXT NOT NULL, -- 'GTiff', 'GTX', 'NTv2', 'JSON', 'GPKG' proj_method TEXT NOT NULL, -- gridshift, hgridshift, vgridshift, geoid_like, geocentricoffset, tinshift or velocity_grid inverse_direction BOOLEAN NOT NULL CHECK (inverse_direction IN (0, 1)), -- whether the PROJ grid direction is reversed w.r.t to the authority one (TRUE in that case) package_name TEXT, -- no longer used. Must be NULL @@ -906,14 +906,14 @@ CREATE TABLE grid_alternatives( directory TEXT, -- optional directory where the file might be located CONSTRAINT fk_grid_alternatives_grid_packages FOREIGN KEY (package_name) REFERENCES grid_packages(package_name) ON DELETE CASCADE, - CONSTRAINT check_grid_alternatives_grid_fromat CHECK (proj_grid_format IN ('GTiff', 'GTX', 'NTv2', 'JSON')), + CONSTRAINT check_grid_alternatives_grid_fromat CHECK (proj_grid_format IN ('GTiff', 'GTX', 'NTv2', 'JSON', 'GPKG')), CONSTRAINT check_grid_alternatives_proj_method CHECK (proj_method IN ('gridshift', 'hgridshift', 'vgridshift', 'geoid_like', 'geocentricoffset', 'tinshift', 'velocity_grid', 'defmodel')), CONSTRAINT check_grid_alternatives_inverse_direction CHECK (NOT(proj_method = 'geoid_like' AND inverse_direction = 1)), CONSTRAINT check_grid_alternatives_package_name CHECK (package_name IS NULL), CONSTRAINT check_grid_alternatives_direct_download_url CHECK (NOT(direct_download IS NULL AND url IS NOT NULL)), CONSTRAINT check_grid_alternatives_open_license_url CHECK (NOT(open_license IS NULL AND url IS NOT NULL)), CONSTRAINT check_grid_alternatives_constraint_cdn CHECK (NOT(url LIKE 'https://cdn.proj.org/%' AND (direct_download = 0 OR open_license = 0 OR url != 'https://cdn.proj.org/' || proj_grid_name))), - CONSTRAINT check_grid_alternatives_tinshift CHECK ((proj_grid_format != 'JSON' AND proj_method != 'tinshift') OR (proj_grid_format = 'JSON' AND proj_method = 'tinshift')) + CONSTRAINT check_grid_alternatives_tinshift CHECK ((proj_grid_format NOT IN ('JSON','GPKG') AND proj_method != 'tinshift') OR (proj_grid_format IN ('JSON', 'GPKG') AND proj_method = 'tinshift')) ) WITHOUT ROWID; CREATE INDEX idx_grid_alternatives_proj_grid_name ON grid_alternatives(proj_grid_name); diff --git a/data/tests/tinshift_crs_implicit.gpkg b/data/tests/tinshift_crs_implicit.gpkg new file mode 100644 index 0000000000000000000000000000000000000000..52859958875be1e912ff89a80ed70097dd7c6fd0 GIT binary patch literal 22528 zcmeHP-EZ606(>a}l9RYiySjepJl}E%WG%F<4?D7)Wwm0{vRYYkCAp2085}LrHfxDe zNGh?@7MPPRMS*=_dm6BZZV&q}1`OD+x9)9^D==U~p062lA(qunx+;Bp(ttuelNn$abxf>rvzrMT>| zT*_qUQ<*H6&dpD!=d$T^X6{mk%Wu?IE7f9s6UqUSS&q#nGf8#}%5f~@QYL)~L)NRs zH68#eXtJ!{RvLROm+ZqnpG&5Y-NH?hcYz%^r<1Ad>|7?DOV6cpx!Fv=Jmga8w3!DJ z;F4%4MpuH!7M3>drM9|j({1FB<#^W0ans?PY1*QQhrUh}l}K_!X^Bm^V=i^lR*A7E z6B~QIvnOk9Mfx@sy+Vz~CR$#zfoj>+#Lk}5!1h(HlzF#?-3P2YhsjI#Sgqgq)N#?v z)F?f%@1?RNYKE*UqUt2_sVm%lARf~RiD3U{m??_+6{Eu~bRhyE0)r4(X3kEGPmEll z!7+vg{cpD^>pQY&bTwIb3F%~dZYq_UO6KaR+(IU|Fq@gqCg-n5QcN0Vn?mGl{N?z_ z6-*`OMiO9J5)JTfccE(s$b;lvrO%F!kK5U~C2KGnc)4->gCg6RIRlwsm_Ja=KbSwk zEp#CQAp%c00%w`=VW;~)$HeKtU=kT&COirg_Wvimvcv3$2=pNEIQ<{j{}JZ56!T}e zgf2uNMBqt7;A@fS#OO;*%+n=cf1$NrDQy-iWzaWXqNC%ZV-W7F39R5^xekFi|3{cV zQq14r61ot95P_!#f%8M-qpl)=o&NcuiBWIukCcej{~w_K|KO=9cvz4Sfe-;70)CAk zp#MikBJaTY-^WzsU)0BuUq{}fenJ)CZ-6VFK0UUTr6~nYtq$~?O)wo0#IDg+2JnBu zpZFFua| z>ZDwO-$toK@Z2k3-6XH^o17#XqF@|$dU!ZRg2a>^Mb@myaQOB52G2>d-q4hep}@hC zDRh}%$#0bEBp&CQaDFH>_F!{aCUB0&*Wq-{OVnYUggxRNLutvA@e7+%7g|#nr22)` zg$rv77j~}4xi}nE#uEvyrKrLIS*~nA!8n4xL&w(=5Bx#kbAxnE7Zi!W5|JQ6mDVx|7;8;_D15Jl7 z0j&cMp{d-^>9MtVpH_vfObFD5kPdiH<(Nm0nhHc6y35efvAH?=Ucze}JvlWa^w1*5 zn@SId1;p|8!Ua+Fs?^_(+S;zDDt7=dr)4NY1Jd9JLZm+$2E@t+hdvAIWwPZuY#%^Ncc+MUN#R(#lX9)4y zlE2dsI8v?@%Mf>sFV{KXYe|CMX!mv^;V#>_#hX^FGT6_n#K$Im2+(E+2M*G_RwB^%ly)K9b9Qgyq%_L5@*2htNnS?#x1 zRP;USRCFwxrQe_MDyf&rtgK$*>Eo5w%QATIbo7OpnU2{-&4*+$1Y!{?JO#K+W)Tq(74*(~l!@lsOw6 z8y~0d=FIY8a%hj|+@B4W$|3)U1wyrlCF`I)`;;YEEr{B#YzPM|?7`$qegy20!=8|Y zwf+bDztHy2*H**F5P>Hbfp5^Umm&98hA9{(I$~pA+?D&A0cWpc|7RG9Vm@Mi3b)XO z2!sfH&IqIK=^w6IRmAKs_&1zXoGT2IBrMd~ATzh5k#Gw*47P@I?-Cx)!PIMkpx zlojW@t6n;*-#6o0^ZzFl^NE#{@ERfzBJlJha2m#n*f1Tj`u{_UdHD2~IV@p_K!|{g z07EgW36_Gx8u(frl^D-TT)0=;*{m?0)5zBHwm*h1>^w59N&_|0Z#XBwr&n zzQixpJsR#6O~XDBhQXWh!b)+O11kh~7s&5P`#l+sEWvC0lh@+2*|{XRRBmPdQZj{4 z@FKkB(1&i;y2Z&m?-0Ga4Hh|*w^Gx%H^Y(CY{E}Y_a|rill?-V1XoLys|hR$iY6=7 z%383bq{Oe($*Zs%M4Wvkaw|nDWq*;$6`+}z^y5wkxzos<9>|>ua%Ye`gWN|iJR2Qb zTB46;%uO)OfVbxbx2CMqdTUGagrFCKTV*G3vXd%&1wn21+YRKbbdVu&6niE*_R1^t zk!2&%S0XYF4EV|;WC;bFm?1m?OdOpXiH_y-^wHM`BLo}SiqtnRX2^4okQtMaoiuNK z(XoKD_W!6H7-2nT$A5pB9-@APUlMWEnMLq`f=( zK#0JX6amw;2o@)%Uq?$KFc^|m-DFwP;EN~9EtxE}TdlUL zlf3=i60}*KolcS>nAmIZB`ffYO^oa_cSq}_IvmH_D_GmI37OzUg1xq;7>B^}wrk0h zcc{d)Z34co)bug?tv|WBOI-?e5^C#cDrj#?-?z#NSbU4G>Mon8oOdl)2P^Cg6~v;!4~EP z&ZY4eJn1 z(Pr#A)@Qa4XP_QueA#fS|1(Sac#VTNa{K1dkG^C7FbxLVzI-`k69u~p`~cb9Zbgt{ ztwK_{RBDb~TPqQ;Zh%$0L>w{j5w2Y@D2ETVncw>YEzL`zUt8MoGU7 r0CxkdTHnZ*?Dy6BF_Mj1v0NZ#5SWaauYMiCRXo6}q8^&hhlzg!1LriI literal 0 HcmV?d00001 diff --git a/data/tests/tinshift_empty_file.gpkg b/data/tests/tinshift_empty_file.gpkg new file mode 100644 index 0000000000..e69de29bb2 diff --git a/data/tests/tinshift_fallback_nearest_centroid.gpkg b/data/tests/tinshift_fallback_nearest_centroid.gpkg new file mode 100644 index 0000000000000000000000000000000000000000..7f2069cdb3d3591d1fe154325c9fb2803658af80 GIT binary patch literal 21504 zcmeHP&2QVt73YU+$#%SPwlQ?bs-9(FNNVU0N&b@4G*N6?QC60`lH9~e4T6y8zkk} zCOPlDd7Syp%$qkKbou?fq=~5BR$79FQbd3tN#Ytp1VN0$-z5B*Aq)=zGXi&#AYwno z@$A-nS$y{h9h{(QlH6ceE>o_U z#6nxq_H=5YSgPFE-CgUdTDyh1R(7iGwo;R10j#ki)t0GblRZVOipZ{5X6HbBEpOf6O2yo|EwYqa zUP+{u*<^ZUHo3T*Or{o>Qfy|sa-&$vRd%2ppqXKr<@wY+vkT=I;4Gz*OPI4&%58E0 zP*xE|`H@uJW7v5Q@s;#^0*f2C)x-vfgXHXdVtHXPl}slW6Y2Cq$}11)L^7%8K?gW5 zXp+{gL1F{T9`5k!R3RIE6-hFxiB+j3G~teOsiL-wv^^2^ zuva^KqSBUX-y(unh{@1Y%WgJMEe%EJ>`7JJz6!+xXVtLvfa}f_vfn#U>vumlL+~;& zNlxwCZ&?!*O;ja8Hed3&rQCWzw$Ldr!Tq16Qw04H{Sa<`@DcD4I30luJvJ4Y8ov?> zkt4LGKrd`GMYS#pT2~QOi;$d8F3uzpGxO<6BE6bgT1_Qq7m}$vH5hQR)L8W0==c@< zqL7t>hmlFpAbPJ0z1u}TN!e9$EE0+I^72qrU_h|5V#ZHOjB_#yIiTs^67=8b-@(lf zJ_0@huQ>u^bY#Tr`{(E=>9Qo$I6Y-kXutoz=9TSd-$$SifiKekvHwrezar>Az{L+f z0zLw-8vV-8C?>*9ldrh&?*>CXGwycS=rn(~ZT8Yb4I8@18&vVGdLNPrpFAFUZRk+&~R4BsV zc0P~rzE`HSgWltI*qWdTymr{>lVPPAyiBS~qGBY5&9C>jIkqOMRYmG(5^OB#O4qrK z%yzzlqEWU9`-gmW50;fh1p9bg1$Ni$mzshxSRk%zQcH|Qr*~$iTQk$O%JhxZ>CM&Y z`kg2ng{{hHJkGWxnLj}5#ce1Ui*W4FjJ1RVCktX0OI1`}sv&R@xhz*=*NTO5rIdjv zs_*lSDzET8iH%~3%Vi4~0u(c7aa7_qxDr=bzGXwy5;Xxgb{!itGByZ`$wIMFgzUED>@*0rl#96n zq}}8S6&B=LHD0Z@`#TXb+k6ziSQWLXDXR%`O;BolSCRS(ij_g<@aU4ZHcdRx<_=m- zLtU?2ZaubebTpm|o?!j#c8Yx}QBw|~ewy8utoGV#KRI-wAw5Nv(SG|?1)mW^!SM1j z`N5Q3N&P~4W%Uz>2P>^###tn?o;w6Z!%B7*%qG9c*Q+Hs2Z~IgC*U zj1d5i!B$-j%VD!w*xJWcp=TuSG-X}glsT-sVVj+0cF7y2Pm__$=-KH94WBvIII`oR z5#VsEd`UZM8C+~gFJnXcB?4amBdOmL@aqR3fj19umBzD@=J*lqo2f&_5FJT&ZYc|awSz`E|GYOXzCSi))ZQ)Go#{_`1^c~#) z!2$4tkARQB{~-e7)Lq^Gr$&efeUrLR-E;W=$I@speEBl@T(VtD)33*~IpdJ3YbIcm z-x_o;2h-}a151;a^N(G9+sD$qb4+qLbrwx&KZ5BH|4?XlMQiSglfy(n)r3~(Kc7zH zX+=#R1w1ES2D<~+Ue+7}IKvLNwW_ofzEzj_`o7!ovsDlJk|uUBsTTLt*YVV6<>hpZ zqw%eg;qawPq_Pgej_9Ox^Qtdz=_QZn(gEqPK-1eqxn&Fg4 z{#%KZmqh$1H5LpB^xxs)2Oj|+ zfj1a|5LgY8BPMg%4oZ>=a2RCc}S0klS9XfvoW!ED%2my%h{!y-FS#XCyo|V!;y`UwMfv zVFjmV2yZ~fkIsw-!g4zj>l#0B5ZK5otJt^#wcr z2kT^j_#ysJB9t@b)$bZH^A0ic8qlu)W`J%G{Hu?EkH9NK0FVFv{{Jh3-+!@>z$--n z_kWT;AmG;zJ_0@h|04vBhR)%y4??F*l)NSY)&#gz*@g{%<9ih^L9|`Y6|zWY!I2UD S(~kqV3J3Tq$%p#$VfE{*(9OzuvJY}3hy?R+^U+w~?B*Sn6LwH?~DM5C!acBZk% z?##GJN|A6=x=(=Mm4}4~{staFApQd$cO@iHB_t3F4=V%@%egb2FME>hE?TOFxk+Z` zo^$TG_x|p^XU_eYl-|opnuuyGr737ANemGrNxX#+K@j8cHvxZU2*JaU8G$=V5YeCE zcxLO}48D7Wabo7tgLA+bqW^-~6ZF67Pw2lJ5f6EvIWmto;o-S{1-eo%SzT7&H88(?%T3p&?xAWPugF3g%rf+9K2PdeSB-a_1OP9(f zG2c?OJ)N2_6w5bvch@?q)@q`TRh&wzrBo$Z0Bfv>HBk}eiip%BRTG=Xmw#RxxkgM- zPo-}7R;GlO=TTj3RYgsza9gF!#pfBWkiq3Qm<-!1GDP`!YI&YvOItU&Vj;V3i!3IW zmgC7KHj!GMNh~ZS63K_R;TIE%@|BIaxrvzr_M zlu<-cK9DMV3_Ir`zMPtiV{rqws#phcker!|FU>C`6RE^PJe8VHdX*s+PbBm*=m6&g zP0~75NNix)!#&rM>pgl8`GFG88zpX-lCzDTDrBRtB1uM7u_85v2HbHeRn(G^wkM(< z_G)`iR9aH?+eBcMmam8D__WSRtsAXxa~|KJKup;zx$~f0@sNN za&q6!Wlc~tQI!PQ%;Zx`x%Gf-p;I!!RFb0o?rHBj6+OSqQvEk4=Ur$5(^F zkr7%`U>w#PqFNILt)qylMM%sg7N+Cz>A6%no?1yRt|aGY;)@G+;Ej`^#v*S-##eDp zK`Q|duO>l*=)DdMZx`8H*ZiF~$+2)a+$+jMQGqvtT@o{XQer5`C=`IEe^1c=p#K0j zKlljv2)yD5jM3o{bLd~9Bc#iaP~-HZO`-kq|B5%ZUwj{dJ_J6`_{aV~MgNwd{|Fa9 z_z3t2ylM!1jS5UoT%m*ZJOb7inp=h3PNt9t7yk+w2v3B7Z~6oVJDV@}{C|r6kf8qx z7eDw2_z0XC1TGGRC#;D89`vsdO-|UpKc-OF|EK8Zc>e$V%v9X3h>w7efCB;NjKDSj zr^czru>bcFLH&#Pi24=vl=umefu9pXc6d0nyF?Na>{=bF4;pZEfEPMiOFBXRC6kZ8 zN(DlZ2>Ja3O}N((>+StIUunsjC~K-KLa!CMbeTit^xHWOxmYNwr{!g#DWWoWyNvP$ z_}k9q5Z?Dn7kALR+zwk6G=bNS+I=#tRE5l>nj|VlV%Yq8Z<}MQqFPa;wkE;GlCE@} z+emNc$|w?H8?b-KSN330Swygp$CY7s&Cb*mjKcbGO_Q2pG%~d_J=L6^s+Om2u1syN zOx5m0*a&P@Mq)9xDargHS}$xv#aM)6M`o-k96DJLvskL4@=_Ini^yfTBD+?|m&(O7 zL{V*@uUB}5?@4SFid;66#}J^XL5raxx4{*;{2Es>pj6wi%O*5<*w*A>L(`OtLGQqu z&_rrzIJ6n@=vBCv@vhe36RzmF6#e2!2>8)Y*?(N?nkpM*x`FGtBPBI$6+-VYe8-bTpHt5PBEva7(%hcb6ULcd>E$ zcQF3L1o<2CfMnqz>^SrDgCpB3r9v{*}w_b&+ez#rxG>g2->I7Y07G^zxI=ZCmPaIR2luZUsd2KaV`*A zS|Z<{v>T~kNN=ou;_zUj^^2Tr{8*!c9%hlXlwfn{C2+zxU>~7!Hot?djMe9xgCvJB zYJ)KXz%kgWt7AEAHVa$(xGD64#GIzAqnk2^bvJCYv(7Gg!}Mt~av424{h;A9-!+cn zcxVJT+$vwtj#>s68`A68kbZ%H*Z)ZB&jkGX!AIciBJkug86d7qPap5}ULt#%)>c<$ zXWL5a{tSM%w%ZM1oj-JfWKx_;w~ASnOg=^u#DVk;=<~u1b{UET|EB50q}#5 zfRDicAp+ynUETkuMu;%|4t1Zp=kWiJrO`m>`gQV|WV@E8Uyo;V#vxVLOu#0;Ht1ds zrqyQ$mL@OfAG`XtkEMI(nB;KkESl0ffawtbP-t{SYwn7Z!$d&Ugl7A{pHAaxMOC-{ zo)a&F-3BW!YYqWiV29gUQJMst0{Z6T7HXjd|*8d+M{wa=OOR*w)B! z=-M@MN5{Iiw4;g-I`|lJU-Z1IW&4;>Jn>~ACH6`$I7U==_q;xtd~&2`;X0J!id=N5 z@4a~8d>uX=$Dk;@bfqd|wy(?Wcg-j5Nd#h~>Y@JH6A$AI+T_-=+ncUwhEpQ(nG#7a ziC8x|76^sIXoWs<*Oh48e8L0d3;c1)JQp%9e1h4MQA+sH-HF3VNtA zrC?j1A#i?PLUz7*{8RMb2>PFJ@q>?mkHA?*AOLnFNRkvi47MRSLQ=-~|AfHb<)7t> z`*rgX@DVT(ph=n-p$WL?{y*_Uf_^|frwZ`!Qp53^7XuTMlfh?&hmv?;tqNHWy)~3? zWw~3(BBFE&mAD+YR<>z)t!M_;i7*eojORD9>nxl?uvUSbk%TjnWYHRYwjaG6nO|C% zgAh+`EHBQ*@e_OqZv-5}Vcxh!(BsEQ?c9T{p6J8)3|^aIQG7n;BqzMdNpG@K3Dn^2 zT;c5)t_oI-Hj0H!cTG`_+bE-VU^NJt>q_Wh92N4;Dx+1P>6my)C*0BrES)$}I_Z{9 zV(BE7?q0qS2(7J=-C2DROwr)$dERO%W4GQ|l01cV*}YVDiX;=2;lCirEw9x;*2o77 z#Ja)P0-+l>$gXil!c!v_JdyF07bp@|aB6|@#$&8|emoFLr^)WslM#Z2Y^my57M*kH z1q!1JGEx2JiH-rBvHnM-;1Jg59Qg0AlS9Oh@P87aoGEX7&xo0Ki0RWDGj86UY5Zp& z0Uv>vg8+X2`}_Yd2YWwPAAy&R03QD&eMrEsAAAIS1pY?|bkAMFUmpZdnJ9Tp0IUgc ovAhi%{KoexUV>=5l+9<5&VnN&`llamRg#bN=cCwv0D=SrOaK4? literal 0 HcmV?d00001 diff --git a/data/tests/tinshift_simplified_kkj_etrs.gpkg b/data/tests/tinshift_simplified_kkj_etrs.gpkg new file mode 100644 index 0000000000000000000000000000000000000000..aff12442a6016401294e7f2c57002bfb0496cac5 GIT binary patch literal 24576 zcmeHPUu@gP87D<2vXgjCyLv9#3BEE5WG%ETQkG>oNL$sRZ8^5&N^%>;H8@71Y&sI9 zkW^x)Z2%|zvp?91WiJDQV(WlD?5*!BwxYl=pey#Wy$oG}0Uff3t{66S{gXa)-;p9k z%Cwto{g;+!ndIH~-FJ7tyZi1BzhmK{d0CfGwV~EU9mR39|0eMFFgVSS4w!0rvmR0)brFyxcSaU;H!Xu zKPwHV zU%Z|!t-N?ue`!5E`|_bv7e9DCU3%`VcdowhW?yM0us?inxc?Xiv)7gr7+#d>s+BrQZsk&V5WYqF|M~>Y6C0aPysyfRDf}j=+8@ z)Msgi+o>?=Qa^!yYN!)J`ThSEuWUd2J_20`+)V$+^Zx+#TY~xUkK`LxcJUTz(?Rq zgTTR_(15)VzI83um7)oX$tNa#7DqKz=42sMc`Wh2l@k#!%O`) zh`@E?jlgdLPZC#%8Ti}5W%lkJTuG9ItW=~8?Q9K92L!RDH{>1oKj+}f-wXr?!(sBN z^SXGtCapHltP16ZqDzXdxjfWVo=+EfR7^iO&m$KJ4I6OAThZrOR5=pDVup{iDxR3R#xSvE;kf2{5U_IUYakW zaG0rynl6;rU`tv;$01k6%mUv@)M6ZlUFxbX*QMd`;nK+A`pDr*@$l^A;f2Y=)l*?6 ztTfJtBN3)9E5Zgko?C)~F$nuMEni*SaFW1hlQdNmhaY z)v}-pZI0<&p3lr=F$8GX42z&VKh5X)>=a)xq12k9E{iomm8ycasabpxXdU<^G?D7r zJGc<`XjQnC39i}@IG1-%in{TrtSf5IC5j9VCKBZ3NT+dh<HTuuxM>-JIC0)dg-N453?B5QGMMA^-5bCrgXQx5wrD86V zg}4iRw#Wcqy&`DkMt3JdR+|swgf&SITfDlcR7ABRv{bpPpja4m4v#L`k)}%<`shZz zW{PW<%PptF9Nrhn1b5JWopy@nRFbZ4LjBZQbwz8p*Y4on9Tn*=qRjT&Eh_jVu_ri~ zOp?zGbt^iYe8Ed;AWFNPiZ0doNrGuaZ>GIpDm0BK`y1&iN9U(S!0bOJ021N z4!6h~)T6e-#ftQ(-Dz*`0q_5j)cpkfeCH$JBXDyFd~ci#68DUZTnV*@rZrt}YLla* zO|`Ke#UpmJS(D4ga6Y=OHI&Ebs$4UNC3+Ik=Bv3@&0r8PiDnROTH^|^U8lh`0Q0fY z7$)FJ1@>m~9sCZChuh2$+H9S+(i9+DMhq^-PK?BoBQdrZi%lB8Q&7}F-2Yuxy*F2W zpW8>kN8lDlpg*u|%>M&@M2MOTtOrgz=Kq)FeZj$_N6Dw<&eYPH*W=Zkq)JNJTvl%n z{@V8G<#r}xbzn>Kn*8IberJwl%+7I;V^U{_sg3in91_lnwU%VBU2$|+2xz)kZ~p(w zX}qkc7<$@s;boGWBH*Gs7;uIi)7G+DSA}|27OH356Tgn=?Q=;ByJ5K!@s!u}lxOE< zcZvHVi+y_sj~pRO2G;GZ9Zfph!ZqCGJ>aU@F5)mAqZo*ZXQj6}ddAku0b?&ea{0?~411COS#;A3gDz#q>nEx48$ zP7dy4IpSU%k*)at;9w|3UP>9|!^vS&J?nm(yHr;2|5zZr?4i*jtezcZ+s* z8#L^}r01Oov?-fil!!V12m3$Y_RrB){YM{xTNiF`+y+-{MZocyo@Dccg5#R!Qi0NQ(W|1IEeX!jI2WN>xuYYjgYd@S#ufL>U zeQ~9fPVal6fAsF3{W@Jb_3lOf`QLNo_Bhusu0Oh@&!*pe;p#hYJqL06`hPe4pI4qu zm#!Z9*4sZ`U z>iti+(7rey0Uv=cECPE$SM>Ig0ki*qNKhYsVYBLI%tydS;FCvyA}FE{_x~Q^M+9{? za4nF7htJhrzVBdgU}&iK>D)P4I&bd^*$?eKl!r6?BghU!=>jV7^ZZn?6NdMSqOebd zX}OsxVLEf10V@Q17s%=1oSry?rqYGt@X7F4G7$sECQ`8!eu%NLxN*l%rW-SFVf6Up zNNb%2i=5$eY!vU!Fo+$CIDS5nL238cpYO3+|Gld49Tx=3qAn zS^G-p9E);UXOYn{2s1G8;^y4k9OmYB0#sMiH-@Jx&KF`zzFLm zJN`4r$sXca{3Q{*oH5+o)GvNHH%;C%_v+i1=U)5gciO@J%}2mT;0uZX9{*ti;5#1y zAAzqR0=WNE#4`kSl{yBu&&xfU)&ph)ngR3rr3HvAs$f^FH_o)vfKRl@ z>!P|2Rs&M4E(&aFe4J(qi?e(_mpQ&nXW$5d3RwZ;EU|{z6gNgQ~#<`>!IvL|QdS%6fGR|^wOupP9`}fvVt2~Jn4O2*94wvh;Yxdy2p%H_dWY&uV73h=_fiF~&{kEkxh9uU zAz#SmvV6;;*$+-8!;ar<*0vc80=s@_;L9IbE_jp}AcxL$QkjAVD`^0!mB>f7aQgvu zn086T*#8Tx5%BY!kHDvfz~wuiul5}Q1^=YH6?$Z*fF{OpH>gNpB!hM2!E2UcHi@=w z1~2H}3DB?|{cEvcoy*3PpmQgpTq4PF@x-{7=1aoPjnV11>E-b_H;yTvq4qU2%zRsU z(u_;-%r}%(L#fitPO4|j%yG@)9JTZ?n8aBlq=sY#Riv_9hvCT_BQyl(iS#uIwXxTl zYm(ZKjW_<-k{AqzB}Fq#g{M^b44QmSf+4qFZzvi{x4&41k0gynW5_l9;mD5DY}h)v z2;*jZgQ{)WjG5q^;2Mk*`X;bEW?M3+tx#mxHUnQvtzjY9kQ(f$r46v#t4Yx^*gIh% zl?r^HO~VZ7_9KK!bx}4V;g704;EeJFP;0Vs2FF-_j9vpBYYE3z+78xH@pPl5N7wXv z4JV{ZhCm|_`|(*89Md>6WHG@22nO3YFE}m@iS3?5rEy-VHN=V|0j<$e%Z4qtm4ee^ z{DRyx61UT7qh;1lJFaK4<#r_kZnzGyDf(!;j_KBV_9~>CBWoKrnu->B9Ss&Jch+6t-4_di7gXBPW%aD{HHbu`1#I9z(?RS zLSSpp?fCR@?=CM&yj}pCjN^%>=3kX`IWX2Mu zkW^x)$uK8bvX^#Sfejl54DEnnz+Tpey|>t~w{8!^9{R8X!#3opeb^p`qIu}PBSngo zX*ce+WXX7OMBd$Z-`)N0`##0An5JwiV*e(mxB$mpRy#M#eL_o;tUe>Wo@^nQ|$fRDg`4}oWn z0q_1J@N=jae)@W_Zj#g|@Jl>vUhu6Sd<1+1o_GWX-m2h-JVgx=)UT)wxcI?Gz(>GB zAVVD<9vbdHMd6^UvY=FJl2(xfy`@T;oe+&gr^Z-zERra(iP_lnY;0yc7L8tr{Ajsy z^h0*!?aEK?l6T%r1rCQ^2=$-B+9RF}ermXrLizgt zm=Ct!d>?@>1oorW)tenU`yg&#lo2>1v*ZU}rS5F8$Op6czaDXLPVzLd|br1Cjf zcRWu9hXw{A-dlJK7t=V_z3t2JTV9y>lqrb7Xn!6pXnJM=*<0b zO2E|rdwBhS?}_QS-w_`H9{~pf&J}@c{U7KLT!)wX9}t00hz|nq1a1>=5GnY#OGqC) zICwEm60#yno7&|XxDN0_OK-@#&_CkjTVD(W2SXw9TUT{qttM5Q8&$sCP;^PrHCKe1 zD{#pYhf2va84kI!(5O+CSA@ERO5FJp%H`qzawdb83h8XJu!7EUD~u@U0k+ZWS+1ZiGg%A@`8dpQ%FjJQmeiO~-m!V?}!m%wY zRu?v%SrD_cG*#ne5y2FZOK}BeE}tuw3Q35f$_8I8^D5t_SjZQ+bSj4-K%-__7!|k$ zuE6EyxS|QA))aJEsPU>);k7Nzq6>p`;Fr)uqUYdXHsm2ySju?UXz)>2^gx1o@S?1% zXwOZG3=U3Bk+;Gf;pl3qQJ^a=uv1jJC`=%B$_wL0wKJqnIclg?L6NTkV3uT9;>I|$ zK#Uuz2w*yf$z}-BRgo3hVBr?o9dZ}$r+y3iKSYqfCa;nV+zb%|w$|P}(rf6?tyek_ z)g@iP!fs$=MGo%-#mYjXhY;$>l2d6=a=DaG=b&ts%as_QtBbr=ZgeXVvSdDlE7l}E zWYOw^A_}U=w^X^SqnH^ehlfgbxarcSKCxM^ne5vAa`WjBhYp3)!CkChM^5paO48LW z7@u0Ju4rv}?N099)sXJN%9P)3R>9lEf#6^~PTn2v3{p3dF<9M+2X_uyH_7hCk1;#w zVHO#i2{wmb1iR<~a}AZ!xfN_>Y?-g`OmY~bCKw|C9D}X71j}Kw8RTeb*lEgIhADGc zchfdI`|P4OO`onHm(lB#-)Z=)XN{vd9vlG|GlYptT0?o2uE;g> zNk`8j+I-vBsu=`4&Y~Gan-0D_YyoO8v%qxhICOOo>%YsZ_mTJCZ-9@0kHEu1pg*u` z%>M&@#1OR@xD;4(%>Qr6hk}DAPm*`!&eYPH*W=Zkq)JNJTvqQ)ero6Rawn0oIu&GArpsqf(6i4)|C zfpvRpN0TnM@CbKBkGMv*D{~Ytd>JT-XQf9SV!VCvh_RTwxn)%0a!PSUj=9*kpJ(8A zhh5@mGz!l$S){G;bxr$S>&ea{0@1R#i64Ei;A3f&<>u4NS=TbdNfG@_ikKHgcsq7D zI5;#!-b@(%vtG7Y*Q4%US;?QXLwMOkqa|29I~Gfb)*3CfEb-T9*n>%CoCxTGT9x!2 zIbn1D5B`6?@1H|f{aYV_#}|RGlD#iL>reF&U?!SEc|)j5p6?mkHCKzfoR|Wu@DTVmI(6AJ0w0}YlipVIGea<>fT@3bGSXvWSXY zhMOyO((qoXWw?bBT7=yoWbG@VD=f*T{FJ$08PRS_G}P-pBLUf(jOd5CdutD?~V{`WXqywTa28iAD}S?Aw6nb zJkc?MGxz_91UO;sQ}OT4lRdts{24+0m^uX)Km0EuaO?P6mA(_iK;X8#J@iVdh-N17 zQb3g8%P-tgF?H--ii1^7ag{bRX$H;>EN6`<&n9LjX{NYzmMi4b^Bt1uSbTeRu4AC! zH1b?NyEM;ZLH;~UXpqGG;yy~#jHNeFovk$I09-hN)7WUY7{1OVUlPH;X# z^t4l<(r}f5e0R_$`zZL0p*H&XhJ;F7UvzxHqAz34_T#lv|-Y2Lw*r9 zRnqy~DqT_qMT5NsP2#hX3eKkul?KI67Lr-6#1&R)OOg0Sg>R^GRaQW?T@+>TiBjPd zP+d}G5k!|}qH!02Ru;gGwIY-;bB1;0RfHN;Ey1wJN|k0}#@39Rze%DrRDIo`fn3wl zfajD7wu1E-hy?8AbnHsbF6*)@oK@>?RwWus7=!jQS11|6+mUOUF%sM6_^edZWMxAJ zA(@Uq|09tpk48R54#?qFl|cF0dtPnNW=sS}57$B1^erHH z)$T#kN<~JWfUc$1Fca)Z4Sv+pCb$CDr13JiP+=yb2p@CPFhR0?3#n3Fkc~q4qim<0 zMQ%WD4Mq`{SbddV2eWN4$6neF(g|U$(bC7)^?D6g1OtN!jYJ&B$767UaSX{Kg2xd! zY~!}zvNTk-cNJpes#0qRqN4)S3Cvl8(_{Ra+%&T6dfI52)MVr^>z#&N>Y3~Ts5v#0wL4rKB5Y@-u>%VTmJj(9oRWhL(hP@CEw2$U}7L19^T z3Z2hp5cosF#7jgL8+dG2;RqU-9Ak!{%U~(1mW0EHTuVo=~#-SBb59|MHk*o_Bz zoEy)7YLkFJKlljv2s|VNwhugwPapU0@uI}*1%UkmtWa7`X4>b~y$I2AF`Y{xBMV%{ UjicY2@GETMUs>5Q?zh7K0(USEMF0Q* literal 0 HcmV?d00001 diff --git a/docs/source/install.rst b/docs/source/install.rst index 0b8df3427e..92ce54c103 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -459,7 +459,7 @@ All cached entries can be viewed using ``cmake -LAH`` from a build directory. .. versionadded:: 9.6 - Embed files from ending with .tif, .json or .pol in the PROJ library itself. + Embed files from ending with .tif, .json, .gpkg or .pol in the PROJ library itself. The pointed directory can potentially be the full PROJ-data package (uncompressed). In that case, about 6 GB of free disk and 16 GB of RAM are required to build PROJ. @@ -604,3 +604,8 @@ Run PROJ tests cd c:\dev\PROJ cd _build.vs2019 ctest -V --build-config Release + + +.. spelling:word-list:: + + gpkg diff --git a/docs/source/operations/transformations/tinshift.rst b/docs/source/operations/transformations/tinshift.rst index 07700dc740..08f7efade3 100644 --- a/docs/source/operations/transformations/tinshift.rst +++ b/docs/source/operations/transformations/tinshift.rst @@ -20,10 +20,11 @@ Triangulated Irregular Network based transformation The ``tinshift`` transformation takes one mandatory -argument, ``file``, that points to a JSON file, which contains the +argument, ``file``, that points to a :ref:`TIN JSON ` or +:ref:`TIN GeoPackage ` file, which contains the triangulation and associated metadata. Input and output coordinates must be geographic or projected coordinates. -Depending on the content of the JSON file, horizontal, vertical or both +Depending on the content of the TIN file, horizontal, vertical or both components of the coordinates may be transformed. The transformation is invertible, with the same computational complexity than @@ -37,7 +38,7 @@ Required .. option:: +file= - Filename to the JSON file for the TIN. + Filename or URL to the JSON or GeoPackage file for the TIN. Example @@ -116,7 +117,14 @@ the inverse transformation), otherwise which triangle will be selected is unspecified. Besides that, the triangulation does not need to have particular properties (like being a Delaunay triangulation) -JSON File format -++++++++++++++++ +JSON TIN File format +++++++++++++++++++++ See :ref:`tin_json`. + +GeoPackage TIN File format +++++++++++++++++++++++++++ + +.. versionadded:: 9.8.0 + +See :ref:`tin_gpkg`. diff --git a/docs/source/specifications/index.rst b/docs/source/specifications/index.rst index 2dfff9f7b1..e39b734dde 100644 --- a/docs/source/specifications/index.rst +++ b/docs/source/specifications/index.rst @@ -12,5 +12,6 @@ for the sake of broader interoperability. projjson geodetictiffgrids + tin_gpkg tin_json diff --git a/docs/source/specifications/tin_gpkg.rst b/docs/source/specifications/tin_gpkg.rst new file mode 100644 index 0000000000..e00f0186fe --- /dev/null +++ b/docs/source/specifications/tin_gpkg.rst @@ -0,0 +1,241 @@ +.. _tin_gpkg: + +================================================================================ +Triangulated irregular network (TIN) GeoPackage format +================================================================================ + +Introduction +------------ + +The TIN GeoPackage format defines the format for a triangulation model stored in a +`GeoPackage `__ file, that defines a horizontal and/or +vertical transformation between two CRS. + +That file format may be used by the :ref:`tinshift ` operation since +PROJ 9.8 + +Note: this format has similar capabilities than the :ref:`TIN JSON ` file format, +but can handle arbitrarily large triangulations. + +The `tin_json_to_gpkg_tin.py `__ +Python script can be used to convert a TIN JSON file into a TIN GeoPackage file. + +Specification +------------- + +The file MUST be a GeoPackage v1.2, v1.3 or v1.4 file. + +``vertices`` table +++++++++++++++++++ + +The file MUST contain a ``vertices`` GeoPackage table, of geometry type ``POINT``, +registered in the ``gpkg_contents`` table with ``data_type`` set to ``features``. + +The ``x_min``, ``y_min``, ``x_max`` and ``y_max`` columns record in the +record of ``gpkg_contents`` MUST contain the bounding box of all vertices. + +The CRS referenced by the table MUST be the source CRS of the transformation. + +The ``vertices`` table MUST have at least the following columns: + +- ``fid``: the integer primary key column. + +- ``geom``: the GeoPackage POINT geometry column. There are constraints on the + format of the GeoPackage BLOB. Both the GeoPackage header and the WKB MUST use Intel + byte order. The GeoPackage header MUST NOT have an envelope. + The point encodes the coordinates in the source CRS (equivalent of the + ``source_x`` and ``source_y`` columns in the TIN JSON format. Such columns + may also be added to the ``vertices`` table, but they are not used by the + PROJ implementation) + +If the TIN defines a horizontal shift, the following columns MUST be present: + +- ``target_x``: target X coordinate, in a column of type REAL, if the TIN defines a horizontal shift. + +- ``target_y``: target Y coordinate, in a column of type REAL, if the TIN defines a horizontal shift. + +If the TIN defines a vertical shift, the following columns MUST be present: + +- ``source_z`` and ``target_z``: source and target Z coordinate, in columns of type REAL + +- (or) ``offset_z``: value of ``target_z`` - ``source_z``, in a column of type REAL + +Additional columns may be present, but are not used by PROJ. + +The ``vertices`` table may have a RTree spatial index, but this is not used +by PROJ. + + +``triangles_def`` table ++++++++++++++++++++++++ + +The file MUST contain a ``triangles_def`` GeoPackage table. + +Its ``data_type`` registered in the ``gpkg_contents`` may be ``attributes``, +or ``features``. This is not used by PROJ. + +The ``triangles_def`` table MUST have at least the following columns: + +- ``fid``: the integer primary key column. + +- ``idx_vertex1``: foreign key to the ``fid`` column of ``vertices``, for the 1st vertex of the triangle. + +- ``idx_vertex2``: foreign key to the ``fid`` column of ``vertices``, for the 2nd vertex of the triangle. + +- ``idx_vertex3``: foreign key to the ``fid`` column of ``vertices``, for the 3rd vertex of the triangle. + +Additional columns may be present, but are not used by PROJ. + + +``gpkg_metadata`` and ``gpkg_metadata_reference`` tables +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +The file MUST contain the ``gpkg_metadata`` and ``gpkg_metadata_reference`` tables, +as defined by the `GeoPackage metadata `__ +extension. + +The ``gpkg_metadata`` table MUST contain the following record: + +- ``id`` = 1 + +- ``md_scope`` = ``dataset`` + +- ``md_uri`` = ``https://proj.org`` + +- ``mime_type`` = ``application/json`` + +- ``metadata``: serialized JSON content of the metadata described below. + +The ``gpkg_metadata_reference`` table MUST contain the following record: + +- ``reference_scope`` = ``geopackage`` + +- ``table_name`` = NULL + +- ``column_name`` = NULL + +- ``timestamp``: any ISO-8601 valid value + +- ``row_id_value`` = NULL + +- ``md_file_id`` = 1 + +- ``md_parent_id`` = NULL + + +JSON content in the ``gpkg_metadata`` table ++++++++++++++++++++++++++++++++++++++++++++ + +Required members +~~~~~~~~~~~~~~~~ + +The content of the record of id 1 in the ``gpkg_metadata`` table MUST be a +serialized JSON object, with the following required members: + +- ``file_type`` = ``triangulation_file`` + +- ``format_version`` = ``1.0`` or ``1.1`` + +- ``transformed_components``: an array which may contain one or two strings: + ``horizontal`` when horizontal components of the coordinates are transformed + and/or ``vertical`` when the vertical component is transformed. + +Optional members used by PROJ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- ``fallback_strategy`` + + String identifying how to treat points that do not fall into any of the + specified triangles. When this item is set, ``format_version`` must be set to >= 1.1. + Possible values are ``none``, ``nearest_side`` and ``nearest_centroid``. The + default is ``none`` and signifies, that points that fall outside the + specified triangles are not transformed. This is also the behavior for + ``format_version`` before 1.1. If ``fallback_strategy`` is set to + ``nearest_side``, then points that do not fall into any triangle are + transformed according to the triangle closest to them by euclidean distance. + If ``fallback_strategy`` is set to ``nearest_centroid``, then points that do + not fall into any triangle are transformed according to the triangle with the + closest centroid (intersection of its medians). + +Conditional members +~~~~~~~~~~~~~~~~~~~ + +When ``transformed_components`` contains the ``horizontal`` string, the JSON +metadata MUST contain the additional members: + +- ``min_shift_x``: minimum value of ``target_x`` - ``source_x`` over all vertices of the ``vertices`` table. + +- ``max_shift_x``: maximum value of ``target_x`` - ``source_x`` over all vertices of the ``vertices`` table. + +- ``min_shift_y``: minimum value of ``target_y`` - ``source_y`` over all vertices of the ``vertices`` table. + +- ``max_shift_y``: maximum value of ``target_y`` - ``source_y`` over all vertices of the ``vertices`` table. + +Note: those fields are used when performing inverse transformations + + +If the above described ``fallback_strategy`` field is set to ``nearest_side`` or ``nearest_centroid``, +a ``num_vertices`` member MUST be set with the number of vertices of the ``vertices`` table. + +Note: this is used to compute the initial search radius in the inverse transformation for those modes. + +Other optional members +~~~~~~~~~~~~~~~~~~~~~~ + +All other members defined in the :ref:`TIN JSON ` file format (except ``vertices``, ``vertices_columns``, +``triangles`` and ``triangles_columns`` which are redundant with the GeoPackage ``vertices`` and ``triangles_def`` +tables) may be set. + +``rtree_triangles_geom`` RTree virtual table +++++++++++++++++++++++++++++++++++++++++++++ + +A RTree virtual table MUST be created as following: + +.. code-block:: sql + + CREATE VIRTUAL TABLE rtree_triangles_geom USING rtree(id, minx, maxx, miny, maxy) + + +It MUST be filled with as many records as there are triangles in the ``triangles_def`` +table, with the ``id`` column of ``rtree_triangles_geom`` containing the +corresponding value of the ``fid`` column of ``triangles_def``, and the +``minx``, ``maxx``, ``miny``, ``maxy`` columns containing the bounding box of +each triangle. + +Optional ``triangles`` spatial view ++++++++++++++++++++++++++++++++++++ + +The file MAY contain an optional ``triangles`` spatial view, extending the +``triangles_def`` table with the GeoPackage geometry blob of each triangle. + +Such spatial view can be used to display the triangles in a GIS software that +supports GeoPackage (and spatial views), such as QGIS. + +The `tin_json_to_gpkg_tin.py `__ +script creates such view with: + +.. code-block:: python + + srs_id_i32le = as_i32le_hex(srs_id) + wkb_polygon_i32le = as_i32le_hex(3) + number_rings_i32le = as_i32le_hex(1) + number_vertices_i32le = as_i32le_hex(4) + triangle_gpkg_prefix = f"47500001{srs_id_i32le}01{wkb_polygon_i32le}{number_rings_i32le}{number_vertices_i32le}" + + # other_fields is typically ",v1.target_x, v1.target_y, v2.target_x, v2.target_y, v3.target_x, v3.target_y" for a horizontal transformation + + # 14 = GPKG_header_size_without_envelope (8) + WKB point header (5) + base_one_index (1) + ds.ExecuteSQL( + f"CREATE VIEW triangles AS SELECT triangles_def.fid AS OGC_FID {other_fields}, CAST(X'{triangle_gpkg_prefix}' || substr(v1.geom, 14) || substr(v2.geom, 14) || substr(v3.geom, 14) || substr(v1.geom, 14) AS BLOB) AS geom FROM triangles_def LEFT JOIN vertices v1 ON idx_vertex1 = v1.fid LEFT JOIN vertices v2 ON idx_vertex2 = v2.fid LEFT JOIN vertices v3 ON idx_vertex3 = v3.fid" + ) + ds.ExecuteSQL( + f"INSERT INTO gpkg_contents (table_name, identifier, data_type, srs_id, min_x, min_y, max_x, max_y) VALUES ('triangles', 'triangles', 'features', {srs_id}, {min_x}, {min_y}, {max_x}, {max_y})" + ) + ds.ExecuteSQL( + f"INSERT INTO gpkg_geometry_columns (table_name, column_name, geometry_type_name, srs_id, z, m) values ('triangles', 'geom', 'POLYGON', {srs_id}, 0, 0)" + ) + + +.. spelling:word-list:: + + RTree diff --git a/docs/source/specifications/tin_json.rst b/docs/source/specifications/tin_json.rst index e898126d39..910421f7c8 100644 --- a/docs/source/specifications/tin_json.rst +++ b/docs/source/specifications/tin_json.rst @@ -10,6 +10,9 @@ file. That file format may be used by the :ref:`tinshift ` operation since PROJ 7.2 +Note: the :ref:`TIN GeoPackage ` file format can also be used since +PROJ 9.8. It can handle arbitrarily large triangulations. + Below a minimal example, from the KKJ to ETRS89 transformation, with just a single triangle: diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index e6aa37c4a5..a0e5b3ed86 100644 --- a/scripts/reference_exported_symbols.txt +++ b/scripts/reference_exported_symbols.txt @@ -831,7 +831,7 @@ pj_ctx::pj_ctx(pj_ctx const&) pj_ctx::set_ca_bundle_path(std::string const&) pj_ctx::set_search_paths(std::vector > const&) pj_ell_set(pj_ctx*, ARG_list*, double*, double*) -pj_find_file(pj_ctx*, char const*, char*, unsigned long) +pj_find_file(pj_ctx*, char const*, char*, unsigned long, bool) pj_fwd(PJ_LP, PJconsts*) pj_get_datums_ref() pj_get_default_ctx() diff --git a/src/filemanager.cpp b/src/filemanager.cpp index bc36c18ec5..1b18251fca 100644 --- a/src/filemanager.cpp +++ b/src/filemanager.cpp @@ -1936,10 +1936,12 @@ NS_PROJ::FileManager::open_resource_file(PJ_CONTEXT *ctx, const char *name, * will receive the full filename on success. * Will be zero-terminated. * @param out_full_filename_size size of out_full_filename. + * @param disable_network whether network access must be disabled * @return 1 if the file was found, 0 otherwise. */ int pj_find_file(PJ_CONTEXT *ctx, const char *short_filename, - char *out_full_filename, size_t out_full_filename_size) { + char *out_full_filename, size_t out_full_filename_size, + bool disable_network) { const auto iter = ctx->lookupedFiles.find(short_filename); if (iter != ctx->lookupedFiles.end()) { if (iter->second.empty()) { @@ -1952,11 +1954,11 @@ int pj_find_file(PJ_CONTEXT *ctx, const char *short_filename, } const bool old_network_enabled = proj_context_is_network_enabled(ctx) != FALSE; - if (old_network_enabled) + if (disable_network && old_network_enabled) proj_context_set_enable_network(ctx, false); auto file = NS_PROJ::FileManager::open_resource_file( ctx, short_filename, out_full_filename, out_full_filename_size); - if (old_network_enabled) + if (disable_network && old_network_enabled) proj_context_set_enable_network(ctx, true); if (file) { ctx->lookupedFiles[short_filename] = out_full_filename; diff --git a/src/grids.cpp b/src/grids.cpp index ca426d2445..00c468aa21 100644 --- a/src/grids.cpp +++ b/src/grids.cpp @@ -3946,7 +3946,8 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) { /* full path of grid */ if (!pj_find_file(ctx, gridname, grinfo.filename, - sizeof(grinfo.filename) - 1)) { + sizeof(grinfo.filename) - 1, + /* disable_network = */ true)) { // Can happen when using a remote grid grinfo.filename[0] = 0; } @@ -4021,7 +4022,8 @@ PJ_INIT_INFO proj_init_info(const char *initname) { memset(&ininfo, 0, sizeof(PJ_INIT_INFO)); file_found = - pj_find_file(ctx, initname, ininfo.filename, sizeof(ininfo.filename)); + pj_find_file(ctx, initname, ininfo.filename, sizeof(ininfo.filename), + /* disable_network = */ true); if (!file_found || strlen(initname) > 64) { if (strcmp(initname, "epsg") == 0 || strcmp(initname, "EPSG") == 0) { const char *val; diff --git a/src/init.cpp b/src/init.cpp index 6567af1701..119eb4877f 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -247,11 +247,13 @@ static paralist *get_init(PJ_CONTEXT *ctx, const char *key, if (strncmp(xkey, "epsg:", 5) == 0) { exists = ctx->epsg_file_exists; if (exists < 0) { - exists = pj_find_file(ctx, initname, unused, sizeof(unused)); + exists = pj_find_file(ctx, initname, unused, sizeof(unused), + /* disable_network = */ true); ctx->epsg_file_exists = exists; } } else { - exists = pj_find_file(ctx, initname, unused, sizeof(unused)); + exists = pj_find_file(ctx, initname, unused, sizeof(unused), + /* disable_network = */ true); } if (!exists) { diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 8ac664b98d..655ea81c17 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -1261,7 +1261,8 @@ void DatabaseContext::Private::open(const std::string &databasePath, #ifndef USE_ONLY_EMBEDDED_RESOURCE_FILES path.resize(2048); const bool found = - pj_find_file(pjCtxt(), "proj.db", &path[0], path.size() - 1) != 0; + pj_find_file(pjCtxt(), "proj.db", &path[0], path.size() - 1, + /* disable_network = */ true) != 0; path.resize(strlen(path.c_str())); if (!found) #endif @@ -3498,7 +3499,8 @@ bool DatabaseContext::lookForGridInfo( bool gridAvailableWithNewName = pj_find_file(ctxt, proj_grid_name.c_str(), &fullFilenameNewName[0], - fullFilenameNewName.size() - 1) != 0; + fullFilenameNewName.size() - 1, + /* disable_network = */ true) != 0; proj_context_errno_set(ctxt, errno_before); fullFilenameNewName.resize(strlen(fullFilenameNewName.c_str())); if (gridAvailableWithNewName) { diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index c3c8f0c54c..cd5e00da89 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -12358,7 +12358,8 @@ PROJStringParser::createFromPROJString(const std::string &projString) { std::string initname(stepName); initname.resize(initname.find(':')); int file_found = - pj_find_file(ctx, initname.c_str(), unused, sizeof(unused)); + pj_find_file(ctx, initname.c_str(), unused, sizeof(unused), + /* disable_network = */ true); if (!file_found) { auto obj = createFromUserInput(stepName, d->dbContext_, true); diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index bd72c097e7..bdcd8caca4 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -159,6 +159,7 @@ set(SRC_LIBPROJ_TRANSFORMATIONS transformations/xyzgridshift.cpp transformations/defmodel.cpp transformations/tinshift.cpp + transformations/tinshift_gpkg.cpp transformations/vertoffset.cpp ) @@ -430,7 +431,7 @@ if (EMBED_RESOURCE_FILES) endif() endif() -set(EMBED_RESOURCE_DIRECTORY "" CACHE PATH "Directory that contains .tif, .json or .pol files to embed into libproj") +set(EMBED_RESOURCE_DIRECTORY "" CACHE PATH "Directory that contains .tif, .json, .gpkg or .pol files to embed into libproj") set(FILES_TO_EMBED) if (EMBED_RESOURCE_DIRECTORY) if (NOT EMBED_RESOURCE_FILES) @@ -441,9 +442,13 @@ if (EMBED_RESOURCE_DIRECTORY) message(FATAL_ERROR "${EMBED_RESOURCE_DIRECTORY} is not a valid directory") endif() - file(GLOB FILES_TO_EMBED "${EMBED_RESOURCE_DIRECTORY}/*.tif" "${EMBED_RESOURCE_DIRECTORY}/*.json" "${EMBED_RESOURCE_DIRECTORY}/*.pol") + file(GLOB FILES_TO_EMBED + "${EMBED_RESOURCE_DIRECTORY}/*.tif" + "${EMBED_RESOURCE_DIRECTORY}/*.json" + "${EMBED_RESOURCE_DIRECTORY}/*.gpkg" + "${EMBED_RESOURCE_DIRECTORY}/*.pol") if (NOT FILES_TO_EMBED) - message(FATAL_ERROR "No .tif, .json or .pol files found in ${EMBED_RESOURCE_DIRECTORY}") + message(FATAL_ERROR "No .tif, .json, .gpkg or .pol files found in ${EMBED_RESOURCE_DIRECTORY}") endif() endif() diff --git a/src/proj_internal.h b/src/proj_internal.h index 61cfbfd03c..b989a636d3 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -1065,7 +1065,7 @@ void pj_stderr_logger(void *, int, const char *); // PROJ_DLL for tests int PROJ_DLL pj_find_file(PJ_CONTEXT *ctx, const char *short_filename, char *out_full_filename, - size_t out_full_filename_size); + size_t out_full_filename_size, bool disable_network); // To remove when PROJ_LIB definitely goes away void PROJ_DLL pj_stderr_proj_lib_deprecation_warning(); diff --git a/src/sqlite3_utils.cpp b/src/sqlite3_utils.cpp index 8f578a92ae..6d80321cbc 100644 --- a/src/sqlite3_utils.cpp +++ b/src/sqlite3_utils.cpp @@ -40,6 +40,8 @@ #include "memvfs.h" #endif +#include "filemanager.hpp" + #include #include #include // std::ostringstream @@ -233,6 +235,200 @@ std::unique_ptr SQLite3VFS::create(bool fakeSync, bool fakeLock, // --------------------------------------------------------------------------- +static int VFSNetworkClose(sqlite3_file *file) { + File *pj_file; + memcpy(&pj_file, reinterpret_cast(file) + sizeof(sqlite3_file), + sizeof(pj_file)); + delete pj_file; + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkRead(sqlite3_file *file, void *pBuffer, int iAmt, + sqlite3_int64 iOffset) { + File *pj_file; + memcpy(&pj_file, reinterpret_cast(file) + sizeof(sqlite3_file), + sizeof(pj_file)); + if (!pj_file->seek(iOffset)) + return SQLITE_IOERR_READ; + const int nRead = static_cast(pj_file->read(pBuffer, iAmt)); + if (nRead < iAmt) { + memset(static_cast(pBuffer) + nRead, 0, iAmt - nRead); + return SQLITE_IOERR_SHORT_READ; + } + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkWrite(sqlite3_file *, const void *, int, sqlite3_int64) { + return SQLITE_IOERR_WRITE; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkFileSize(sqlite3_file *file, sqlite3_int64 *pSize) { + File *pj_file; + memcpy(&pj_file, reinterpret_cast(file) + sizeof(sqlite3_file), + sizeof(pj_file)); + pj_file->seek(0, SEEK_END); + *pSize = pj_file->tell(); + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkTruncate(sqlite3_file *, sqlite3_int64) { + return SQLITE_IOERR_TRUNCATE; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkSync(sqlite3_file *, int) { return SQLITE_OK; } + +// --------------------------------------------------------------------------- + +static int VFSNetworkLock(sqlite3_file *, int) { return SQLITE_OK; } + +// --------------------------------------------------------------------------- + +static int VFSNetworkUnlock(sqlite3_file *, int) { return SQLITE_OK; } + +// --------------------------------------------------------------------------- + +static int VFSNetworkCheckReservedLock(sqlite3_file *, int *pResOut) { + *pResOut = 0; + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkFileControl(sqlite3_file *, int, void *) { + return SQLITE_NOTFOUND; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkSectorSize(sqlite3_file *) { return 0; } + +// --------------------------------------------------------------------------- + +static int VFSNetworkDeviceCharacteristics(sqlite3_file *) { return 0; } + +// --------------------------------------------------------------------------- + +static const sqlite3_io_methods VFSNetworkIOMethods = { + 1, + VFSNetworkClose, + VFSNetworkRead, + VFSNetworkWrite, + VFSNetworkTruncate, + VFSNetworkSync, + VFSNetworkFileSize, + VFSNetworkLock, + VFSNetworkUnlock, + VFSNetworkCheckReservedLock, + VFSNetworkFileControl, + VFSNetworkSectorSize, + VFSNetworkDeviceCharacteristics, + nullptr, // xShmMap + nullptr, // xShmLock + nullptr, // xShmBarrier + nullptr, // xShmUnmap + nullptr, // xFetch + nullptr, // xUnfetch +}; + +static int VFSNetworkOpen(sqlite3_vfs *vfs, const char *name, + sqlite3_file *file, int /*flags*/, + int * /*outFlags*/) { + if (!name) + return SQLITE_CANTOPEN; + PJ_CONTEXT *ctx = static_cast(vfs->pAppData); + auto pjFile = pj_network_file_open(ctx, name); + if (!pjFile) + return SQLITE_CANTOPEN; + file->pMethods = &VFSNetworkIOMethods; + File *pj_file_raw = pjFile.release(); + memcpy(reinterpret_cast(file) + sizeof(sqlite3_file), &pj_file_raw, + sizeof(pj_file_raw)); + + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkFullPathname(sqlite3_vfs *, const char *zName, int nOut, + char *zOut) { + if (static_cast(strlen(zName)) >= nOut) { + fprintf(stderr, "Maximum pathname length reserved for SQLite3 VFS " + "isn't large enough.\n"); + return SQLITE_CANTOPEN; + } + strncpy(zOut, zName, nOut); + zOut[nOut - 1] = '\0'; + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkAccess(sqlite3_vfs *, const char *, int, int *pResOut) { + *pResOut = 0; + return SQLITE_OK; +} + +// --------------------------------------------------------------------------- + +static int VFSNetworkDelete(sqlite3_vfs *, const char *, int) { + return SQLITE_IOERR_DELETE; +} + +// --------------------------------------------------------------------------- + +/* static */ std::unique_ptr +SQLite3VFS::createNetwork(PJ_CONTEXT *ctx) { + // Install SQLite3 logger if PROJ_LOG_SQLITE3 env var is defined + InstallSqliteLogger::GetSingleton(); + + // Call to sqlite3_initialize() is normally not needed, except for + // people building SQLite3 with -DSQLITE_OMIT_AUTOINIT + sqlite3_initialize(); + sqlite3_vfs *defaultVFS = sqlite3_vfs_find(nullptr); + assert(defaultVFS); + + auto vfs = new pj_sqlite3_vfs(); + + std::ostringstream buffer; + buffer << vfs; + vfs->namePtr = buffer.str(); + + vfs->base.iVersion = 1; + vfs->base.szOsFile = sizeof(sqlite3_file) + sizeof(File *); + vfs->base.mxPathname = defaultVFS->mxPathname; + vfs->base.zName = vfs->namePtr.c_str(); + vfs->base.pAppData = ctx; + vfs->base.xOpen = VFSNetworkOpen; + vfs->base.xDelete = VFSNetworkDelete; + vfs->base.xAccess = VFSNetworkAccess; + vfs->base.xFullPathname = VFSNetworkFullPathname; + vfs->base.xDlOpen = defaultVFS->xDlOpen; + vfs->base.xDlError = defaultVFS->xDlError; + vfs->base.xDlSym = defaultVFS->xDlSym; + vfs->base.xDlClose = defaultVFS->xDlClose; + vfs->base.xRandomness = defaultVFS->xRandomness; + vfs->base.xSleep = defaultVFS->xSleep; + vfs->base.xCurrentTime = defaultVFS->xCurrentTime; + vfs->base.xGetLastError = defaultVFS->xGetLastError; + vfs->base.xCurrentTimeInt64 = defaultVFS->xCurrentTimeInt64; + if (sqlite3_vfs_register(&(vfs->base), false) == SQLITE_OK) { + return std::unique_ptr(new SQLite3VFS(vfs)); + } + return nullptr; +} + +// --------------------------------------------------------------------------- + #ifdef EMBED_RESOURCE_FILES struct pj_sqlite3_memvfs : public pj_sqlite3_vfs { diff --git a/src/sqlite3_utils.hpp b/src/sqlite3_utils.hpp index e592e7604d..6ae1b942ac 100644 --- a/src/sqlite3_utils.hpp +++ b/src/sqlite3_utils.hpp @@ -69,6 +69,8 @@ class SQLite3VFS { size_t bufferSize); #endif + static std::unique_ptr createNetwork(PJ_CONTEXT *ctx); + const char *name() const; sqlite3_vfs *raw() { return &(vfs_->base); } }; diff --git a/src/transformations/tinshift.cpp b/src/transformations/tinshift.cpp index 7b044b20bf..95e9af2ea0 100644 --- a/src/transformations/tinshift.cpp +++ b/src/transformations/tinshift.cpp @@ -26,19 +26,29 @@ *****************************************************************************/ #define PROJ_COMPILATION +#define FROM_PROJ_CPP -#include "tinshift.hpp" #include "filemanager.hpp" +#include "proj/internal/internal.hpp" #include "proj_internal.h" +#include "tinshift_gpkg.hpp" +#include "tinshift_iface.hpp" +#include "tinshift_json.hpp" PROJ_HEAD(tinshift, "Triangulation based transformation"); -using namespace TINSHIFT_NAMESPACE; +using namespace TINSHIFT_JSON_NAMESPACE; + +// --------------------------------------------------------------------------- + +TINShiftEvaluator::~TINShiftEvaluator() = default; + +// --------------------------------------------------------------------------- namespace { struct tinshiftData { - std::unique_ptr evaluator{}; + std::unique_ptr evaluator{}; tinshiftData() = default; @@ -85,6 +95,52 @@ PJ *PJ_TRANSFORMATION(tinshift, 1) { return pj_tinshift_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG); } + P->fwd4d = tinshift_forward_4d; + P->inv4d = tinshift_reverse_4d; + P->left = PJ_IO_UNITS_WHATEVER; + P->right = PJ_IO_UNITS_WHATEVER; + + if (NS_PROJ::internal::ends_with(filename, ".gpkg") || + NS_PROJ::internal::ends_with(filename, ".GPKG")) { + std::string path; + if (NS_PROJ::internal::starts_with(filename, "http://") || + NS_PROJ::internal::starts_with(filename, "https://")) { + path = filename; + } else { + path.resize(2048); + const bool found = + pj_find_file(P->ctx, filename, &path[0], path.size() - 1, + /* disable_network = */ false) != 0; + path.resize(strlen(path.c_str())); + if (!found) { + proj_log_error(P, _("non existing file")); + return pj_tinshift_destructor( + P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID); + } + } + + std::unique_ptr file; + try { + file = TINShiftGeopackageFile::open(P->ctx, path); + + } catch (const std::exception &e) { + proj_log_error(P, _("invalid file: %s"), e.what()); + return pj_tinshift_destructor( + P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID); + } + + auto Q = new tinshiftData(); + P->opaque = (void *)Q; + P->destructor = pj_tinshift_destructor; + + Q->evaluator = + std::make_unique(std::move(file)); + + return P; + } + + // else JSON file + auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, filename); if (nullptr == file) { proj_log_error(P, _("Cannot open %s"), filename); @@ -120,17 +176,13 @@ PJ *PJ_TRANSFORMATION(tinshift, 1) { P->destructor = pj_tinshift_destructor; try { - Q->evaluator.reset(new Evaluator(TINShiftFile::parse(jsonStr))); + Q->evaluator = std::make_unique( + TINShiftJSONFile::parse(jsonStr)); } catch (const std::exception &e) { proj_log_error(P, _("invalid model: %s"), e.what()); return pj_tinshift_destructor( P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID); } - P->fwd4d = tinshift_forward_4d; - P->inv4d = tinshift_reverse_4d; - P->left = PJ_IO_UNITS_WHATEVER; - P->right = PJ_IO_UNITS_WHATEVER; - return P; } diff --git a/src/transformations/tinshift_gpkg.cpp b/src/transformations/tinshift_gpkg.cpp new file mode 100644 index 0000000000..0713cb30e7 --- /dev/null +++ b/src/transformations/tinshift_gpkg.cpp @@ -0,0 +1,945 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to TIN based transformations using a + *GeoPackage file Author: Even Rouault, + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault, + * + * SPDX-License-Identifier: MIT + *****************************************************************************/ + +#define FROM_PROJ_CPP + +#include "tinshift_gpkg.hpp" +#include "proj/internal/include_nlohmann_json.hpp" +#include "proj/internal/internal.hpp" + +#ifdef EMBED_RESOURCE_FILES +#include "embedded_resources.h" +#endif + +#include +#include +#include +#include +#include +#include +#include // std::ostringstream + +#include + +using json = nlohmann::json; +using namespace NS_PROJ; + +// --------------------------------------------------------------------------- + +const char *TINShiftGeopackageException::what() const noexcept { + return msg_.c_str(); +} + +// --------------------------------------------------------------------------- + +std::unique_ptr +TINShiftGeopackageFile::open(PJ_CONTEXT *ctx, const std::string &filename) { + std::string sqliteURI(filename); + + auto poTinshift = + std::unique_ptr(new TINShiftGeopackageFile()); + if (internal::starts_with(filename, "http://") || + internal::starts_with(filename, "https://")) { + poTinshift->sqlite_vfs_ = SQLite3VFS::createNetwork(ctx); + } + + while (true) { + if (sqlite3_open_v2(sqliteURI.c_str(), &poTinshift->sqlite_handle_, + SQLITE_OPEN_READONLY | SQLITE_OPEN_URI, + poTinshift->sqlite_vfs_ + ? poTinshift->sqlite_vfs_->name() + : nullptr) != SQLITE_OK || + !poTinshift->sqlite_handle_) { + +#ifdef EMBED_RESOURCE_FILES + if (!poTinshift->sqlite_vfs_) { + unsigned int size = 0; + const unsigned char *in_memory_data = + pj_get_embedded_resource(filename.c_str(), &size); + if (!in_memory_data) { + throw TINShiftGeopackageException(std::string("Open of ") + .append(filename) + .append(" failed")); + } + + poTinshift->sqlite_vfs_ = + SQLite3VFS::createMem(in_memory_data, size); + if (poTinshift->sqlite_vfs_ == nullptr) { + throw TINShiftGeopackageException(std::string("Open of ") + .append(filename) + .append(" failed")); + } + + std::ostringstream buffer; + buffer << "file:/"; + buffer << filename; + buffer << "?immutable=1&ptr="; + buffer << reinterpret_cast(in_memory_data); + buffer << "&sz="; + buffer << size; + buffer << "&max="; + buffer << size; + buffer << "&vfs="; + buffer << poTinshift->sqlite_vfs_->name(); + sqliteURI = buffer.str(); + + continue; + } +#endif + + throw TINShiftGeopackageException( + std::string("Open of ").append(sqliteURI).append(" failed")); + } else { + break; + } + } + + poTinshift->getMetadata(); + poTinshift->prepareDatabaseQuery(); + + return poTinshift; +} + +// --------------------------------------------------------------------------- + +TINShiftGeopackageFile::~TINShiftGeopackageFile() { + if (select_stmt_) + sqlite3_finalize(select_stmt_); + if (sqlite_handle_) + sqlite3_close(sqlite_handle_); +} + +// --------------------------------------------------------------------------- + +namespace { + +static std::string getString(const json &j, const char *key, bool optional) { + if (!j.contains(key)) { + if (optional) { + return std::string(); + } + throw TINShiftGeopackageException(std::string("Missing \"") + key + + "\" key"); + } + const json v = j[key]; + if (!v.is_string()) { + throw TINShiftGeopackageException(std::string("The value of \"") + key + + "\" should be a string"); + } + return v.get(); +} + +static std::string getReqString(const json &j, const char *key) { + return getString(j, key, false); +} + +static double getReqDouble(const json &j, const char *key) { + if (!j.contains(key)) { + throw TINShiftGeopackageException(std::string("Missing \"") + key + + "\" key"); + } + + const json v = j[key]; + if (!v.is_number()) { + throw TINShiftGeopackageException(std::string("The value of \"") + key + + "\" should be a number"); + } + return v.get(); +} + +static int getReqInt(const json &j, const char *key) { + if (!j.contains(key)) { + throw TINShiftGeopackageException(std::string("Missing \"") + key + + "\" key"); + } + + const json v = j[key]; + if (!v.is_number()) { + throw TINShiftGeopackageException(std::string("The value of \"") + key + + "\" should be an integer"); + } + return v.get(); +} + +// --------------------------------------------------------------------------- + +static json getArrayMember(const json &j, const char *key) { + if (!j.contains(key)) { + throw TINShiftGeopackageException(std::string("Missing \"") + key + + "\" key"); + } + const json obj = j[key]; + if (!obj.is_array()) { + throw TINShiftGeopackageException(std::string("The value of \"") + key + + "\" should be a array"); + } + return obj; +} + +} // namespace +// --------------------------------------------------------------------------- + +void TINShiftGeopackageFile::getMetadata() { + std::string metadata; + // Get JSON metadata entry from gpkg_metadata table + { + char **papszResult = nullptr; + int nRows = 0; + char *pszErrMsg = nullptr; + if (sqlite3_get_table(sqlite_handle_, + "SELECT metadata FROM gpkg_metadata WHERE " + "md_standard_uri = 'https://proj.org'", + &papszResult, &nRows, nullptr, + &pszErrMsg) != SQLITE_OK) { + sqlite3_free_table(papszResult); + auto ex = TINShiftGeopackageException( + std::string("Cannot get metadata: ").append(pszErrMsg)); + sqlite3_free(pszErrMsg); + throw ex; + } + if (nRows != 1 || !papszResult[1]) { + sqlite3_free_table(papszResult); + throw TINShiftGeopackageException( + "Cannot get metadata in gpkg_metadata table"); + } + metadata = papszResult[1]; + sqlite3_free_table(papszResult); + } + + json j; + try { + j = json::parse(metadata); + } catch (const std::exception &e) { + throw TINShiftGeopackageException( + std::string("Cannot parse JSON metadata: ").append(e.what())); + } + if (!j.is_object()) { + throw TINShiftGeopackageException("JSON metadata is not an object"); + } + + if (getReqString(j, "file_type") != "triangulation_file") { + throw TINShiftGeopackageException("File is not of the expected type"); + } + + const std::string version = getReqString(j, "format_version"); + if (!internal::starts_with(version, "1.")) { + throw TINShiftGeopackageException(std::string("format_version = ") + .append(version) + .append(" is not supported")); + } + + const auto jTransformedComponents = + getArrayMember(j, "transformed_components"); + for (const json &jComp : jTransformedComponents) { + if (!jComp.is_string()) { + throw TINShiftGeopackageException( + "transformed_components[] item is not a string"); + } + const auto jCompStr = jComp.get(); + if (jCompStr == "horizontal") { + horizontalTransformation_ = true; + } else if (jCompStr == "vertical") { + verticalTransformation_ = true; + } else { + throw TINShiftGeopackageException( + "transformed_components[] = " + jCompStr + " is not handled"); + } + } + + if (horizontalTransformation_) { + min_shift_x_ = getReqDouble(j, "min_shift_x"); + min_shift_y_ = getReqDouble(j, "min_shift_y"); + max_shift_x_ = getReqDouble(j, "max_shift_x"); + max_shift_y_ = getReqDouble(j, "max_shift_y"); + } + + if (j.contains("fallback_strategy")) { + const auto fallback_strategy = getReqString(j, "fallback_strategy"); + if (fallback_strategy == "nearest_side") { + fallbackStrategy_ = FALLBACK_NEAREST_SIDE; + } else if (fallback_strategy == "nearest_centroid") { + fallbackStrategy_ = FALLBACK_NEAREST_CENTROID; + } else if (fallback_strategy == "none") { + fallbackStrategy_ = FALLBACK_NONE; + } else { + throw TINShiftGeopackageException("invalid fallback_strategy"); + } + } + + if (fallbackStrategy_ != FALLBACK_NONE) { + numVertices_ = getReqInt(j, "num_vertices"); + if (numVertices_ <= 0) { + throw TINShiftGeopackageException("invalid value for num_vertices"); + } + } + + // Get bounding box of vertices table + { + char **papszResult = nullptr; + int nRows = 0; + char *pszErrMsg = nullptr; + if (sqlite3_get_table(sqlite_handle_, + "SELECT min_x, min_y, max_x, max_y FROM " + "gpkg_contents WHERE table_name = 'vertices'", + &papszResult, &nRows, nullptr, + &pszErrMsg) != SQLITE_OK) { + sqlite3_free_table(papszResult); + auto ex = TINShiftGeopackageException( + std::string("Cannot get bounding box of vertices table: ") + .append(pszErrMsg)); + sqlite3_free(pszErrMsg); + throw ex; + } + if (nRows != 1 || !papszResult[4] || !papszResult[5] || + !papszResult[6] || !papszResult[7]) { + sqlite3_free_table(papszResult); + throw TINShiftGeopackageException( + "Cannot get bounding box of vertices table"); + } + min_x_ = internal::c_locale_stod(papszResult[4]); + min_y_ = internal::c_locale_stod(papszResult[5]); + max_x_ = internal::c_locale_stod(papszResult[6]); + max_y_ = internal::c_locale_stod(papszResult[7]); + sqlite3_free_table(papszResult); + if (!(min_x_ < max_x_) || !(min_y_ < max_y_)) { + throw TINShiftGeopackageException( + "Invalid bounding box of vertices table"); + } + } +} + +// --------------------------------------------------------------------------- + +void TINShiftGeopackageFile::prepareDatabaseQuery() { + // Get field names of vertices table + { + char **papszResult = nullptr; + int nRows = 0; + int nCols = 0; + char *pszErrMsg = nullptr; + if (sqlite3_get_table(sqlite_handle_, "PRAGMA table_info(vertices)", + &papszResult, &nRows, &nCols, + &pszErrMsg) != SQLITE_OK) { + sqlite3_free_table(papszResult); + auto ex = TINShiftGeopackageException( + std::string("Cannot get definition of table vertices: ") + .append(pszErrMsg)); + sqlite3_free(pszErrMsg); + throw ex; + } + if (nRows == 0 || nCols != 6) { + sqlite3_free_table(papszResult); + throw TINShiftGeopackageException( + "Cannot get definition of table vertices"); + } + for (int i = 0; i < nRows; ++i) { + const char *pszColName = papszResult[(i + 1) * nCols + 1]; + if (pszColName && strcmp(pszColName, "fid") != 0 && + strcmp(pszColName, "geom") != 0) { + std::string osColName(pszColName); + if (osColName == "source_z") + sourceZ_ = true; + else if (osColName == "target_x") + targetX_ = true; + else if (osColName == "target_y") + targetY_ = true; + else if (osColName == "target_z") + targetZ_ = true; + else if (osColName == "offset_z") + offsetZ_ = true; + vertices_cols_.push_back(std::move(osColName)); + } + } + sqlite3_free_table(papszResult); + } + + if (horizontalTransformation_ && !targetX_) + throw TINShiftGeopackageException( + "target_x field missing in table vertices"); + if (horizontalTransformation_ && !targetY_) + throw TINShiftGeopackageException( + "target_y field missing in table vertices"); + if (verticalTransformation_ && !((sourceZ_ && targetZ_) || offsetZ_)) + throw TINShiftGeopackageException("(source_z and target_z) or offset_z " + "fields missing in table vertices"); + + std::vector colsToRequest; + if (horizontalTransformation_) { + colsToRequest.push_back("target_x"); + colsToRequest.push_back("target_y"); + } + if (verticalTransformation_) { + if (sourceZ_) { + colsToRequest.push_back("source_z"); + colsToRequest.push_back("target_z"); + } else { + colsToRequest.push_back("offset_z"); + } + } + + std::string sql( + "SELECT v1.geom AS v1_geom, v2.geom AS v2_geom, v3.geom AS v3_geom"); + for (const std::string &col : colsToRequest) { + sql += ", v1."; + sql += col; + sql += " AS v1_"; + sql += col; + + sql += ", v2."; + sql += col; + sql += " AS v2_"; + sql += col; + + sql += ", v3."; + sql += col; + sql += " AS v3_"; + sql += col; + } + sql += " FROM triangles_def " + "LEFT JOIN vertices v1 ON idx_vertex1 = v1.fid " + "LEFT JOIN vertices v2 ON idx_vertex2 = v2.fid " + "LEFT JOIN vertices v3 ON idx_vertex3 = v3.fid " + "WHERE triangles_def.fid IN (" + "SELECT id FROM rtree_triangles_geom " + "WHERE maxx >= ? AND minx <= ? AND maxy >= ? AND miny <= ?)"; + + sqlite3_prepare_v2(sqlite_handle_, sql.c_str(), -1, &select_stmt_, nullptr); + if (!select_stmt_) { + throw TINShiftGeopackageException(sqlite3_errmsg(sqlite_handle_)); + } +} + +// --------------------------------------------------------------------------- + +/************************************************************************/ +/* swap_words() */ +/* */ +/* Convert the byte order of the given word(s) in place. */ +/************************************************************************/ + +static const int byte_order_test = 1; +#define IS_LSB \ + (1 == (reinterpret_cast(&byte_order_test))[0]) + +static void swap_words(void *dataIn, size_t word_size, size_t word_count) + +{ + unsigned char *data = static_cast(dataIn); + for (size_t word = 0; word < word_count; word++) { + for (size_t i = 0; i < word_size / 2; i++) { + unsigned char t; + + t = data[i]; + data[i] = data[word_size - i - 1]; + data[word_size - i - 1] = t; + } + + data += word_size; + } +} + +// --------------------------------------------------------------------------- + +static inline double sqr(double x) { return x * x; } +static inline double squared_distance(double x1, double y1, double x2, + double y2) { + return sqr(x1 - x2) + sqr(y1 - y2); +} +static double sq_distance_point_segment(double x, double y, double x1, + double y1, double x2, double y2, + double dist12) { + // squared distance of point x/y to line segment x1/y1 -- x2/y2 + const double t = ((x - x1) * (x2 - x1) + (y - y1) * (y2 - y1)) / dist12; + if (t <= 0.0) { + // closest to x1/y1 + return squared_distance(x, y, x1, y1); + } + if (t >= 1.0) { + // closest to y2/y2 + return squared_distance(x, y, x2, y2); + } + + // closest to line segment x1/y1 -- x2/y2 + return squared_distance(x, y, x1 + t * (x2 - x1), y1 + t * (y2 - y1)); +} + +// --------------------------------------------------------------------------- + +bool TINShiftGeopackageFile::findTriangle( + double x, double y, bool forwardDirection, double &lambda1, double &lambda2, + double &lambda3, std::vector &valsVertex1, + std::vector &valsVertex2, std::vector &valsVertex3) { + + const bool useForwardDirection = + forwardDirection || !horizontalTransformation_; + + // Compute barycentric coefficients + const auto computeLambdas = [x, y, &lambda1, &lambda2, + &lambda3](const std::array &vX, + const std::array &vY, + bool checkWithinTriangle) { + constexpr double EPSILON = 1e-10; + const double det_T = (vY[1] - vY[2]) * (vX[0] - vX[2]) + + (vX[2] - vX[1]) * (vY[0] - vY[2]); + if (std::fabs(det_T) < EPSILON) { + // Degenerate triangle + return false; + } + lambda1 = + ((vY[1] - vY[2]) * (x - vX[2]) + (vX[2] - vX[1]) * (y - vY[2])) / + det_T; + lambda2 = + ((vY[2] - vY[0]) * (x - vX[2]) + (vX[0] - vX[2]) * (y - vY[2])) / + det_T; + lambda3 = 1 - lambda1 - lambda2; + return !checkWithinTriangle || + (lambda1 >= -EPSILON && lambda1 <= 1 + EPSILON && + lambda2 >= -EPSILON && lambda2 <= 1 + EPSILON && + lambda3 >= -EPSILON && lambda3 <= 1 + EPSILON); + }; + + // Get the current triangle's vertices coordinates + const auto getXY = [this, useForwardDirection](std::array &vX, + std::array &vY) { + if (useForwardDirection) { + const unsigned char *v1_geom = static_cast( + sqlite3_column_blob(select_stmt_, 0)); + const unsigned char *v2_geom = static_cast( + sqlite3_column_blob(select_stmt_, 1)); + const unsigned char *v3_geom = static_cast( + sqlite3_column_blob(select_stmt_, 2)); + constexpr int SIZEOF_GPKG_POINT_PREFIX = 8 + 1 + 4; + constexpr int SIZEOF_GPKG_POINT = + SIZEOF_GPKG_POINT_PREFIX + 2 * static_cast(sizeof(double)); + if (sqlite3_column_bytes(select_stmt_, 0) != SIZEOF_GPKG_POINT || + sqlite3_column_bytes(select_stmt_, 1) != SIZEOF_GPKG_POINT || + sqlite3_column_bytes(select_stmt_, 2) != SIZEOF_GPKG_POINT) { + return false; + } + + memcpy(&vX[0], v1_geom + SIZEOF_GPKG_POINT_PREFIX, sizeof(double)); + memcpy(&vY[0], v1_geom + SIZEOF_GPKG_POINT_PREFIX + sizeof(double), + sizeof(double)); + memcpy(&vX[1], v2_geom + SIZEOF_GPKG_POINT_PREFIX, sizeof(double)); + memcpy(&vY[1], v2_geom + SIZEOF_GPKG_POINT_PREFIX + sizeof(double), + sizeof(double)); + memcpy(&vX[2], v3_geom + SIZEOF_GPKG_POINT_PREFIX, sizeof(double)); + memcpy(&vY[2], v3_geom + SIZEOF_GPKG_POINT_PREFIX + sizeof(double), + sizeof(double)); + + if (!IS_LSB) { + swap_words(vX.data(), sizeof(double), vX.size()); + swap_words(vY.data(), sizeof(double), vY.size()); + } + } else { + vX[0] = sqlite3_column_double(select_stmt_, 3); + vX[1] = sqlite3_column_double(select_stmt_, 4); + vX[2] = sqlite3_column_double(select_stmt_, 5); + vY[0] = sqlite3_column_double(select_stmt_, 6); + vY[1] = sqlite3_column_double(select_stmt_, 7); + vY[2] = sqlite3_column_double(select_stmt_, 8); + } + return true; + }; + + // Get the values at the vertices of the current triangle + const auto collectValsVertex = [this, useForwardDirection, &valsVertex1, + &valsVertex2, &valsVertex3]() { + if (useForwardDirection) { + int nextAttrIdx = 3; + + if (horizontalTransformation_) { + // target_x + valsVertex1.push_back(sqlite3_column_double(select_stmt_, 3)); + // target_y + valsVertex1.push_back(sqlite3_column_double(select_stmt_, 6)); + + // target_x + valsVertex2.push_back(sqlite3_column_double(select_stmt_, 4)); + // target_y + valsVertex2.push_back(sqlite3_column_double(select_stmt_, 7)); + + // target_x + valsVertex3.push_back(sqlite3_column_double(select_stmt_, 5)); + // target_y + valsVertex3.push_back(sqlite3_column_double(select_stmt_, 8)); + + nextAttrIdx = 9; + } + + if (verticalTransformation_) { + if (sourceZ_) { + // source_z + valsVertex1.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx)); + // target_z + valsVertex1.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 3)); + + // source_z + valsVertex2.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 1)); + // target_z + valsVertex2.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 4)); + + // source_z + valsVertex3.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 2)); + // target_z + valsVertex3.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 5)); + } else { + // offset_z + valsVertex1.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx)); + valsVertex2.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 1)); + valsVertex3.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 2)); + } + } + } else { + int nextAttrIdx = 3; + + if (horizontalTransformation_) { + + const unsigned char *v1_geom = + static_cast( + sqlite3_column_blob(select_stmt_, 0)); + const unsigned char *v2_geom = + static_cast( + sqlite3_column_blob(select_stmt_, 1)); + const unsigned char *v3_geom = + static_cast( + sqlite3_column_blob(select_stmt_, 2)); + constexpr int SIZEOF_GPKG_POINT_PREFIX = 8 + 1 + 4; + constexpr int SIZEOF_GPKG_POINT = + SIZEOF_GPKG_POINT_PREFIX + + 2 * static_cast(sizeof(double)); + if (sqlite3_column_bytes(select_stmt_, 0) != + SIZEOF_GPKG_POINT || + sqlite3_column_bytes(select_stmt_, 1) != + SIZEOF_GPKG_POINT || + sqlite3_column_bytes(select_stmt_, 2) != + SIZEOF_GPKG_POINT) { + return false; + } + double vX[3], vY[3]; + memcpy(&vX[0], v1_geom + SIZEOF_GPKG_POINT_PREFIX, + sizeof(double)); + memcpy(&vY[0], + v1_geom + SIZEOF_GPKG_POINT_PREFIX + sizeof(double), + sizeof(double)); + memcpy(&vX[1], v2_geom + SIZEOF_GPKG_POINT_PREFIX, + sizeof(double)); + memcpy(&vY[1], + v2_geom + SIZEOF_GPKG_POINT_PREFIX + sizeof(double), + sizeof(double)); + memcpy(&vX[2], v3_geom + SIZEOF_GPKG_POINT_PREFIX, + sizeof(double)); + memcpy(&vY[2], + v3_geom + SIZEOF_GPKG_POINT_PREFIX + sizeof(double), + sizeof(double)); + + if (!IS_LSB) { + swap_words(vX, sizeof(double), 3); + swap_words(vY, sizeof(double), 3); + } + + // source_x + valsVertex1.push_back(vX[0]); + // source_y + valsVertex1.push_back(vY[0]); + + // source_x + valsVertex2.push_back(vX[1]); + // source_y + valsVertex2.push_back(vY[1]); + + // source_x + valsVertex3.push_back(vX[2]); + // source_y + valsVertex3.push_back(vY[2]); + + nextAttrIdx = 9; + } + + if (verticalTransformation_) { + if (sourceZ_) { + // source_z + valsVertex1.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx)); + // target_z + valsVertex1.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 3)); + + // source_z + valsVertex2.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 1)); + // target_z + valsVertex2.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 4)); + + // source_z + valsVertex3.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 2)); + // target_z + valsVertex3.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 5)); + } else { + // offset_z + valsVertex1.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx)); + valsVertex2.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 1)); + valsVertex3.push_back( + sqlite3_column_double(select_stmt_, nextAttrIdx + 2)); + } + } + } + return true; + }; + + // Try to use last triangle + if (cachedForwardDirection_ == forwardDirection && + !cachedValsVertex1_.empty() && + computeLambdas(cachedVerticesX_, cachedVerticesY_, + /* checkWithinTriangle = */ true)) { + valsVertex1 = cachedValsVertex1_; + valsVertex2 = cachedValsVertex2_; + valsVertex3 = cachedValsVertex3_; + return true; + } + + constexpr double EPSILON = 1e-10; + + sqlite3_reset(select_stmt_); + if (useForwardDirection) { + sqlite3_bind_double(select_stmt_, 1, x - EPSILON); + sqlite3_bind_double(select_stmt_, 2, x + EPSILON); + sqlite3_bind_double(select_stmt_, 3, y - EPSILON); + sqlite3_bind_double(select_stmt_, 4, y + EPSILON); + } else { + // Take into account the shift between source and target horizontal + // coordinates when searching in the reverse direction. + sqlite3_bind_double(select_stmt_, 1, x - max_shift_x_ - EPSILON); + sqlite3_bind_double(select_stmt_, 2, x - min_shift_x_ + EPSILON); + sqlite3_bind_double(select_stmt_, 3, y - max_shift_y_ - EPSILON); + sqlite3_bind_double(select_stmt_, 4, y - min_shift_y_ + EPSILON); + } + + // Iterate over triangles whose bounding box intersecs the point of interest + while (sqlite3_step(select_stmt_) == SQLITE_ROW) { + std::array vX, vY; + if (!getXY(vX, vY)) + return false; + + if (computeLambdas(vX, vY, /* checkWithinTriangle = */ true)) { + const bool status = collectValsVertex(); + if (status) { + cachedForwardDirection_ = forwardDirection; + cachedVerticesX_ = vX; + cachedVerticesY_ = vY; + cachedValsVertex1_ = valsVertex1; + cachedValsVertex2_ = valsVertex2; + cachedValsVertex3_ = valsVertex3; + } + return status; + } + } + + if (fallbackStrategy_ == FALLBACK_NONE) { + return false; + } + + const double x_search = std::clamp( + x + (useForwardDirection ? 0 : -(min_shift_x_ + max_shift_x_) * 0.5), + min_x_, max_x_); + const double y_search = std::clamp( + y + (useForwardDirection ? 0 : -(min_shift_y_ + max_shift_y_) * 0.5), + min_y_, max_y_); + double closest_sq_dist = std::numeric_limits::infinity(); + double closest_dist = std::numeric_limits::infinity(); + std::array closest_vX{0, 0, 0}; + std::array closest_vY{0, 0, 0}; + + // Initiate the search radius from the density of points + // and double it if there is no match. + const double bbox_w = max_x_ - min_x_; + const double bbox_h = max_y_ - min_y_; + double radius = sqrt(bbox_w * bbox_h / numVertices_); + constexpr int MAX_ITERS = 20; + for (int iter = 0; iter <= MAX_ITERS && radius <= std::max(bbox_w, bbox_h); + ++iter) { + sqlite3_reset(select_stmt_); + sqlite3_bind_double(select_stmt_, 1, x_search - radius); + sqlite3_bind_double(select_stmt_, 2, x_search + radius); + sqlite3_bind_double(select_stmt_, 3, y_search - radius); + sqlite3_bind_double(select_stmt_, 4, y_search + radius); + + // Iterate over triangles in the current area of search + while (sqlite3_step(select_stmt_) == SQLITE_ROW) { + std::array vX, vY; + if (!getXY(vX, vY)) + return false; + + // don't check this triangle if the query point plusminus the + // currently closest found distance is outside the triangle's + // bounding box + if (x + closest_dist < std::min(vX[0], std::min(vX[1], vX[2])) || + x - closest_dist > std::max(vX[0], std::max(vX[1], vX[2])) || + y + closest_dist < std::min(vY[0], std::min(vY[1], vY[2])) || + y - closest_dist > std::max(vY[0], std::min(vY[1], vY[2]))) { + continue; + } + const double dist12 = squared_distance(vX[0], vY[0], vX[1], vY[1]); + const double dist23 = squared_distance(vX[1], vY[1], vX[2], vY[2]); + const double dist13 = squared_distance(vX[0], vY[0], vX[2], vY[2]); + if (dist12 < EPSILON || dist23 < EPSILON || dist13 < EPSILON) { + // do not use degenerate triangles + continue; + } + + bool isClosest = false; + if (fallbackStrategy_ == FALLBACK_NEAREST_SIDE) { + // we don't know whether the points of the triangle are given + // clockwise or counter-clockwise, so we have to check the + // distance of the point to all three sides of the triangle + double sq_dist = sq_distance_point_segment( + x, y, vX[0], vY[0], vX[1], vY[1], dist12); + if (sq_dist < closest_sq_dist) { + closest_sq_dist = sq_dist; + isClosest = true; + } + sq_dist = sq_distance_point_segment(x, y, vX[1], vY[1], vX[2], + vY[2], dist23); + if (sq_dist < closest_sq_dist) { + closest_sq_dist = sq_dist; + isClosest = true; + } + sq_dist = sq_distance_point_segment(x, y, vX[0], vY[0], vX[2], + vY[2], dist13); + if (sq_dist < closest_sq_dist) { + closest_sq_dist = sq_dist; + isClosest = true; + } + } else { + assert(fallbackStrategy_ == FALLBACK_NEAREST_CENTROID); + const double c_x = (vX[0] + vX[1] + vX[2]) / 3.0; + const double c_y = (vY[0] + vY[1] + vY[2]) / 3.0; + const double sq_dist = squared_distance(x, y, c_x, c_y); + if (sq_dist < closest_sq_dist) { + closest_sq_dist = sq_dist; + isClosest = true; + } + } + if (isClosest) { + closest_dist = sqrt(closest_sq_dist); + closest_vX = vX; + closest_vY = vY; + valsVertex1.clear(); + valsVertex2.clear(); + valsVertex3.clear(); + if (!collectValsVertex()) + return false; + } + } + + if (std::isinf(closest_sq_dist)) { + // No match: increase search radius + radius *= 2; + } else { + break; + } + } + if (std::isinf(closest_sq_dist)) { + // nothing was found due to empty triangle list or only degenerate + // triangles + return false; + } + + return computeLambdas(closest_vX, closest_vY, + /* checkWithinTriangle = */ false); +} + +// --------------------------------------------------------------------------- + +TINShiftGeopackageEvaluator::TINShiftGeopackageEvaluator( + std::unique_ptr fileIn) + : file_(std::move(fileIn)) {} + +// --------------------------------------------------------------------------- + +bool TINShiftGeopackageEvaluator::transform(bool forwardDirection, double x, + double y, double z, double &x_out, + double &y_out, double &z_out) { + double lambda1 = 0; + double lambda2 = 0; + double lambda3 = 0; + std::vector valsVertex1; + std::vector valsVertex2; + std::vector valsVertex3; + if (!file_->findTriangle(x, y, forwardDirection, lambda1, lambda2, lambda3, + valsVertex1, valsVertex2, valsVertex3)) + return false; + + x_out = x; + y_out = y; + z_out = z; + + if (file_->hasHorizontalTransformation()) { + x_out = lambda1 * valsVertex1[0] + lambda2 * valsVertex2[0] + + lambda3 * valsVertex3[0]; + y_out = lambda1 * valsVertex1[1] + lambda2 * valsVertex2[1] + + lambda3 * valsVertex3[1]; + } + if (file_->hasVerticalTransformation()) { + const int nextAttrIdx = file_->hasHorizontalTransformation() ? 2 : 0; + const double sign = forwardDirection ? 1 : -1; + if (file_->hasOffsetZ()) { + z_out += sign * (lambda1 * valsVertex1[nextAttrIdx] + + lambda2 * valsVertex2[nextAttrIdx] + + lambda3 * valsVertex3[nextAttrIdx]); + } else { + const double offset1 = + valsVertex1[nextAttrIdx + 1] - valsVertex1[nextAttrIdx]; + const double offset2 = + valsVertex2[nextAttrIdx + 1] - valsVertex2[nextAttrIdx]; + const double offset3 = + valsVertex3[nextAttrIdx + 1] - valsVertex3[nextAttrIdx]; + z_out += sign * (lambda1 * offset1 + lambda2 * offset2 + + lambda3 * offset3); + } + } + + return true; +} + +// --------------------------------------------------------------------------- + +bool TINShiftGeopackageEvaluator::forward(double x, double y, double z, + double &x_out, double &y_out, + double &z_out) { + return transform(/* forwardDirection = */ true, x, y, z, x_out, y_out, + z_out); +} + +// --------------------------------------------------------------------------- + +bool TINShiftGeopackageEvaluator::inverse(double x, double y, double z, + double &x_out, double &y_out, + double &z_out) { + return transform(/* forwardDirection = */ false, x, y, z, x_out, y_out, + z_out); +} + +// --------------------------------------------------------------------------- diff --git a/src/transformations/tinshift_gpkg.hpp b/src/transformations/tinshift_gpkg.hpp new file mode 100644 index 0000000000..85d628f904 --- /dev/null +++ b/src/transformations/tinshift_gpkg.hpp @@ -0,0 +1,135 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to TIN based transformations using a + *GeoPackage file Author: Even Rouault, + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault, + * + * SPDX-License-Identifier: MIT + *****************************************************************************/ + +#ifndef TINSHIFT_GEOPACKAGE_H +#define TINSHIFT_GEOPACKAGE_H + +#include +#include +#include +#include +#include + +#include "sqlite3_utils.hpp" +#include "tinshift_iface.hpp" + +struct sqlite3; +struct sqlite3_stmt; + +// --------------------------------------------------------------------------- + +class TINShiftGeopackageException : public std::exception { + public: + explicit TINShiftGeopackageException(const std::string &msg) : msg_(msg) {} + const char *what() const noexcept override; + + private: + std::string msg_; +}; + +// --------------------------------------------------------------------------- + +class TINShiftGeopackageFile { + public: + /** Open the provided GeoPackage and return an object. + * + * @throws TINShiftGeopackageException in case of error. + */ + static std::unique_ptr + open(PJ_CONTEXT *ctx, const std::string &filename); + + /** Destructor */ + ~TINShiftGeopackageFile(); + + /** Find the triangle into which (x, y) is located */ + bool findTriangle(double x, double y, bool forwardDirection, + double &lambda1, double &lambda2, double &lambda3, + std::vector &valsVertex1, + std::vector &valsVertex2, + std::vector &valsVertex3); + + bool hasHorizontalTransformation() const { + return horizontalTransformation_; + } + bool hasVerticalTransformation() const { return verticalTransformation_; } + bool hasOffsetZ() const { return offsetZ_; } + + private: + std::unique_ptr sqlite_vfs_{}; + sqlite3 *sqlite_handle_ = nullptr; + sqlite3_stmt *select_stmt_ = nullptr; + std::vector vertices_cols_{}; + int numVertices_ = 0; + bool horizontalTransformation_ = false; + bool verticalTransformation_ = false; + bool sourceZ_ = false; + bool targetX_ = false; + bool targetY_ = false; + bool targetZ_ = false; + bool offsetZ_ = false; + double min_x_ = 0; + double min_y_ = 0; + double max_x_ = 0; + double max_y_ = 0; + double min_shift_x_ = 0; + double min_shift_y_ = 0; + double max_shift_x_ = 0; + double max_shift_y_ = 0; + + enum FallbackStrategy { + FALLBACK_NONE, + FALLBACK_NEAREST_SIDE, + FALLBACK_NEAREST_CENTROID, + }; + FallbackStrategy fallbackStrategy_ = FALLBACK_NONE; + + // Cached results + bool cachedForwardDirection_ = false; + std::array cachedVerticesX_{0, 0, 0}; + std::array cachedVerticesY_{0, 0, 0}; + std::vector cachedValsVertex1_{}; + std::vector cachedValsVertex2_{}; + std::vector cachedValsVertex3_{}; + + private: + TINShiftGeopackageFile() = default; + TINShiftGeopackageFile(const TINShiftGeopackageFile &) = delete; + TINShiftGeopackageFile &operator=(const TINShiftGeopackageFile &) = delete; + + void getMetadata(); + void prepareDatabaseQuery(); +}; + +// --------------------------------------------------------------------------- + +/** Class to evaluate the transformation of a coordinate */ +class TINShiftGeopackageEvaluator : public TINShiftEvaluator { + public: + /** Constructor. */ + explicit TINShiftGeopackageEvaluator( + std::unique_ptr fileIn); + + bool forward(double x, double y, double z, double &x_out, double &y_out, + double &z_out) override; + + bool inverse(double x, double y, double z, double &x_out, double &y_out, + double &z_out) override; + + private: + std::unique_ptr file_; + + bool transform(bool forwardDirection, double x, double y, double z, + double &x_out, double &y_out, double &z_out); +}; + +// --------------------------------------------------------------------------- + +#endif diff --git a/src/transformations/tinshift_iface.hpp b/src/transformations/tinshift_iface.hpp new file mode 100644 index 0000000000..e2dc75d6c0 --- /dev/null +++ b/src/transformations/tinshift_iface.hpp @@ -0,0 +1,33 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to TIN based transformations + * Author: Even Rouault, + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault, + * + * SPDX-License-Identifier: MIT + *****************************************************************************/ + +#ifndef TINSHIFT_IFACE_H +#define TINSHIFT_IFACE_H + +// --------------------------------------------------------------------------- + +/** Interface to evaluate the transformation of a coordinate */ +class TINShiftEvaluator { + public: + virtual ~TINShiftEvaluator(); + + /** Evaluate displacement of a position given by (x,y,z,t) and + * return it in (x_out,y_out_,z_out). + */ + virtual bool forward(double x, double y, double z, double &x_out, + double &y_out, double &z_out) = 0; + + /** Apply inverse transformation. */ + virtual bool inverse(double x, double y, double z, double &x_out, + double &y_out, double &z_out) = 0; +}; + +#endif // TINSHIFT_IFACE_H diff --git a/src/transformations/tinshift.hpp b/src/transformations/tinshift_json.hpp similarity index 91% rename from src/transformations/tinshift.hpp rename to src/transformations/tinshift_json.hpp index 919b7063af..5beb433574 100644 --- a/src/transformations/tinshift.hpp +++ b/src/transformations/tinshift_json.hpp @@ -25,8 +25,8 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ -#ifndef TINSHIFT_HPP -#define TINSHIFT_HPP +#ifndef TINSHIFT_JSON_HPP +#define TINSHIFT_JSON_HPP #ifdef PROJ_COMPILATION #include "proj/internal/include_nlohmann_json.hpp" @@ -44,14 +44,15 @@ #include #include "quadtree.hpp" +#include "tinshift_iface.hpp" -#ifndef TINSHIFT_NAMESPACE -#define TINSHIFT_NAMESPACE TINShift +#ifndef TINSHIFT_JSON_NAMESPACE +#define TINSHIFT_JSON_NAMESPACE TINShiftJSON #endif -#include "tinshift_exceptions.hpp" +#include "tinshift_json_exceptions.hpp" -namespace TINSHIFT_NAMESPACE { +namespace TINSHIFT_JSON_NAMESPACE { enum FallbackStrategy { FALLBACK_NONE, @@ -64,13 +65,13 @@ using json = nlohmann::json; // --------------------------------------------------------------------------- /** Content of a TINShift file. */ -class TINShiftFile { +class TINShiftJSONFile { public: /** Parse the provided serialized JSON content and return an object. * * @throws ParsingException in case of error. */ - static std::unique_ptr parse(const std::string &text); + static std::unique_ptr parse(const std::string &text); /** Get file type. Should always be "triangulation_file" */ const std::string &fileType() const { return mFileType; } @@ -205,7 +206,7 @@ class TINShiftFile { const std::vector &triangles() const { return mTriangles; } private: - TINShiftFile() = default; + TINShiftJSONFile() = default; std::string mFileType{}; std::string mFormatVersion{}; @@ -229,26 +230,26 @@ class TINShiftFile { // --------------------------------------------------------------------------- /** Class to evaluate the transformation of a coordinate */ -class Evaluator { +class TINShiftJSONEvaluator : public TINShiftEvaluator { public: /** Constructor. */ - explicit Evaluator(std::unique_ptr &&fileIn); + explicit TINShiftJSONEvaluator(std::unique_ptr &&fileIn); /** Get file */ - const TINShiftFile &file() const { return *(mFile.get()); } + const TINShiftJSONFile &file() const { return *(mFile.get()); } /** Evaluate displacement of a position given by (x,y,z,t) and * return it in (x_out,y_out_,z_out). */ bool forward(double x, double y, double z, double &x_out, double &y_out, - double &z_out); + double &z_out) override; /** Apply inverse transformation. */ bool inverse(double x, double y, double z, double &x_out, double &y_out, - double &z_out); + double &z_out) override; private: - std::unique_ptr mFile; + std::unique_ptr mFile; // Reused between invocations to save memory allocations std::vector mTriangleIndices{}; @@ -259,10 +260,10 @@ class Evaluator { // --------------------------------------------------------------------------- -} // namespace TINSHIFT_NAMESPACE +} // namespace TINSHIFT_JSON_NAMESPACE // --------------------------------------------------------------------------- -#include "tinshift_impl.hpp" +#include "tinshift_json_impl.hpp" -#endif // TINSHIFT_HPP +#endif // TINSHIFT_JSON_HPP diff --git a/src/transformations/tinshift_exceptions.hpp b/src/transformations/tinshift_json_exceptions.hpp similarity index 95% rename from src/transformations/tinshift_exceptions.hpp rename to src/transformations/tinshift_json_exceptions.hpp index 1d13c2b160..8b43768226 100644 --- a/src/transformations/tinshift_exceptions.hpp +++ b/src/transformations/tinshift_json_exceptions.hpp @@ -25,13 +25,13 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ -#ifndef TINSHIFT_NAMESPACE +#ifndef TINSHIFT_JSON_NAMESPACE #error "Should be included only by tinshift.hpp" #endif #include -namespace TINSHIFT_NAMESPACE { +namespace TINSHIFT_JSON_NAMESPACE { // --------------------------------------------------------------------------- @@ -49,4 +49,4 @@ const char *ParsingException::what() const noexcept { return msg_.c_str(); } // --------------------------------------------------------------------------- -} // namespace TINSHIFT_NAMESPACE +} // namespace TINSHIFT_JSON_NAMESPACE diff --git a/src/transformations/tinshift_impl.hpp b/src/transformations/tinshift_json_impl.hpp similarity index 96% rename from src/transformations/tinshift_impl.hpp rename to src/transformations/tinshift_json_impl.hpp index 9508ed0a2b..42a192fb6c 100644 --- a/src/transformations/tinshift_impl.hpp +++ b/src/transformations/tinshift_json_impl.hpp @@ -25,7 +25,7 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ -#ifndef TINSHIFT_NAMESPACE +#ifndef TINSHIFT_JSON_NAMESPACE #error "Should be included only by tinshift.hpp" #endif @@ -33,7 +33,7 @@ #include #include -namespace TINSHIFT_NAMESPACE { +namespace TINSHIFT_JSON_NAMESPACE { // --------------------------------------------------------------------------- @@ -76,8 +76,9 @@ static json getArrayMember(const json &j, const char *key) { // --------------------------------------------------------------------------- -std::unique_ptr TINShiftFile::parse(const std::string &text) { - std::unique_ptr tinshiftFile(new TINShiftFile()); +std::unique_ptr +TINShiftJSONFile::parse(const std::string &text) { + std::unique_ptr tinshiftFile(new TINShiftJSONFile()); json j; try { j = json::parse(text); @@ -361,7 +362,7 @@ std::unique_ptr TINShiftFile::parse(const std::string &text) { // --------------------------------------------------------------------------- -static NS_PROJ::QuadTree::RectObj GetBounds(const TINShiftFile &file, +static NS_PROJ::QuadTree::RectObj GetBounds(const TINShiftJSONFile &file, bool forward) { NS_PROJ::QuadTree::RectObj rect; rect.minx = std::numeric_limits::max(); @@ -386,7 +387,7 @@ static NS_PROJ::QuadTree::RectObj GetBounds(const TINShiftFile &file, // --------------------------------------------------------------------------- static std::unique_ptr> -BuildQuadTree(const TINShiftFile &file, bool forward) { +BuildQuadTree(const TINShiftJSONFile &file, bool forward) { auto quadtree = std::unique_ptr>( new NS_PROJ::QuadTree::QuadTree(GetBounds(file, forward))); const auto &triangles = file.triangles(); @@ -425,7 +426,8 @@ BuildQuadTree(const TINShiftFile &file, bool forward) { // --------------------------------------------------------------------------- -Evaluator::Evaluator(std::unique_ptr &&fileIn) +TINShiftJSONEvaluator::TINShiftJSONEvaluator( + std::unique_ptr &&fileIn) : mFile(std::move(fileIn)) {} // --------------------------------------------------------------------------- @@ -452,8 +454,8 @@ static double distance_point_segment(double x, double y, double x1, double y1, return squared_distance(x, y, x1 + t * (x2 - x1), y1 + t * (y2 - y1)); } -static const TINShiftFile::VertexIndices * -FindTriangle(const TINShiftFile &file, +static const TINShiftJSONFile::VertexIndices * +FindTriangle(const TINShiftJSONFile &file, const NS_PROJ::QuadTree::QuadTree &quadtree, std::vector &triangleIndices, double x, double y, bool forward, double &lambda1, double &lambda2, double &lambda3) { @@ -594,8 +596,8 @@ FindTriangle(const TINShiftFile &file, // --------------------------------------------------------------------------- -bool Evaluator::forward(double x, double y, double z, double &x_out, - double &y_out, double &z_out) { +bool TINShiftJSONEvaluator::forward(double x, double y, double z, double &x_out, + double &y_out, double &z_out) { if (!mQuadTreeForward) mQuadTreeForward = BuildQuadTree(*(mFile.get()), true); @@ -638,8 +640,8 @@ bool Evaluator::forward(double x, double y, double z, double &x_out, // --------------------------------------------------------------------------- -bool Evaluator::inverse(double x, double y, double z, double &x_out, - double &y_out, double &z_out) { +bool TINShiftJSONEvaluator::inverse(double x, double y, double z, double &x_out, + double &y_out, double &z_out) { NS_PROJ::QuadTree::QuadTree *quadtree; if (!mFile->transformHorizontalComponent() && mFile->transformVerticalComponent()) { @@ -688,4 +690,4 @@ bool Evaluator::inverse(double x, double y, double z, double &x_out, return true; } -} // namespace TINSHIFT_NAMESPACE +} // namespace TINSHIFT_JSON_NAMESPACE diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 19de7df692..86c2ff16d7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -78,8 +78,13 @@ proj_add_gie_test("adams_ws2" "gie/adams_ws2.gie") proj_add_gie_test("guyou" "gie/guyou.gie") proj_add_gie_test("peirce_q" "gie/peirce_q.gie") proj_add_gie_test("tinshift" "gie/tinshift.gie") +proj_add_gie_test("tinshift_gpkg" "gie/tinshift_gpkg.gie") proj_add_gie_test("spilhaus" "gie/spilhaus.gie") +if(CURL_ENABLED AND RUN_NETWORK_DEPENDENT_TESTS) +proj_add_gie_network_dependent_test("tinshift_gpkg_network" "gie/tinshift_gpkg_network.gie") +endif() + if(TIFF_ENABLED) proj_add_gie_test("Deformation" "gie/deformation.gie") proj_add_gie_test("geotiff_grids" "gie/geotiff_grids.gie") diff --git a/test/gie/tinshift_gpkg.gie b/test/gie/tinshift_gpkg.gie new file mode 100644 index 0000000000..978e62b648 --- /dev/null +++ b/test/gie/tinshift_gpkg.gie @@ -0,0 +1,58 @@ + +------------------------------------------------------------------------------- +=============================================================================== +Test +proj=tinshift with GeoPackage files +=============================================================================== + + + +# +file doesn't point to an existing file +operation +proj=tinshift +file=i_do_not_exist.gpkg +expect failure errno invalid_op_file_not_found_or_invalid + +# Not a GPKG file +operation +proj=tinshift +file=tests/tinshift_empty_file.gpkg +expect failure errno invalid_op_file_not_found_or_invalid + + +# Tests on a file without explicit CRS +operation +proj=tinshift +file=tests/tinshift_crs_implicit.gpkg +accept 2 49 +expect 2.1 49.1 +roundtrip 10 + +accept 0 0 +expect failure + +direction inverse +accept 0 0 +expect failure + + +# Tests on a file with explicit CRS +operation +proj=tinshift +file=tests/tinshift_simplified_kkj_etrs.gpkg +tolerance 0.1 mm +# Verified with https://kartta.paikkatietoikkuna.fi/?lang=en with EPSG:2393 to EPSG:3067 +accept 3210000.0000 6700000.0000 +expect 209948.3217 6697187.0009 +roundtrip 10 + +operation +proj=tinshift +file=tests/tinshift_simplified_n60_n2000.gpkg +tolerance 0.1 mm +accept 3210000.0000 6700000.0000 10.0 +expect 3210000.0000 6700000.0000 10.2886 +roundtrip 10 + +# Test fallback strategy nearest_side +operation +proj=tinshift +file=tests/tinshift_fallback_nearest_side.gpkg +accept 2 3 +expect 4 6 +roundtrip 1 + +# Test fallback strategy nearest_centroid +operation +proj=tinshift +file=tests/tinshift_fallback_nearest_centroid.gpkg +accept 3 0 +expect 3 0 +roundtrip 1 + + diff --git a/test/gie/tinshift_gpkg_network.gie b/test/gie/tinshift_gpkg_network.gie new file mode 100644 index 0000000000..6c0cd93996 --- /dev/null +++ b/test/gie/tinshift_gpkg_network.gie @@ -0,0 +1,27 @@ + +------------------------------------------------------------------------------- +=============================================================================== +Test +proj=tinshift with GeoPackage files on network +=============================================================================== + + + +# +file doesn't point to an existing file +operation +proj=tinshift +file=https://cdn.proj.org/i_do_not_exist.gpkg +expect failure errno invalid_op_file_not_found_or_invalid + +# Tests on a file with explicit CRS +operation +proj=tinshift +file=https://cdn.proj.org/fi_nls_ykj_etrs35fin.gpkg +tolerance 0.1 mm +# Verified with https://kartta.paikkatietoikkuna.fi/?lang=en with EPSG:2393 to EPSG:3067 +accept 3210000.0000 6700000.0000 +expect 209948.3217 6697187.0009 +roundtrip 1 + +operation +proj=tinshift +file=fi_nls_n60_n2000.gpkg +tolerance 0.1 mm +accept 3210000.0000 6700000.0000 10.0 +expect 3210000.0000 6700000.0000 10.2886 +roundtrip 1 + + diff --git a/test/unit/test_network.cpp b/test/unit/test_network.cpp index ecced83351..d4802b9dfa 100644 --- a/test/unit/test_network.cpp +++ b/test/unit/test_network.cpp @@ -1897,14 +1897,16 @@ TEST(networking, download_whole_files) { char out_full_filename[1024]; EXPECT_FALSE(pj_find_file(ctx, "dk_sdfe_dvr90.tif", out_full_filename, - sizeof(out_full_filename))); + sizeof(out_full_filename), + /* disable_network = */ true)); EXPECT_STREQ(out_full_filename, ""); ASSERT_TRUE( proj_download_file(ctx, "dk_sdfe_dvr90.tif", false, nullptr, nullptr)); EXPECT_TRUE(pj_find_file(ctx, "dk_sdfe_dvr90.tif", out_full_filename, - sizeof(out_full_filename))); + sizeof(out_full_filename), + /* disable_network = */ true)); EXPECT_NE(out_full_filename[0], 0); // lookForGridInfo() returns false because the grid is not known in the DB, diff --git a/test/unit/test_tinshift.cpp b/test/unit/test_tinshift.cpp index 2b42294fda..86b8d3840c 100644 --- a/test/unit/test_tinshift.cpp +++ b/test/unit/test_tinshift.cpp @@ -29,10 +29,12 @@ #include "gtest_include.h" #define PROJ_COMPILATION -#define TINSHIFT_NAMESPACE TestTINShift -#include "transformations/tinshift.hpp" +#define TINSHIFT_JSON_NAMESPACE TestTINShift +#include "transformations/tinshift_json.hpp" -using namespace TINSHIFT_NAMESPACE; +using namespace TINSHIFT_JSON_NAMESPACE; + +TINShiftEvaluator::~TINShiftEvaluator() = default; namespace { @@ -54,20 +56,20 @@ static json getMinValidContent() { // --------------------------------------------------------------------------- TEST(tinshift, basic) { - EXPECT_THROW(TINShiftFile::parse("foo"), ParsingException); - EXPECT_THROW(TINShiftFile::parse("null"), ParsingException); - EXPECT_THROW(TINShiftFile::parse("{}"), ParsingException); + EXPECT_THROW(TINShiftJSONFile::parse("foo"), ParsingException); + EXPECT_THROW(TINShiftJSONFile::parse("null"), ParsingException); + EXPECT_THROW(TINShiftJSONFile::parse("{}"), ParsingException); const auto jMinValid(getMinValidContent()); { - auto f = TINShiftFile::parse(jMinValid.dump()); + auto f = TINShiftJSONFile::parse(jMinValid.dump()); EXPECT_EQ(f->fileType(), "triangulation_file"); EXPECT_EQ(f->formatVersion(), "1.0"); EXPECT_EQ(f->inputCRS(), "EPSG:2393"); EXPECT_EQ(f->outputCRS(), "EPSG:3067"); EXPECT_EQ(f->fallbackStrategy(), FALLBACK_NONE); - auto eval = Evaluator(std::move(f)); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0; @@ -121,8 +123,8 @@ TEST(tinshift, basic) { {0, 0, 10.5, 10.6}, {0, 1, 15.0, 15.2}, {1, 1, 17.5, 18.0}}; j["triangles"] = {{0, 1, 2}}; - auto f = TINShiftFile::parse(j.dump()); - auto eval = Evaluator(std::move(f)); + auto f = TINShiftJSONFile::parse(j.dump()); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0; @@ -152,8 +154,8 @@ TEST(tinshift, basic) { j["vertices"] = {{0, 0, 0.1}, {0, 1, 0.2}, {1, 1, 0.5}}; j["triangles"] = {{0, 1, 2}}; - auto f = TINShiftFile::parse(j.dump()); - auto eval = Evaluator(std::move(f)); + auto f = TINShiftJSONFile::parse(j.dump()); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0; @@ -186,8 +188,8 @@ TEST(tinshift, basic) { {1, 1, 100, 100, 0.5}}; j["triangles"] = {{0, 1, 2}}; - auto f = TINShiftFile::parse(j.dump()); - auto eval = Evaluator(std::move(f)); + auto f = TINShiftJSONFile::parse(j.dump()); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0; @@ -212,7 +214,7 @@ TEST(tinshift, basic) { { auto j(jMinValid); j["fallback_strategy"] = "none"; - EXPECT_THROW(TINShiftFile::parse(j.dump()), ParsingException); + EXPECT_THROW(TINShiftJSONFile::parse(j.dump()), ParsingException); } // invalid fallback_strategy field with 1.1 version @@ -220,7 +222,7 @@ TEST(tinshift, basic) { auto j(jMinValid); j["format_version"] = "1.1"; j["fallback_strategy"] = "invalid"; - EXPECT_THROW(TINShiftFile::parse(j.dump()), ParsingException); + EXPECT_THROW(TINShiftJSONFile::parse(j.dump()), ParsingException); } // fail with no triangles and fallback nearest_side @@ -230,8 +232,8 @@ TEST(tinshift, basic) { j["fallback_strategy"] = "nearest_side"; j["triangles"] = json::array(); // empty - auto f = TINShiftFile::parse(j.dump()); - auto eval = Evaluator(std::move(f)); + auto f = TINShiftJSONFile::parse(j.dump()); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0; @@ -247,8 +249,8 @@ TEST(tinshift, basic) { j["fallback_strategy"] = "nearest_side"; j["vertices"] = {{0, 0, 101, 101}, {0, 1, 100, 101}, {0, 1, 100, 100}}; - auto f = TINShiftFile::parse(j.dump()); - auto eval = Evaluator(std::move(f)); + auto f = TINShiftJSONFile::parse(j.dump()); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0; @@ -265,8 +267,8 @@ TEST(tinshift, basic) { j["vertices"] = { {0, 0, 101, 101}, {0, 0.5, 100, 101}, {0, 1, 100, 100}}; - auto f = TINShiftFile::parse(j.dump()); - auto eval = Evaluator(std::move(f)); + auto f = TINShiftJSONFile::parse(j.dump()); + auto eval = TINShiftJSONEvaluator(std::move(f)); double x_out = 0; double y_out = 0; double z_out = 0;