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 5db0fe5b5e..df846e6631 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 95aebfcff6..9fcab009cf 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 "ac049b1ce88561ff3a47f87c922df3a2") +set(PROJ_DB_SQL_EXPECTED_MD5 "e559aa4bb2ae349dd527dc64b43f28d2") 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 0000000000..5285995887 Binary files /dev/null and b/data/tests/tinshift_crs_implicit.gpkg differ 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 0000000000..7f2069cdb3 Binary files /dev/null and b/data/tests/tinshift_fallback_nearest_centroid.gpkg differ diff --git a/data/tests/tinshift_fallback_nearest_side.gpkg b/data/tests/tinshift_fallback_nearest_side.gpkg new file mode 100644 index 0000000000..cc55230ac1 Binary files /dev/null and b/data/tests/tinshift_fallback_nearest_side.gpkg differ diff --git a/data/tests/tinshift_simplified_kkj_etrs.gpkg b/data/tests/tinshift_simplified_kkj_etrs.gpkg new file mode 100644 index 0000000000..aff12442a6 Binary files /dev/null and b/data/tests/tinshift_simplified_kkj_etrs.gpkg differ diff --git a/data/tests/tinshift_simplified_n60_n2000.gpkg b/data/tests/tinshift_simplified_n60_n2000.gpkg new file mode 100644 index 0000000000..dd48e1de26 Binary files /dev/null and b/data/tests/tinshift_simplified_n60_n2000.gpkg differ diff --git a/docs/source/install.rst b/docs/source/install.rst index 73eb5a40f0..76b9f4c2f9 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -460,7 +460,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. @@ -618,3 +618,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 30d191d8ef..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,100 +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) -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 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 1853ca53a5..e39b734dde 100644 --- a/docs/source/specifications/index.rst +++ b/docs/source/specifications/index.rst @@ -12,3 +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 new file mode 100644 index 0000000000..910421f7c8 --- /dev/null +++ b/docs/source/specifications/tin_json.rst @@ -0,0 +1,107 @@ +.. _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 + +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: + +.. 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. 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 8361a85003..3fb7b70ea8 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 feaf2c4b9b..16ca0bc283 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -12446,7 +12446,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 771d4a1328..c3e8cb15cc 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;