diff options
44 files changed, 5552 insertions, 200 deletions
diff --git a/data/Makefile.am b/data/Makefile.am index c4848de1..b48f53fa 100644 --- a/data/Makefile.am +++ b/data/Makefile.am @@ -3,7 +3,8 @@ DATAPATH = $(top_srcdir)/data pkgdata_DATA = proj.ini GL27 nad.lst nad27 nad83 world other.extra \ CH \ ITRF2000 ITRF2008 ITRF2014 proj.db \ - projjson.schema.json + projjson.schema.json \ + deformation_model.schema.json SQL_ORDERED_LIST = sql/begin.sql \ sql/proj_db_table_defs.sql \ @@ -43,6 +44,7 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \ CH \ ITRF2000 ITRF2008 ITRF2014 \ projjson.schema.json \ + deformation_model.schema.json \ CMakeLists.txt \ tests/test_nodata.gtx \ tests/test_vgrid_bigendian_bigtiff.tif \ @@ -89,6 +91,7 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \ tests/us_noaa_geoid06_ak_subset_at_antimeridian.tif \ tests/test_hgrid_little_endian.gsb \ tests/test_hgrid_big_endian.gsb \ + tests/test_3d_grid_projected.tif \ tests/BETA2007.gsb \ tests/MD \ tests/alaska \ @@ -97,6 +100,22 @@ EXTRA_DIST = proj.ini GL27 nad.lst nad27 nad83 \ tests/ntv1_can.dat \ tests/ntv2_0_downsampled.gsb \ tests/ntf_r93.gsb \ + tests/simple_model_degree_3d_grid.tif \ + tests/simple_model_degree_horizontal.json \ + tests/simple_model_degree_3d.json \ + tests/simple_model_metre_3d_grid.tif \ + tests/simple_model_metre_horizontal.json \ + tests/simple_model_metre_3d.json \ + tests/simple_model_metre_3d_geocentric.json \ + tests/simple_model_metre_vertical_grid.tif \ + tests/simple_model_metre_vertical.json \ + tests/simple_model_polar.json \ + tests/simple_model_polar.tif \ + tests/simple_model_wrap_east.json \ + tests/simple_model_wrap_east.tif \ + tests/simple_model_wrap_west.json \ + tests/simple_model_wrap_west.tif \ + tests/simple_model_projected.json \ generate_all_sql_in.cmake sql_filelist.cmake \ $(SQL_ORDERED_LIST) diff --git a/data/deformation_model.schema.json b/data/deformation_model.schema.json new file mode 100644 index 00000000..d7a6d162 --- /dev/null +++ b/data/deformation_model.schema.json @@ -0,0 +1,582 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "description": "Schema for deformation models", + "type": "object", + "properties": { + "file_type": { + "type": "string", + "enum": [ + "deformation_model_master_file" + ], + "description": "File type. Always \"deformation_model_master_file\"" + }, + "format_version": { + "type": "string", + "enum": [ + "1.0" + ] + }, + "name": { + "type": "string", + "description": "A brief descriptive name of the deformation model" + }, + "version": { + "type": "string", + "description": "A string identifying the version of the deformation model. The format for specifying version will be defined by the agency responsible for the deformation model" + }, + "publication_date": { + "$ref": "#/definitions/datetime", + "description": "The date on which this version of the deformation model was published (or possibly the date on which it takes effect?)" + }, + "license": { + "type": "string", + "description": "License under which the model is published" + }, + "description": { + "type": "string", + "description": "A text description of the model" + }, + "authority": { + "type": "object", + "description": "Basic information about the agency responsible for the data set", + "properties": { + "name": { + "type": "string", + "description": "The name of the agency" + }, + "url": { + "type": "string", + "description": "The url of the agency website", + "format": "uri" + }, + "address": { + "type": "string", + "description": "The postal address of the agency" + }, + "email": { + "type": "string", + "description": "An email contact address for the agency", + "format": "email" + } + }, + "required": [ + "name" + ], + "additionalProperties": false + }, + "links": { + "type": "array", + "description": "Links to related information", + "items": { + "type": "object", + "properties": { + "href": { + "type": "string", + "description": "The URL holding the information", + "format": "uri" + }, + "rel": { + "type": "string", + "description": "The relationship to the dataset. Proposed relationships are:\n- \"about\": a web page for human consumption describing the model\n- \"source\": the authoritative source data from which the deformation model is built.\n- \"metadata\": ISO 19115 XML metadata regarding the deformation model." + }, + "type": { + "type": "string", + "description": "MIME type" + }, + "title": { + "type": "string", + "description": "Description of the link" + } + }, + "required": [ + "href" + ], + "additionalProperties": false + } + }, + "source_crs": { + "$ref": "#/definitions/crs", + "description": "The coordinate reference system to which the deformation model applies" + }, + "target_crs": { + "$ref": "#/definitions/crs", + "description": "For a time dependent coordinate transformation the coordinate reference system resulting from applying the deformation" + }, + "definition_crs": { + "$ref": "#/definitions/crs", + "description": "The coordinate reference system used to define the component spatial models. This proposal only supports using the same value for the source and definition coordinate reference system." + }, + "reference_epoch": { + "$ref": "#/definitions/datetime", + "description": "A nominal reference epoch of the deformation model. This is not necessarily used to calculate the deformation model - each component defines its own time function." + }, + "uncertainty_reference_epoch": { + "$ref": "#/definitions/datetime", + "description": "The uncertainties of the deformation model are calculated in terms of this epoch. This is described below in the Time functions section." + }, + "horizontal_offset_unit": { + "type": "string", + "enum": [ + "metre", + "degree" + ] + }, + "vertical_offset_unit": { + "type": "string", + "enum": [ + "metre" + ] + }, + "horizontal_uncertainty_type": { + "type": "string", + "enum": [ + "circular 95% confidence limit" + ] + }, + "horizontal_uncertainty_unit": { + "type": "string", + "enum": [ + "metre" + ] + }, + "vertical_uncertainty_type": { + "type": "string", + "enum": [ + "95% confidence limit" + ] + }, + "vertical_uncertainty_unit": { + "type": "string", + "enum": [ + "metre" + ] + }, + "horizontal_offset_method": { + "type": "string", + "description": "Defines how the horizontal offsets are applied to geographic coordinates", + "enum": [ + "addition", + "geocentric" + ] + }, + "extent": { + "$ref": "#/definitions/extent", + "description": "Defines the region within which the deformation model is defined. It cannot be calculated outside this region. The region is specified by a type and value. This proposal only supports using a bounding box as an array of [west,south,east,north] coordinate values" + }, + "time_extent": { + "type": "object", + "description": "Defines the range of times for which the model is valid, specified by a first and a last value. The deformation model is undefined for dates outside this range.", + "properties": { + "first": { + "$ref": "#/definitions/datetime" + }, + "last": { + "$ref": "#/definitions/datetime" + } + }, + "required": [ + "first", + "last" + ], + "additionalProperties": false + }, + "components": { + "type": "array", + "items": { + "$ref": "#/definitions/component" + } + } + }, + "required": [ + "file_type", + "format_version", + "source_crs", + "target_crs", + "definition_crs", + "extent", + "time_extent", + "components" + ], + "additionalProperties": false, + "definitions": { + "component": { + "type": "object", + "definition": "A component describes an aspect of the deformation, such as glacial isostatic adjustment, secular deformation, earthquakes, etc.", + "properties": { + "description": { + "type": "string", + "description": "A text description of this component of the model" + }, + "extent": { + "$ref": "#/definitions/extent", + "description": "The region within the component is defined. Outside this region the component evaluates to 0. The region is specified by a type and value. This proposal only supports using a bounding box as an array of [west,south,east,north] coordinate values" + }, + "displacement_type": { + "type": "string", + "description": "The displacement parameters defined by the model. The \"none\" option allows for a component which defines uncertainty with different grids to those defining displacement", + "enum": [ + "none", + "horizontal", + "vertical", + "3d" + ] + }, + "uncertainty_type": { + "type": "string", + "description": "The uncertainty parameters defined by the model", + "enum": [ + "none", + "horizontal", + "vertical", + "3d" + ] + }, + "horizontal_uncertainty": { + "type": "number", + "description": "The horizontal uncertainty to use if it is not defined explicitly in the spatial model" + }, + "vertical_uncertainty": { + "type": "number", + "description": "The vertical uncertainty to use if it is not defined explicitly in the spatial model" + }, + "spatial_model": { + "type": "object", + "description": "Defines the spatial model", + "properties": { + "type": { + "type": "string", + "description": "Specifies the type of the spatial model data file. Initially it is proposed that only GeoTIFF is supported", + "enum": [ + "GeoTIFF" + ] + }, + "interpolation_method": { + "type": "string", + "description": "Interpolation method", + "enum": [ + "bilinear", + "geocentric_bilinear" + ] + }, + "filename": { + "type": "string", + "description": "Specifies location of the spatial model GeoTIFF file relative to this JSON file" + }, + "md5_checksum": { + "type": "string", + "description": "A hex encoded MD5 checksum of the grid file that can be used to validate that it is the correct version of the file" + } + }, + "required": [ + "type", + "interpolation_method", + "filename" + ], + "additionalProperties": false + }, + "time_function": { + "$ref": "#/definitions/time_function" + } + }, + "required": [ + "description", + "extent", + "displacement_type", + "spatial_model", + "time_function" + ], + "additionalProperties": false + }, + "crs": { + "type": "string", + "pattern": "^[a-zA-Z]+:[a-zA-Z0-9]+$" + }, + "datetime": { + "type": "string", + "format": "date-time", + "pattern": "^[0-9]{4}-[0-9]{2}-[0-9]{2}T[0-9]{2}:[0-9]{2}:[0-9]{2}Z$" + }, + "extent": { + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "bbox" + ] + }, + "parameters": { + "type": "object", + "properties": { + "bbox": { + "type": "array", + "minItems": 4, + "maxItems": 4, + "items": { + "type": "number" + } + } + } + } + }, + "required": [ + "type", + "parameters" + ], + "additionalProperties": false + }, + "time_function": { + "description": "Function describing a multiplicative factor to apply to the spatial_model depending on the time", + "oneOf": [ + { + "$ref": "#/definitions/time_function_constant" + }, + { + "$ref": "#/definitions/time_function_velocity" + }, + { + "$ref": "#/definitions/time_function_step" + }, + { + "$ref": "#/definitions/time_function_reverse_step" + }, + { + "$ref": "#/definitions/time_function_piecewise" + }, + { + "$ref": "#/definitions/time_function_exponential" + } + ] + }, + "time_function_constant": { + "description": "The valuation of this function is 1 at any epoch", + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "constant" + ] + }, + "parameters": { + "type": "object", + "properties": { + }, + "additionalProperties": false + } + }, + "required": [ + "type" + ], + "additionalProperties": false + }, + "time_function_velocity": { + "description": "The valuation of this function is 0 at reference_epoch, and proportional to the time difference to it at other times", + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "velocity" + ] + }, + "parameters": { + "type": "object", + "properties": { + "reference_epoch": { + "$ref": "#/definitions/datetime" + } + }, + "required": [ + "reference_epoch" + ], + "additionalProperties": false + } + }, + "required": [ + "type", + "parameters" + ], + "additionalProperties": false + }, + "time_function_step": { + "description": "The valuation of this function is 0 before step_epoch, and 1 starting from it", + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "step" + ] + }, + "parameters": { + "type": "object", + "properties": { + "step_epoch": { + "$ref": "#/definitions/datetime" + } + }, + "required": [ + "step_epoch" + ], + "additionalProperties": false + } + }, + "required": [ + "type", + "parameters" + ], + "additionalProperties": false + }, + "time_function_reverse_step": { + "description": "The valuation of this function is 1 before step_epoch, and 0 starting from it", + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "reverse_step" + ] + }, + "parameters": { + "type": "object", + "properties": { + "step_epoch": { + "$ref": "#/definitions/datetime" + } + }, + "required": [ + "step_epoch" + ], + "additionalProperties": false + } + }, + "required": [ + "type", + "parameters" + ], + "additionalProperties": false + }, + "time_function_piecewise": { + "description": "Piecewise time function", + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "piecewise" + ] + }, + "parameters": { + "type": "object", + "properties": { + "before_first": { + "type": "string", + "description": "Defines the behaviour of the function before the first defined epoch", + "enum": [ + "zero", + "constant", + "linear" + ] + }, + "after_last": { + "type": "string", + "description": "Defines the behaviour of the function after the last defined epoch", + "enum": [ + "zero", + "constant", + "linear" + ] + }, + "model": { + "type": "array", + "description": "A sorted array data points each defined by two elements, \"epoch\" defines the date/time of the data point, and \"scale_factor\" is the corresponding function value. The array is sorted in order of increasing epoch. Note: where the time function includes a step it is represented by two consecutive data points with the same epoch. The first defines the scale factor that applies before the epoch and the second the scale factor that applies after the epoch", + "items": { + "type": "object", + "properties": { + "epoch": { + "$ref": "#/definitions/datetime" + }, + "scale_factor": { + "type": "number" + } + }, + "required": [ + "epoch", + "scale_factor" + ], + "additionalProperties": false + }, + "minItems": 2 + } + }, + "required": [ + "before_first", + "after_last", + "model" + ], + "additionalProperties": false + } + }, + "required": [ + "type", + "parameters" + ], + "additionalProperties": false + }, + "time_function_exponential": { + "description": "The valuation of this function is an exponential function with a time-based relaxation constant", + "type": "object", + "properties": { + "type": { + "type": "string", + "enum": [ + "exponential" + ] + }, + "parameters": { + "type": "object", + "properties": { + "reference_epoch": { + "$ref": "#/definitions/datetime", + "description": "The date/time at which the exponential decay starts" + }, + "end_epoch": { + "$ref": "#/definitions/datetime", + "description": "The date/time at which the exponential decay ends (optional)" + }, + "relaxation_constant": { + "type": "number", + "description": "Relaxation constant in years" + }, + "before_scale_factor": { + "type": "number", + "description": "The scale factor that applies before the reference epoch" + }, + "initial_scale_factor": { + "type": "number", + "description": "The initial scale factor" + }, + "final_scale_factor": { + "type": "number", + "description": "The scale factor the exponential function approaches" + } + }, + "required": [ + "reference_epoch", + "relaxation_constant", + "before_scale_factor", + "initial_scale_factor", + "final_scale_factor" + ], + "additionalProperties": false + } + }, + "required": [ + "type", + "parameters" + ], + "additionalProperties": false + } + } +}
\ No newline at end of file diff --git a/data/tests/simple_model_degree_3d.json b/data/tests/simple_model_degree_3d.json new file mode 100644 index 00000000..6a41ce5c --- /dev/null +++ b/data/tests/simple_model_degree_3d.json @@ -0,0 +1,54 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "horizontal_offset_unit": "degree", + "horizontal_offset_method": "addition", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "3d", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/simple_model_degree_3d_grid.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +}
\ No newline at end of file diff --git a/data/tests/simple_model_degree_3d_grid.tif b/data/tests/simple_model_degree_3d_grid.tif Binary files differnew file mode 100644 index 00000000..3bbff0a6 --- /dev/null +++ b/data/tests/simple_model_degree_3d_grid.tif diff --git a/data/tests/simple_model_degree_horizontal.json b/data/tests/simple_model_degree_horizontal.json new file mode 100644 index 00000000..533c0054 --- /dev/null +++ b/data/tests/simple_model_degree_horizontal.json @@ -0,0 +1,53 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "horizontal_offset_unit": "degree", + "horizontal_offset_method": "addition", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "horizontal", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/simple_model_degree_3d_grid.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +}
\ No newline at end of file diff --git a/data/tests/simple_model_metre_3d.json b/data/tests/simple_model_metre_3d.json new file mode 100644 index 00000000..201aaebb --- /dev/null +++ b/data/tests/simple_model_metre_3d.json @@ -0,0 +1,54 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "horizontal_offset_unit": "metre", + "horizontal_offset_method": "addition", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "3d", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/simple_model_metre_3d_grid.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +}
\ No newline at end of file diff --git a/data/tests/simple_model_metre_3d_geocentric.json b/data/tests/simple_model_metre_3d_geocentric.json new file mode 100644 index 00000000..1328ff5e --- /dev/null +++ b/data/tests/simple_model_metre_3d_geocentric.json @@ -0,0 +1,54 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "horizontal_offset_unit": "metre", + "horizontal_offset_method": "geocentric", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "3d", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/simple_model_metre_3d_grid.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +}
\ No newline at end of file diff --git a/data/tests/simple_model_metre_3d_grid.tif b/data/tests/simple_model_metre_3d_grid.tif Binary files differnew file mode 100644 index 00000000..40cf2d70 --- /dev/null +++ b/data/tests/simple_model_metre_3d_grid.tif diff --git a/data/tests/simple_model_metre_horizontal.json b/data/tests/simple_model_metre_horizontal.json new file mode 100644 index 00000000..d0ac477c --- /dev/null +++ b/data/tests/simple_model_metre_horizontal.json @@ -0,0 +1,53 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "horizontal_offset_unit": "metre", + "horizontal_offset_method": "addition", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "horizontal", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/simple_model_metre_3d_grid.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +}
\ No newline at end of file diff --git a/data/tests/simple_model_metre_vertical.json b/data/tests/simple_model_metre_vertical.json new file mode 100644 index 00000000..70574340 --- /dev/null +++ b/data/tests/simple_model_metre_vertical.json @@ -0,0 +1,49 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "vertical", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + -180, + -90, + 180, + 90 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/simple_model_metre_vertical_grid.tif" + }, + "time_function": { + "type": "constant" + } + } + ] +}
\ No newline at end of file diff --git a/data/tests/simple_model_metre_vertical_grid.tif b/data/tests/simple_model_metre_vertical_grid.tif Binary files differnew file mode 100644 index 00000000..058b8081 --- /dev/null +++ b/data/tests/simple_model_metre_vertical_grid.tif diff --git a/data/tests/simple_model_polar.json b/data/tests/simple_model_polar.json new file mode 100644 index 00000000..ef99a0cb --- /dev/null +++ b/data/tests/simple_model_polar.json @@ -0,0 +1,54 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "horizontal_offset_unit": "metre", + "horizontal_offset_method": "geocentric", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + 0, + -90, + 360, + -89 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "3d", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + 0, + -90, + 360, + -89 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "geocentric_bilinear", + "filename": "tests/simple_model_polar.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +} diff --git a/data/tests/simple_model_polar.tif b/data/tests/simple_model_polar.tif Binary files differnew file mode 100644 index 00000000..7371ca1e --- /dev/null +++ b/data/tests/simple_model_polar.tif diff --git a/data/tests/simple_model_projected.json b/data/tests/simple_model_projected.json new file mode 100644 index 00000000..c97a7c11 --- /dev/null +++ b/data/tests/simple_model_projected.json @@ -0,0 +1,51 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:2193", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:2193", + "horizontal_offset_unit": "metre", + "horizontal_offset_method": "addition", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + 1500000.0, + 5400000.0, + 1501000.0, + 5401000.0 + ] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "3d", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [ + 1500000.0, + 5400000.0, + 1501000.0, + 5401000.0 + ] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "bilinear", + "filename": "tests/test_3d_grid_projected.tif" + }, + "time_function": { + "type": "constant" + } + } + ] +} diff --git a/data/tests/simple_model_wrap_east.json b/data/tests/simple_model_wrap_east.json new file mode 100644 index 00000000..2a0a2c0e --- /dev/null +++ b/data/tests/simple_model_wrap_east.json @@ -0,0 +1,42 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [-194.2, -37.5, -193.8, -37.2] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "vertical", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [-194.2, -37.5, -193.8, -37.2] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "geocentric_bilinear", + "filename": "tests/simple_model_wrap_east.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +} diff --git a/data/tests/simple_model_wrap_east.tif b/data/tests/simple_model_wrap_east.tif Binary files differnew file mode 100644 index 00000000..816d8a7a --- /dev/null +++ b/data/tests/simple_model_wrap_east.tif diff --git a/data/tests/simple_model_wrap_west.json b/data/tests/simple_model_wrap_west.json new file mode 100644 index 00000000..54a04bd2 --- /dev/null +++ b/data/tests/simple_model_wrap_west.json @@ -0,0 +1,42 @@ +{ + "file_type": "deformation_model_master_file", + "format_version": "1.0", + "source_crs": "EPSG:4326", + "target_crs": "foo:ignored", + "definition_crs": "EPSG:4326", + "vertical_offset_unit": "metre", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [525.8,-37.5,526.2,-37.2] + } + }, + "time_extent": { + "first": "1900-01-01T00:00:00Z", + "last": "2050-01-01T00:00:00Z" + }, + "components": [ + { + "description": "test", + "displacement_type": "vertical", + "uncertainty_type": "none", + "extent": { + "type": "bbox", + "parameters": { + "bbox": [525.8,-37.5,526.2,-37.2] + } + }, + "spatial_model": { + "type": "GeoTIFF", + "interpolation_method": "geocentric_bilinear", + "filename": "tests/simple_model_wrap_west.tif" + }, + "time_function": { + "type": "step", + "parameters": { + "step_epoch": "1900-01-01T00:00:00Z" + } + } + } + ] +} diff --git a/data/tests/simple_model_wrap_west.tif b/data/tests/simple_model_wrap_west.tif Binary files differnew file mode 100644 index 00000000..3a8da6f6 --- /dev/null +++ b/data/tests/simple_model_wrap_west.tif diff --git a/data/tests/test_3d_grid_projected.tif b/data/tests/test_3d_grid_projected.tif Binary files differnew file mode 100644 index 00000000..56138417 --- /dev/null +++ b/data/tests/test_3d_grid_projected.tif diff --git a/docs/source/operations/transformations/defmodel.rst b/docs/source/operations/transformations/defmodel.rst new file mode 100644 index 00000000..fe2ac66a --- /dev/null +++ b/docs/source/operations/transformations/defmodel.rst @@ -0,0 +1,61 @@ +.. _defmodel: + +================================================================================ +Multi-component time-based deformation model +================================================================================ + +.. versionadded:: 7.1.0 + ++---------------------+--------------------------------------------------------------------+ +| **Input type** | Geodetic or projected coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++---------------------+--------------------------------------------------------------------+ +| **Output type** | Geodetic or projected coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++---------------------+--------------------------------------------------------------------+ +| **Domain** | 4D | ++---------------------+--------------------------------------------------------------------+ +| **Available forms** | Forward and inverse | ++---------------------+--------------------------------------------------------------------+ + +The defmodel transformation can be used to represent most deformation models +currently in use, in particular for areas subject to complex deformation, including +large scale secular crustal deformation near plate boundaries and vertical deformation +due to Glacial Isostatic Adjustment (GIA). These can often be represented by a constant +velocity model. Additionally, many areas suffer episodic deformation events such as +earthquakes and aseismic slow slip event. + +The transformation relies on a "master" JSON file, describing general metadata on +the deformation model, its spatial and temporal extent, and listing spatial +components whose values are stored in :ref:`Geodetic TIFF grids (GTG) <geodetictiffgrids>`. +The valuation of each component is modulated by a time function (constant, step, +reverse step, velocity, piecewise, exponential). + +All details on the content of this JSON file are given in the `Proposal for encoding +of a Deformation Model <https://github.com/linz/deformation-model-format/blob/master/doc/JsonGeotiffDeformationModelFormat_20200501.pdf>`__ + +If input coordinates are given in the geographic domain (resp. projected domain), +the output will also be in the geographic domain (resp. projected domain). +The domain should be the corresponding to the source_crs metadata of the model. + +This transformation is a generalization of the :ref:`Kinematic datum shifting utilizing a deformation model <deformation>` transformation. + +Parameters +------------------------------------------------------------------------------- + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +model=<filename> + + Filename to the JSON master file for the deformation model. + + +Example +------------------------------------------------------------------------------- + +Transformating a point with the LINZ NZGD2000 deformation model: + +:: + + echo 166.7133850980 -44.5105886020 293.3700 2007.689 | cct +proj=defmodel +model=nzgd2000-20180701.json diff --git a/docs/source/operations/transformations/deformation.rst b/docs/source/operations/transformations/deformation.rst index 02924a25..c28e9404 100644 --- a/docs/source/operations/transformations/deformation.rst +++ b/docs/source/operations/transformations/deformation.rst @@ -38,6 +38,8 @@ construction of new grids is recommended. Starting with PROJ 7.0, use of a GeoTIFF format is recommended to store both the horizontal and vertical velocities. +More complex deformations can be done with the :ref:`Multi-component time-based deformation model <defmodel>` transformation. + Example ------------------------------------------------------------------------------- diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst index 6bd503d4..9f7900f1 100644 --- a/docs/source/operations/transformations/index.rst +++ b/docs/source/operations/transformations/index.rst @@ -11,6 +11,7 @@ systems are based on different datums. :maxdepth: 1 affine + defmodel deformation geogoffset helmert diff --git a/docs/source/specifications/geodetictiffgrids.rst b/docs/source/specifications/geodetictiffgrids.rst index 861e8b82..7c05867b 100644 --- a/docs/source/specifications/geodetictiffgrids.rst +++ b/docs/source/specifications/geodetictiffgrids.rst @@ -137,9 +137,9 @@ is an easy way to inspect such grid files: * SamplesPerPixel = 1 for vertical shift grids. - In the future, different values of SamplesPerPixel may be used to accommodate - for other needs. For example for deformation models, SamplesPerPixel = 3 to combine - horizontal and vertical adjustments. + * SamplesPerPixel = 3 for deformation models combining + horizontal and vertical adjustments. + And even for the current identified needs of horizontal or vertical shifts, more samples may be present (to indicate for example uncertainties), but will be ignored by PROJ. @@ -215,28 +215,32 @@ is an easy way to inspect such grid files: - ``HORIZONTAL_OFFSET``: implies the presence of at least two samples. The first sample must contain the latitude offset and the second sample must contain the longitude offset. - Corresponds to PROJ ``hgridshift`` method. + Corresponds to PROJ :ref:`hgridshift` method. - ``VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL``: implies the presence of at least one sample. The first sample must contain the vertical adjustment. Must be used when the source/interpolation CRS is a Geographic CRS and the target CRS a Vertical CRS. - Corresponds to PROJ ``vgridshift`` method. + Corresponds to PROJ :ref:`vgridshift` method. - ``VERTICAL_OFFSET_VERTICAL_TO_VERTICAL``: implies the presence of at least one sample. The first sample must contain the vertical adjustment. Must be used when the source and target CRS are Vertical CRS. - Corresponds to PROJ ``vgridshift`` method. + Corresponds to PROJ :ref:`vgridshift` method. - ``GEOCENTRIC_TRANSLATION``: implies the presence of at least 3 samples. The first 3 samples must be respectively the geocentric adjustments along the X, Y and Z axis. Must be used when the source and target CRS are geocentric CRS. The interpolation CRS must be a geographic CRS. - Corresponds to PROJ ``xyzgridshift`` method. + Corresponds to PROJ :ref:`xyzgridshift` method. - ``VELOCITY``: implies the presence of at least 3 samples. The first 3 samples must be respectively the velocities along the E(ast), N(orth), U(p) axis in the local topocentric coordinate system. - Corresponds to PROJ ``deformation`` method. + Corresponds to PROJ :ref:`deformation` method. + + - ``DEFORMATION_MODEL``: implies the presence of the ``DISPLACEMENT_TYPE`` + and ``UNCERTAINTY_TYPE`` metadata items. + Corresponds to PROJ :ref:`defmodel` method. For example: @@ -284,6 +288,11 @@ is an easy way to inspect such grid files: Sample values should be the velocity in a linear/time unit in a ENU local topocentric coordinate system. + + ``east_offset`` / ``north_offset`` / ``vertical_offset``: valid for + TYPE=DEFORMATION_MODEL. + For east_offset and north_offset, the unit might be degree or metre. + For vertical_offset, the unit must be metre. + For example: .. code-block:: xml @@ -337,6 +346,16 @@ is an easy way to inspect such grid files: <Item name="UNITTYPE" sample="0" role="unittype">arc-second</Item> <Item name="UNITTYPE" sample="1" role="unittype">arc-second</Item> + * For TYPE=DEFORMATION_MODEL, the type of the displacement must be specified + with a `Item` whose ``name`` is set to ``DISPLACEMENT_TYPE``. + + The accepted values are: ``HORIZONTAL``, ``VERTICAL``, ``3D`` or ``NONE`` + + * For TYPE=DEFORMATION_MODEL, the type of the uncertainty must be specified + with a `Item` whose ``name`` is set to ``UNCERTAINTY_TYPE``. + + The accepted values are: ``HORIZONTAL``, ``VERTICAL``, ``3D`` or ``NONE`` + * The ``target_crs_epsg_code`` metadata item should be present. For a horizontal shift grid, this is the EPSG code of the target geographic CRS. For a vertical shift grid, this is the diff --git a/scripts/reformat_cpp.sh b/scripts/reformat_cpp.sh index 89b8e4e6..3b626e82 100755 --- a/scripts/reformat_cpp.sh +++ b/scripts/reformat_cpp.sh @@ -21,7 +21,12 @@ for i in "$TOPDIR"/include/proj/*.hpp "$TOPDIR"/include/proj/internal/*.hpp \ "$TOPDIR"/src/tracing.cpp "$TOPDIR"/src/grids.hpp "$TOPDIR"/src/grids.cpp \ "$TOPDIR"/src/filemanager.hpp "$TOPDIR"/src/filemanager.cpp \ "$TOPDIR"/src/networkfilemanager.cpp \ - "$TOPDIR"/src/sqlite3_utils.hpp "$TOPDIR"/src/sqlite3_utils.cpp ; do + "$TOPDIR"/src/sqlite3_utils.hpp "$TOPDIR"/src/sqlite3_utils.cpp \ + "$TOPDIR"/src/transformations/defmodel.hpp \ + "$TOPDIR"/src/transformations/defmodel_exceptions.hpp \ + "$TOPDIR"/src/transformations/defmodel_impl.hpp \ + "$TOPDIR"/src/transformations/defmodel.cpp \ + ; do if ! echo "$i" | grep -q "lru_cache.hpp"; then "$SCRIPT_DIR"/reformat.sh "$i"; fi diff --git a/src/4D_api.cpp b/src/4D_api.cpp index 2b00011d..8f427412 100644 --- a/src/4D_api.cpp +++ b/src/4D_api.cpp @@ -1681,14 +1681,14 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) { grinfo.n_lat = grid.height(); /* cell size */ - grinfo.cs_lon = extent.resLon; - grinfo.cs_lat = extent.resLat; + grinfo.cs_lon = extent.resX; + grinfo.cs_lat = extent.resY; /* bounds of grid */ - grinfo.lowerleft.lam = extent.westLon; - grinfo.lowerleft.phi = extent.southLat; - grinfo.upperright.lam = extent.eastLon; - grinfo.upperright.phi = extent.northLat; + grinfo.lowerleft.lam = extent.west; + grinfo.lowerleft.phi = extent.south; + grinfo.upperright.lam = extent.east; + grinfo.upperright.phi = extent.north; }; { diff --git a/src/Makefile.am b/src/Makefile.am index 3feaa810..9f24a60d 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -196,6 +196,10 @@ libproj_la_SOURCES = \ transformations/molodensky.cpp \ transformations/vgridshift.cpp \ transformations/xyzgridshift.cpp \ + transformations/defmodel.cpp \ + transformations/defmodel.hpp \ + transformations/defmodel_exceptions.hpp \ + transformations/defmodel_impl.hpp \ \ aasincos.cpp adjlon.cpp \ dmstor.cpp auth.cpp \ diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp index 6ce0ab18..8940afde 100644 --- a/src/apps/gie.cpp +++ b/src/apps/gie.cpp @@ -723,18 +723,27 @@ Attempt to interpret args as a PJ_COORD. /* This could be avoided if proj_dmstor used the same proj_strtod() */ /* as gie, but that is not the case (yet). When we remove projects.h */ /* from the public API we can change that. */ + + // Even Rouault: unsure about the above. Coordinates are not necessarily + // geographic coordinates, and the roundtrip through radians for + // big projected coordinates cause inaccuracies, that can cause + // test failures when testing points at edge of grids. + // For example 1501000.0 becomes 1501000.000000000233 double d = proj_strtod(prev, (char **) &endp); - double dms = PJ_TODEG(proj_dmstor (prev, (char **) &dmsendp)); - /* TODO: When projects.h is removed, call proj_dmstor() in all cases */ - if (d != dms && fabs(d) < fabs(dms) && fabs(dms) < fabs(d) + 1) { - d = dms; - endp = dmsendp; + if( *endp != '\0' && !isspace(*endp) ) + { + double dms = PJ_TODEG(proj_dmstor (prev, (char **) &dmsendp)); + /* TODO: When projects.h is removed, call proj_dmstor() in all cases */ + if (d != dms && fabs(d) < fabs(dms) && fabs(dms) < fabs(d) + 1) { + d = dms; + endp = dmsendp; + } + /* A number like -81d00'00.000 will be parsed correctly by both */ + /* proj_strtod and proj_dmstor but only the latter will return */ + /* the correct end-pointer. */ + if (d == dms && endp != dmsendp) + endp = dmsendp; } - /* A number like -81d00'00.000 will be parsed correctly by both */ - /* proj_strtod and proj_dmstor but only the latter will return */ - /* the correct end-pointer. */ - if (d == dms && endp != dmsendp) - endp = dmsendp; /* Break out if there were no more numerals */ if (prev==endp) diff --git a/src/grids.cpp b/src/grids.cpp index 067d3d45..8065813a 100644 --- a/src/grids.cpp +++ b/src/grids.cpp @@ -77,21 +77,21 @@ static void swap_words(void *dataIn, size_t word_size, size_t word_count) // --------------------------------------------------------------------------- bool ExtentAndRes::fullWorldLongitude() const { - return eastLon - westLon + resLon >= 2 * M_PI - 1e-10; + return isGeographic && east - west + resX >= 2 * M_PI - 1e-10; } // --------------------------------------------------------------------------- bool ExtentAndRes::contains(const ExtentAndRes &other) const { - return other.westLon >= westLon && other.eastLon <= eastLon && - other.southLat >= southLat && other.northLat <= northLat; + return other.west >= west && other.east <= east && other.south >= south && + other.north <= north; } // --------------------------------------------------------------------------- bool ExtentAndRes::intersects(const ExtentAndRes &other) const { - return other.westLon < eastLon && westLon <= other.westLon && - other.southLat < northLat && southLat <= other.northLat; + return other.west < east && west <= other.west && other.south < north && + south <= other.north; } // --------------------------------------------------------------------------- @@ -119,12 +119,13 @@ VerticalShiftGrid::~VerticalShiftGrid() = default; static ExtentAndRes globalExtent() { ExtentAndRes extent; - extent.westLon = -M_PI; - extent.southLat = -M_PI / 2; - extent.eastLon = M_PI; - extent.northLat = M_PI / 2; - extent.resLon = M_PI; - extent.resLat = M_PI / 2; + extent.isGeographic = true; + extent.west = -M_PI; + extent.south = -M_PI / 2; + extent.east = M_PI; + extent.north = M_PI / 2; + extent.resX = M_PI; + extent.resY = M_PI / 2; return extent; } @@ -238,12 +239,13 @@ GTXVerticalShiftGrid *GTXVerticalShiftGrid::open(PJ_CONTEXT *ctx, } ExtentAndRes extent; - extent.westLon = xorigin * DEG_TO_RAD; - extent.southLat = yorigin * DEG_TO_RAD; - extent.resLon = xstep * DEG_TO_RAD; - extent.resLat = ystep * DEG_TO_RAD; - extent.eastLon = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD; - extent.northLat = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD; + extent.isGeographic = true; + extent.west = xorigin * DEG_TO_RAD; + extent.south = yorigin * DEG_TO_RAD; + extent.resX = xstep * DEG_TO_RAD; + extent.resY = ystep * DEG_TO_RAD; + extent.east = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD; + extent.north = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD; return new GTXVerticalShiftGrid(ctx, std::move(fp), name, columns, rows, extent); @@ -931,6 +933,10 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { unsigned short count = 0; unsigned short *geokeys = nullptr; bool pixelIsArea = false; + + ExtentAndRes extent; + extent.isGeographic = true; + if (!TIFFGetField(m_hTIFF, TIFFTAG_GEOKEYDIRECTORY, &count, &geokeys)) { pj_log(m_ctx, PJ_LOG_DEBUG_MINOR, "No GeoKeys tag"); } else { @@ -953,6 +959,7 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { for (unsigned int i = 4; i + 3 < count; i += 4) { constexpr unsigned short GTModelTypeGeoKey = 1024; + constexpr unsigned short ModelTypeProjected = 1; constexpr unsigned short ModelTypeGeographic = 2; constexpr unsigned short GTRasterTypeGeoKey = 1025; @@ -960,10 +967,13 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { // constexpr unsigned short RasterPixelIsPoint = 2; if (geokeys[i] == GTModelTypeGeoKey) { - if (geokeys[i + 3] != ModelTypeGeographic) { - pj_log(m_ctx, PJ_LOG_ERROR, "Only GTModelTypeGeoKey = " - "ModelTypeGeographic is " - "supported"); + if (geokeys[i + 3] == ModelTypeProjected) { + extent.isGeographic = false; + } else if (geokeys[i + 3] != ModelTypeGeographic) { + pj_log(m_ctx, PJ_LOG_ERROR, + "Only GTModelTypeGeoKey = " + "ModelTypeGeographic or ModelTypeProjected are " + "supported"); return nullptr; } } else if (geokeys[i] == GTRasterTypeGeoKey) { @@ -976,8 +986,8 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { double hRes = 0; double vRes = 0; - double westLon = 0; - double northLat = 0; + double west = 0; + double north = 0; double *matrix = nullptr; if (TIFFGetField(m_hTIFF, TIFFTAG_GEOTRANSMATRIX, &count, &matrix) && @@ -991,9 +1001,9 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { return nullptr; } - westLon = matrix[3]; + west = matrix[3]; hRes = matrix[0]; - northLat = matrix[7]; + north = matrix[7]; vRes = -matrix[5]; // negation to simulate GeoPixelScale convention } else { double *geopixelscale = nullptr; @@ -1022,34 +1032,33 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() { return nullptr; } - westLon = geotiepoints[3] - geotiepoints[0] * hRes; - northLat = geotiepoints[4] + geotiepoints[1] * vRes; + west = geotiepoints[3] - geotiepoints[0] * hRes; + north = geotiepoints[4] + geotiepoints[1] * vRes; } if (pixelIsArea) { - westLon += 0.5 * hRes; - northLat -= 0.5 * vRes; + west += 0.5 * hRes; + north -= 0.5 * vRes; } - ExtentAndRes extent; - extent.westLon = westLon * DEG_TO_RAD; - extent.northLat = northLat * DEG_TO_RAD; - extent.resLon = hRes * DEG_TO_RAD; - extent.resLat = fabs(vRes) * DEG_TO_RAD; - extent.eastLon = (westLon + hRes * (width - 1)) * DEG_TO_RAD; - extent.southLat = (northLat - vRes * (height - 1)) * DEG_TO_RAD; + const double mulFactor = extent.isGeographic ? DEG_TO_RAD : 1; + extent.west = west * mulFactor; + extent.north = north * mulFactor; + extent.resX = hRes * mulFactor; + extent.resY = fabs(vRes) * mulFactor; + extent.east = (west + hRes * (width - 1)) * mulFactor; + extent.south = (north - vRes * (height - 1)) * mulFactor; if (vRes < 0) { - std::swap(extent.northLat, extent.southLat); + std::swap(extent.north, extent.south); } - if (!(fabs(extent.westLon) <= 4 * M_PI && - fabs(extent.eastLon) <= 4 * M_PI && - fabs(extent.northLat) <= M_PI + 1e-5 && - fabs(extent.southLat) <= M_PI + 1e-5 && - extent.westLon < extent.eastLon && - extent.southLat < extent.northLat && extent.resLon > 1e-10 && - extent.resLat > 1e-10)) { + if (!((!extent.isGeographic || + (fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI && + fabs(extent.north) <= M_PI + 1e-5 && + fabs(extent.south) <= M_PI + 1e-5)) && + extent.west < extent.east && extent.south < extent.north && + extent.resX > 1e-10 && extent.resY > 1e-10)) { pj_log(m_ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", m_filename.c_str()); return nullptr; @@ -1408,17 +1417,19 @@ bool VerticalShiftGridSet::reopen(PJ_CONTEXT *ctx) { // --------------------------------------------------------------------------- -static bool isPointInExtent(double lon, double lat, const ExtentAndRes &extent, +static bool isPointInExtent(double x, double y, const ExtentAndRes &extent, double eps = 0) { - if (!(lat + eps >= extent.southLat && lat - eps <= extent.northLat)) + if (!(y + eps >= extent.south && y - eps <= extent.north)) return false; if (extent.fullWorldLongitude()) return true; - if (lon + eps < extent.westLon) - lon += 2 * M_PI; - else if (lon - eps > extent.eastLon) - lon -= 2 * M_PI; - if (!(lon + eps >= extent.westLon && lon - eps <= extent.eastLon)) + if (extent.isGeographic) { + if (x + eps < extent.west) + x += 2 * M_PI; + else if (x - eps > extent.east) + x -= 2 * M_PI; + } + if (!(x + eps >= extent.west && x - eps <= extent.east)) return false; return true; } @@ -1584,28 +1595,27 @@ NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, } ExtentAndRes extent; - extent.westLon = -to_double(header + 72) * DEG_TO_RAD; - extent.southLat = to_double(header + 24) * DEG_TO_RAD; - extent.eastLon = -to_double(header + 56) * DEG_TO_RAD; - extent.northLat = to_double(header + 40) * DEG_TO_RAD; - extent.resLon = to_double(header + 104) * DEG_TO_RAD; - extent.resLat = to_double(header + 88) * DEG_TO_RAD; - if (!(fabs(extent.westLon) <= 4 * M_PI && - fabs(extent.eastLon) <= 4 * M_PI && - fabs(extent.northLat) <= M_PI + 1e-5 && - fabs(extent.southLat) <= M_PI + 1e-5 && - extent.westLon < extent.eastLon && - extent.southLat < extent.northLat && extent.resLon > 1e-10 && - extent.resLat > 1e-10)) { + extent.isGeographic = true; + extent.west = -to_double(header + 72) * DEG_TO_RAD; + extent.south = to_double(header + 24) * DEG_TO_RAD; + extent.east = -to_double(header + 56) * DEG_TO_RAD; + extent.north = to_double(header + 40) * DEG_TO_RAD; + extent.resX = to_double(header + 104) * DEG_TO_RAD; + extent.resY = to_double(header + 88) * DEG_TO_RAD; + if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI && + fabs(extent.north) <= M_PI + 1e-5 && + fabs(extent.south) <= M_PI + 1e-5 && extent.west < extent.east && + extent.south < extent.north && extent.resX > 1e-10 && + extent.resY > 1e-10)) { pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", filename.c_str()); pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); return nullptr; } const int columns = static_cast<int>( - fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1); + fabs((extent.east - extent.west) / extent.resX + 0.5) + 1); const int rows = static_cast<int>( - fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) + 1); + fabs((extent.north - extent.south) / extent.resY + 0.5) + 1); return new NTv1Grid(ctx, std::move(fp), filename, columns, rows, extent); } @@ -1695,17 +1705,17 @@ CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, } ExtentAndRes extent; - static_assert(sizeof(extent.westLon) == 8, "wrong sizeof"); - static_assert(sizeof(extent.southLat) == 8, "wrong sizeof"); - static_assert(sizeof(extent.resLon) == 8, "wrong sizeof"); - static_assert(sizeof(extent.resLat) == 8, "wrong sizeof"); - memcpy(&extent.westLon, header + 96, 8); - memcpy(&extent.southLat, header + 104, 8); - memcpy(&extent.resLon, header + 112, 8); - memcpy(&extent.resLat, header + 120, 8); - if (!(fabs(extent.westLon) <= 4 * M_PI && - fabs(extent.southLat) <= M_PI + 1e-5 && extent.resLon > 1e-10 && - extent.resLat > 1e-10)) { + extent.isGeographic = true; + static_assert(sizeof(extent.west) == 8, "wrong sizeof"); + static_assert(sizeof(extent.south) == 8, "wrong sizeof"); + static_assert(sizeof(extent.resX) == 8, "wrong sizeof"); + static_assert(sizeof(extent.resY) == 8, "wrong sizeof"); + memcpy(&extent.west, header + 96, 8); + memcpy(&extent.south, header + 104, 8); + memcpy(&extent.resX, header + 112, 8); + memcpy(&extent.resY, header + 120, 8); + if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.south) <= M_PI + 1e-5 && + extent.resX > 1e-10 && extent.resY > 1e-10)) { pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", filename.c_str()); pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); @@ -1719,8 +1729,8 @@ CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp, pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); return nullptr; } - extent.eastLon = extent.westLon + (width - 1) * extent.resLon; - extent.northLat = extent.southLat + (height - 1) * extent.resLon; + extent.east = extent.west + (width - 1) * extent.resX; + extent.north = extent.south + (height - 1) * extent.resX; return new CTable2Grid(ctx, std::move(fp), filename, width, height, extent); } @@ -1902,8 +1912,8 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx, constexpr int OFFSET_GS_COUNT = 8 + 16 * 10; constexpr int OFFSET_SOUTH_LAT = 8 + 16 * 4; if (must_swap) { - // 6 double values: southLat, northLat, eastLon, westLon, resLat, - // resLon + // 6 double values: south, north, east, west, resY, + // resX for (int i = 0; i < 6; i++) { swap_words(header + OFFSET_SOUTH_LAT + 16 * i, sizeof(double), 1); @@ -1915,42 +1925,40 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx, gridName.append(header + 8, 8); ExtentAndRes extent; - extent.southLat = to_double(header + OFFSET_SOUTH_LAT) * DEG_TO_RAD / - 3600.0; /* S_LAT */ - extent.northLat = to_double(header + OFFSET_SOUTH_LAT + 16) * - DEG_TO_RAD / 3600.0; /* N_LAT */ - extent.eastLon = -to_double(header + OFFSET_SOUTH_LAT + 16 * 2) * - DEG_TO_RAD / 3600.0; /* E_LONG */ - extent.westLon = -to_double(header + OFFSET_SOUTH_LAT + 16 * 3) * - DEG_TO_RAD / 3600.0; /* W_LONG */ - extent.resLat = + extent.isGeographic = true; + extent.south = to_double(header + OFFSET_SOUTH_LAT) * DEG_TO_RAD / + 3600.0; /* S_LAT */ + extent.north = to_double(header + OFFSET_SOUTH_LAT + 16) * DEG_TO_RAD / + 3600.0; /* N_LAT */ + extent.east = -to_double(header + OFFSET_SOUTH_LAT + 16 * 2) * + DEG_TO_RAD / 3600.0; /* E_LONG */ + extent.west = -to_double(header + OFFSET_SOUTH_LAT + 16 * 3) * + DEG_TO_RAD / 3600.0; /* W_LONG */ + extent.resY = to_double(header + OFFSET_SOUTH_LAT + 16 * 4) * DEG_TO_RAD / 3600.0; - extent.resLon = + extent.resX = to_double(header + OFFSET_SOUTH_LAT + 16 * 5) * DEG_TO_RAD / 3600.0; - if (!(fabs(extent.westLon) <= 4 * M_PI && - fabs(extent.eastLon) <= 4 * M_PI && - fabs(extent.northLat) <= M_PI + 1e-5 && - fabs(extent.southLat) <= M_PI + 1e-5 && - extent.westLon < extent.eastLon && - extent.southLat < extent.northLat && extent.resLon > 1e-10 && - extent.resLat > 1e-10)) { + if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI && + fabs(extent.north) <= M_PI + 1e-5 && + fabs(extent.south) <= M_PI + 1e-5 && extent.west < extent.east && + extent.south < extent.north && extent.resX > 1e-10 && + extent.resY > 1e-10)) { pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s", filename.c_str()); pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); return nullptr; } const int columns = static_cast<int>( - fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1); + fabs((extent.east - extent.west) / extent.resX + 0.5) + 1); const int rows = static_cast<int>( - fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) + - 1); + fabs((extent.north - extent.south) / extent.resY + 0.5) + 1); pj_log(ctx, PJ_LOG_DEBUG_MINOR, "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", gridName.c_str(), - columns, rows, extent.westLon * RAD_TO_DEG, - extent.southLat * RAD_TO_DEG, extent.eastLon * RAD_TO_DEG, - extent.northLat * RAD_TO_DEG); + columns, rows, extent.west * RAD_TO_DEG, + extent.south * RAD_TO_DEG, extent.east * RAD_TO_DEG, + extent.north * RAD_TO_DEG); unsigned int gs_count; memcpy(&gs_count, header + OFFSET_GS_COUNT, 4); @@ -2393,8 +2401,8 @@ const HorizontalShiftGrid *HorizontalShiftGrid::gridAt(double lon, double lat) const { for (const auto &child : m_children) { const auto &extentChild = child->extentAndRes(); - const double epsilon = (extentChild.resLon + extentChild.resLat) * - REL_TOLERANCE_HGRIDSHIFT; + const double epsilon = + (extentChild.resX + extentChild.resY) * REL_TOLERANCE_HGRIDSHIFT; if (isPointInExtent(lon, lat, extentChild, epsilon)) { return child->gridAt(lon, lat); } @@ -2411,7 +2419,7 @@ const HorizontalShiftGrid *HorizontalShiftGridSet::gridAt(double lon, } const auto &extent = grid->extentAndRes(); const double epsilon = - (extent.resLon + extent.resLat) * REL_TOLERANCE_HGRIDSHIFT; + (extent.resX + extent.resY) * REL_TOLERANCE_HGRIDSHIFT; if (isPointInExtent(lon, lat, extent, epsilon)) { return grid->gridAt(lon, lat); } @@ -2723,11 +2731,11 @@ bool GenericShiftGridSet::reopen(PJ_CONTEXT *ctx) { // --------------------------------------------------------------------------- -const GenericShiftGrid *GenericShiftGrid::gridAt(double lon, double lat) const { +const GenericShiftGrid *GenericShiftGrid::gridAt(double x, double y) const { for (const auto &child : m_children) { const auto &extentChild = child->extentAndRes(); - if (isPointInExtent(lon, lat, extentChild)) { - return child->gridAt(lon, lat); + if (isPointInExtent(x, y, extentChild)) { + return child->gridAt(x, y); } } return this; @@ -2735,15 +2743,14 @@ const GenericShiftGrid *GenericShiftGrid::gridAt(double lon, double lat) const { // --------------------------------------------------------------------------- -const GenericShiftGrid *GenericShiftGridSet::gridAt(double lon, - double lat) const { +const GenericShiftGrid *GenericShiftGridSet::gridAt(double x, double y) const { for (const auto &grid : m_grids) { if (dynamic_cast<NullGenericShiftGrid *>(grid.get())) { return grid.get(); } const auto &extent = grid->extentAndRes(); - if (isPointInExtent(lon, lat, extent)) { - return grid->gridAt(lon, lat); + if (isPointInExtent(x, y, extent)) { + return grid->gridAt(x, y); } } return nullptr; @@ -2872,9 +2879,9 @@ static PJ_LP pj_hgrid_interpolate(PJ_LP t, const HorizontalShiftGrid *grid, int in; const auto &extent = grid->extentAndRes(); - t.lam /= extent.resLon; + t.lam /= extent.resX; indx.lam = std::isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam)); - t.phi /= extent.resLat; + t.phi /= extent.resY; indx.phi = std::isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi)); frct.lam = t.lam - indx.lam; @@ -2959,13 +2966,13 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in, tb = in; const auto *extent = &(grid->extentAndRes()); const double epsilon = - (extent->resLon + extent->resLat) * REL_TOLERANCE_HGRIDSHIFT; - tb.lam -= extent->westLon; + (extent->resX + extent->resY) * REL_TOLERANCE_HGRIDSHIFT; + tb.lam -= extent->west; if (tb.lam + epsilon < 0) tb.lam += 2 * M_PI; - else if (tb.lam - epsilon > extent->eastLon - extent->westLon) + else if (tb.lam - epsilon > extent->east - extent->west) tb.lam -= 2 * M_PI; - tb.phi -= extent->southLat; + tb.phi -= extent->south; t = pj_hgrid_interpolate(tb, grid, true); if (grid->hasChanged()) { @@ -2995,8 +3002,8 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in, /* to fetch a new grid into which iterate... */ if (del.lam == HUGE_VAL) { PJ_LP lp; - lp.lam = t.lam + extent->westLon; - lp.phi = t.phi + extent->southLat; + lp.lam = t.lam + extent->west; + lp.phi = t.phi + extent->south; auto newGrid = findGrid(grids, lp, gridset); if (newGrid == nullptr || newGrid == grid || newGrid->isNullGrid()) break; @@ -3004,15 +3011,15 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in, grid->name().c_str(), newGrid->name().c_str()); grid = newGrid; extent = &(grid->extentAndRes()); - t.lam = lp.lam - extent->westLon; - t.phi = lp.phi - extent->southLat; + t.lam = lp.lam - extent->west; + t.phi = lp.phi - extent->south; tb = in; - tb.lam -= extent->westLon; + tb.lam -= extent->west; if (tb.lam + epsilon < 0) tb.lam += 2 * M_PI; - else if (tb.lam - epsilon > extent->eastLon - extent->westLon) + else if (tb.lam - epsilon > extent->east - extent->west) tb.lam -= 2 * M_PI; - tb.phi -= extent->southLat; + tb.phi -= extent->south; dif.lam = std::numeric_limits<double>::max(); dif.phi = std::numeric_limits<double>::max(); continue; @@ -3041,8 +3048,8 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in, fprintf(stderr, "Inverse grid shift iteration failed, presumably at " "grid edge.\nUsing first approximation.\n"); - in.lam = adjlon(t.lam + extent->westLon); - in.phi = t.phi + extent->southLat; + in.lam = adjlon(t.lam + extent->west); + in.phi = t.phi + extent->south; return in; } @@ -3097,14 +3104,21 @@ PJ_LP pj_hgrid_value(PJ *P, const ListOfHGrids &grids, PJ_LP lp) { /* normalize input to ll origin */ const auto &extent = grid->extentAndRes(); + if (!extent.isGeographic) { + pj_log(P->ctx, PJ_LOG_ERROR, + "Can only handle grids referenced in a geographic CRS"); + pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return out; + } + const double epsilon = - (extent.resLon + extent.resLat) * REL_TOLERANCE_HGRIDSHIFT; - lp.lam -= extent.westLon; + (extent.resX + extent.resY) * REL_TOLERANCE_HGRIDSHIFT; + lp.lam -= extent.west; if (lp.lam + epsilon < 0) lp.lam += 2 * M_PI; - else if (lp.lam - epsilon > extent.eastLon - extent.westLon) + else if (lp.lam - epsilon > extent.east - extent.west) lp.lam -= 2 * M_PI; - lp.phi -= extent.southLat; + lp.phi -= extent.south; out = pj_hgrid_interpolate(lp, grid, false); if (grid->hasChanged()) { @@ -3151,10 +3165,16 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids, } const auto &extent = grid->extentAndRes(); + if (!extent.isGeographic) { + pj_log(ctx, PJ_LOG_ERROR, + "Can only handle grids referenced in a geographic CRS"); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return HUGE_VAL; + } /* Interpolation of a location within the grid */ - double grid_x = (input.lam - extent.westLon) / extent.resLon; - if (input.lam < extent.westLon) { + double grid_x = (input.lam - extent.west) / extent.resX; + if (input.lam < extent.west) { if (extent.fullWorldLongitude()) { // The first fmod goes to ]-lim, lim[ range // So we add lim again to be in ]0, 2*lim[ and fmod again @@ -3162,9 +3182,9 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids, grid->width(), grid->width()); } else { - grid_x = (input.lam + 2 * M_PI - extent.westLon) / extent.resLon; + grid_x = (input.lam + 2 * M_PI - extent.west) / extent.resX; } - } else if (input.lam > extent.eastLon) { + } else if (input.lam > extent.east) { if (extent.fullWorldLongitude()) { // The first fmod goes to ]-lim, lim[ range // So we add lim again to be in ]0, 2*lim[ and fmod again @@ -3172,10 +3192,10 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids, grid->width(), grid->width()); } else { - grid_x = (input.lam - 2 * M_PI - extent.westLon) / extent.resLon; + grid_x = (input.lam - 2 * M_PI - extent.west) / extent.resX; } } - double grid_y = (input.phi - extent.southLat) / extent.resLat; + double grid_y = (input.phi - extent.south) / extent.resY; int grid_ix = static_cast<int>(lround(floor(grid_x))); if (!(grid_ix >= 0 && grid_ix < grid->width())) { // in the unlikely case we end up here... @@ -3342,11 +3362,9 @@ const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids, // Used by +proj=deformation and +proj=xyzgridshift to do bilinear interpolation // on 3 sample values per node. -bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid, - const PJ_LP &lp, int idx1, - int idx2, int idx3, double &v1, - double &v2, double &v3, - bool &must_retry) { +bool pj_bilinear_interpolation_three_samples( + PJ_CONTEXT *ctx, const GenericShiftGrid *grid, const PJ_LP &lp, int idx1, + int idx2, int idx3, double &v1, double &v2, double &v3, bool &must_retry) { must_retry = false; if (grid->isNullGrid()) { v1 = 0.0; @@ -3356,13 +3374,25 @@ bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid, } const auto &extent = grid->extentAndRes(); - double grid_x = (lp.lam - extent.westLon) / extent.resLon; - if (lp.lam < extent.westLon) { - grid_x = (lp.lam + 2 * M_PI - extent.westLon) / extent.resLon; - } else if (lp.lam > extent.eastLon) { - grid_x = (lp.lam - 2 * M_PI - extent.westLon) / extent.resLon; + if (!extent.isGeographic) { + pj_log(ctx, PJ_LOG_ERROR, + "Can only handle grids referenced in a geographic CRS"); + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return false; + } + + // From a input location lp, determine the grid cell into which it falls, + // by identifying the lower-left x,y of it (ix, iy), and the upper-right + // (ix2, iy2) + + double grid_x = (lp.lam - extent.west) / extent.resX; + // Special case for grids with world extent, and dealing with wrap-around + if (lp.lam < extent.west) { + grid_x = (lp.lam + 2 * M_PI - extent.west) / extent.resX; + } else if (lp.lam > extent.east) { + grid_x = (lp.lam - 2 * M_PI - extent.west) / extent.resX; } - double grid_y = (lp.phi - extent.southLat) / extent.resLat; + double grid_y = (lp.phi - extent.south) / extent.resY; int ix = static_cast<int>(grid_x); int iy = static_cast<int>(grid_y); int ix2 = std::min(ix + 1, grid->width() - 1); @@ -3392,6 +3422,7 @@ bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid, return false; } + // Bilinear interpolation double frct_lam = grid_x - ix; double frct_phi = grid_y - iy; double m10 = frct_lam; diff --git a/src/grids.hpp b/src/grids.hpp index 0fd1b7b0..d060fc95 100644 --- a/src/grids.hpp +++ b/src/grids.hpp @@ -37,12 +37,14 @@ NS_PROJ_START struct ExtentAndRes { - double westLon; // in radian - double southLat; // in radian - double eastLon; // in radian - double northLat; // in radian - double resLon; // in radian - double resLat; // in radian + bool isGeographic; // whether extent and resolutions are in a geographic or + // projected CRS + double west; // in radian for geographic, in CRS units otherwise + double south; // in radian for geographic, in CRS units otherwise + double east; // in radian for geographic, in CRS units otherwise + double north; // in radian for geographic, in CRS units otherwise + double resX; // in radian for geographic, in CRS units otherwise + double resY; // in radian for geographic, in CRS units otherwise bool fullWorldLongitude() const; bool contains(const ExtentAndRes &other) const; @@ -188,7 +190,7 @@ class PROJ_GCC_DLL GenericShiftGrid : public Grid { PROJ_FOR_TEST ~GenericShiftGrid() override; - PROJ_FOR_TEST const GenericShiftGrid *gridAt(double lon, double lat) const; + PROJ_FOR_TEST const GenericShiftGrid *gridAt(double x, double y) const; PROJ_FOR_TEST virtual std::string unit(int sample) const = 0; @@ -228,7 +230,7 @@ class PROJ_GCC_DLL GenericShiftGridSet { grids() const { return m_grids; } - PROJ_FOR_TEST const GenericShiftGrid *gridAt(double lon, double lat) const; + PROJ_FOR_TEST const GenericShiftGrid *gridAt(double x, double y) const; PROJ_FOR_TEST virtual void reassign_context(PJ_CONTEXT *ctx); PROJ_FOR_TEST virtual bool reopen(PJ_CONTEXT *ctx); @@ -253,11 +255,9 @@ PJ_LP pj_hgrid_apply(PJ_CONTEXT *ctx, const ListOfHGrids &grids, PJ_LP lp, const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids, const PJ_LP &input, GenericShiftGridSet *&gridSetOut); -bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid, - const PJ_LP &lp, int idx1, - int idx2, int idx3, double &v1, - double &v2, double &v3, - bool &must_retry); +bool pj_bilinear_interpolation_three_samples( + PJ_CONTEXT *ctx, const GenericShiftGrid *grid, const PJ_LP &lp, int idx1, + int idx2, int idx3, double &v1, double &v2, double &v3, bool &must_retry); NS_PROJ_END diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 71561fcd..c318a266 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -184,6 +184,7 @@ set(SRC_LIBPROJ_TRANSFORMATIONS transformations/molodensky.cpp transformations/vgridshift.cpp transformations/xyzgridshift.cpp + transformations/defmodel.cpp ) set(SRC_LIBPROJ_ISO19111 diff --git a/src/pj_list.h b/src/pj_list.h index b8790b45..cf81ae0e 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -32,6 +32,7 @@ PROJ_HEAD(chamb, "Chamberlin Trimetric") PROJ_HEAD(collg, "Collignon") PROJ_HEAD(comill, "Compact Miller") PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)") +PROJ_HEAD(defmodel, "Deformation model") PROJ_HEAD(deformation, "Kinematic grid shift") PROJ_HEAD(denoy, "Denoyer Semi-Elliptical") PROJ_HEAD(eck1, "Eckert I") diff --git a/src/transformations/defmodel.cpp b/src/transformations/defmodel.cpp new file mode 100644 index 00000000..3d0f2a58 --- /dev/null +++ b/src/transformations/defmodel.cpp @@ -0,0 +1,451 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to deformation model + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#define PJ_LIB__ +#define PROJ_COMPILATION + +#include "defmodel.hpp" +#include "filemanager.hpp" +#include "grids.hpp" +#include "proj_internal.h" + +#include <assert.h> + +#include <map> +#include <memory> +#include <utility> + +PROJ_HEAD(defmodel, "Deformation model"); + +using namespace DeformationModel; + +namespace { + +struct Grid : public GridPrototype { + PJ_CONTEXT *ctx; + const NS_PROJ::GenericShiftGrid *realGrid; + mutable bool checkedHorizontal = false; + mutable bool checkedVertical = false; + mutable int sampleX = 0; + mutable int sampleY = 1; + mutable int sampleZ = 2; + + Grid(PJ_CONTEXT *ctxIn, const NS_PROJ::GenericShiftGrid *realGridIn) + : ctx(ctxIn), realGrid(realGridIn) { + minx = realGridIn->extentAndRes().west; + miny = realGridIn->extentAndRes().south; + resx = realGridIn->extentAndRes().resX; + resy = realGridIn->extentAndRes().resY; + width = realGridIn->width(); + height = realGridIn->height(); + } + + bool checkHorizontal(const std::string &expectedUnit) const { + if (!checkedHorizontal) { + const auto samplesPerPixel = realGrid->samplesPerPixel(); + if (samplesPerPixel < 2) { + pj_log(ctx, PJ_LOG_ERROR, + "defmodel: grid %s has not enough samples", + realGrid->name().c_str()); + return false; + } + bool foundDescX = false; + bool foundDescY = false; + bool foundDesc = false; + for (int i = 0; i < samplesPerPixel; i++) { + const auto desc = realGrid->description(i); + if (desc == "east_offset") { + sampleX = i; + foundDescX = true; + } else if (desc == "north_offset") { + sampleY = i; + foundDescY = true; + } + if (!desc.empty()) { + foundDesc = true; + } + } + if (foundDesc && (!foundDescX || !foundDescY)) { + pj_log(ctx, PJ_LOG_ERROR, + "defmodel: grid %s : Found band description, " + "but not the ones expected", + realGrid->name().c_str()); + return false; + } + const auto unit = realGrid->unit(sampleX); + if (!unit.empty() && unit != expectedUnit) { + pj_log(ctx, PJ_LOG_ERROR, "defmodel: grid %s : Only unit=%s " + "currently handled for this mode", + realGrid->name().c_str(), expectedUnit.c_str()); + return false; + } + checkedHorizontal = true; + } + return true; + } + + bool getLonLatOffset(int ix, int iy, double &lonOffsetRadian, + double &latOffsetRadian) const { + if (!checkHorizontal(STR_DEGREE)) { + return false; + } + float lonOffsetDeg; + float latOffsetDeg; + if (!realGrid->valueAt(ix, iy, sampleX, lonOffsetDeg) || + !realGrid->valueAt(ix, iy, sampleY, latOffsetDeg)) { + return false; + } + lonOffsetRadian = lonOffsetDeg * DEG_TO_RAD; + latOffsetRadian = latOffsetDeg * DEG_TO_RAD; + return true; + } + + bool getZOffset(int ix, int iy, double &zOffset) const { + if (!checkedVertical) { + const auto samplesPerPixel = realGrid->samplesPerPixel(); + if (samplesPerPixel == 1) { + sampleZ = 0; + } else if (samplesPerPixel < 3) { + pj_log(ctx, PJ_LOG_ERROR, + "defmodel: grid %s has not enough samples", + realGrid->name().c_str()); + return false; + } + bool foundDesc = false; + bool foundDescZ = false; + for (int i = 0; i < samplesPerPixel; i++) { + const auto desc = realGrid->description(i); + if (desc == "vertical_offset") { + sampleZ = i; + foundDescZ = true; + } + if (!desc.empty()) { + foundDesc = true; + } + } + if (foundDesc && !foundDescZ) { + pj_log(ctx, PJ_LOG_ERROR, + "defmodel: grid %s : Found band description, " + "but not the ones expected", + realGrid->name().c_str()); + return false; + } + const auto unit = realGrid->unit(sampleZ); + if (!unit.empty() && unit != STR_METRE) { + pj_log(ctx, PJ_LOG_ERROR, + "defmodel: grid %s : Only unit=metre currently " + "handled for this mode", + realGrid->name().c_str()); + return false; + } + checkedVertical = true; + } + float zOffsetFloat = 0.0f; + const bool ret = realGrid->valueAt(ix, iy, sampleZ, zOffsetFloat); + zOffset = zOffsetFloat; + return ret; + } + + bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset, + double &northingOffset) const { + if (!checkHorizontal(STR_METRE)) { + return false; + } + float eastingOffsetFloat = 0.0f; + float northingOffsetFloat = 0.0f; + const bool ret = + realGrid->valueAt(ix, iy, sampleX, eastingOffsetFloat) && + realGrid->valueAt(ix, iy, sampleY, northingOffsetFloat); + eastingOffset = eastingOffsetFloat; + northingOffset = northingOffsetFloat; + return ret; + } + + bool getLonLatZOffset(int ix, int iy, double &lonOffsetRadian, + double &latOffsetRadian, double &zOffset) const { + return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) && + getZOffset(ix, iy, zOffset); + } + + bool getEastingNorthingZOffset(int ix, int iy, double &eastingOffset, + double &northingOffset, + double &zOffset) const { + return getEastingNorthingOffset(ix, iy, eastingOffset, + northingOffset) && + getZOffset(ix, iy, zOffset); + } + +#ifdef DEBUG_DEFMODEL + std::string name() const { return realGrid->name(); } +#endif + + private: + Grid(const Grid &) = delete; + Grid &operator=(const Grid &) = delete; +}; + +struct GridSet : public GridSetPrototype<Grid> { + + PJ_CONTEXT *ctx; + std::unique_ptr<NS_PROJ::GenericShiftGridSet> realGridSet; + std::map<const NS_PROJ::GenericShiftGrid *, std::unique_ptr<Grid>> + mapGrids{}; + + GridSet(PJ_CONTEXT *ctxIn, + std::unique_ptr<NS_PROJ::GenericShiftGridSet> &&realGridSetIn) + : ctx(ctxIn), realGridSet(std::move(realGridSetIn)) {} + + const Grid *gridAt(double x, double y) { + const NS_PROJ::GenericShiftGrid *realGrid = realGridSet->gridAt(x, y); + if (!realGrid) { + return nullptr; + } + auto iter = mapGrids.find(realGrid); + if (iter == mapGrids.end()) { + iter = mapGrids + .insert(std::pair<const NS_PROJ::GenericShiftGrid *, + std::unique_ptr<Grid>>( + realGrid, + std::unique_ptr<Grid>(new Grid(ctx, realGrid)))) + .first; + } + return iter->second.get(); + } + + private: + GridSet(const GridSet &) = delete; + GridSet &operator=(const GridSet &) = delete; +}; + +struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> { + + PJ_CONTEXT *ctx; + PJ *cart; + + EvaluatorIface(PJ_CONTEXT *ctxIn, PJ *cartIn) : ctx(ctxIn), cart(cartIn) {} + + ~EvaluatorIface() { + if (cart) + cart->destructor(cart, 0); + } + + std::unique_ptr<GridSet> open(const std::string &filename) { + auto realGridSet = NS_PROJ::GenericShiftGridSet::open(ctx, filename); + if (!realGridSet) { + pj_log(ctx, PJ_LOG_ERROR, "defmodel: cannot open %s", + filename.c_str()); + return nullptr; + } + return std::unique_ptr<GridSet>( + new GridSet(ctx, std::move(realGridSet))); + } + + bool isGeographicCRS(const std::string &crsDef) { + PJ *P = proj_create(ctx, crsDef.c_str()); + if (P == nullptr) { + return true; // reasonable default value + } + const auto type = proj_get_type(P); + bool ret = (type == PJ_TYPE_GEOGRAPHIC_2D_CRS || + type == PJ_TYPE_GEOGRAPHIC_3D_CRS); + proj_destroy(P); + return ret; + } + + void geographicToGeocentric(double lam, double phi, double height, double a, + double b, double /*es*/, double &X, double &Y, + double &Z) { + (void)a; + (void)b; + assert(cart->a == a); + assert(cart->b == b); + PJ_LPZ lpz; + lpz.lam = lam; + lpz.phi = phi; + lpz.z = height; + PJ_XYZ xyz = cart->fwd3d(lpz, cart); + X = xyz.x; + Y = xyz.y; + Z = xyz.z; + } + + void geocentricToGeographic(double X, double Y, double Z, double a, + double b, double /*es*/, double &lam, + double &phi, double &height) { + (void)a; + (void)b; + assert(cart->a == a); + assert(cart->b == b); + PJ_XYZ xyz; + xyz.x = X; + xyz.y = Y; + xyz.z = Z; + PJ_LPZ lpz = cart->inv3d(xyz, cart); + lam = lpz.lam; + phi = lpz.phi; + height = lpz.z; + } + +#ifdef DEBUG_DEFMODEL + void log(const std::string &msg) { + pj_log(ctx, PJ_LOG_TRACE, "%s", msg.c_str()); + } +#endif + + private: + EvaluatorIface(const EvaluatorIface &) = delete; + EvaluatorIface &operator=(const EvaluatorIface &) = delete; +}; + +struct defmodelData { + std::unique_ptr<Evaluator<Grid, GridSet, EvaluatorIface>> evaluator{}; + EvaluatorIface evaluatorIface; + + explicit defmodelData(PJ_CONTEXT *ctx, PJ *cart) + : evaluatorIface(ctx, cart) {} + + defmodelData(const defmodelData &) = delete; + defmodelData &operator=(const defmodelData &) = delete; +}; + +} // namespace + +static PJ *destructor(PJ *P, int errlev) { + if (nullptr == P) + return nullptr; + + auto Q = static_cast<struct defmodelData *>(P->opaque); + delete Q; + P->opaque = nullptr; + + return pj_default_destructor(P, errlev); +} + +static PJ_COORD forward_4d(PJ_COORD in, PJ *P) { + auto *Q = (struct defmodelData *)P->opaque; + + PJ_COORD out; + out.xyzt.t = in.xyzt.t; + + if (!Q->evaluator->forward(Q->evaluatorIface, in.xyzt.x, in.xyzt.y, + in.xyzt.z, in.xyzt.t, out.xyzt.x, out.xyzt.y, + out.xyzt.z)) { + return proj_coord_error(); + } + + return out; +} + +static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) { + auto *Q = (struct defmodelData *)P->opaque; + + PJ_COORD out; + out.xyzt.t = in.xyzt.t; + + if (!Q->evaluator->inverse(Q->evaluatorIface, in.xyzt.x, in.xyzt.y, + in.xyzt.z, in.xyzt.t, out.xyzt.x, out.xyzt.y, + out.xyzt.z)) { + return proj_coord_error(); + } + + return out; +} + +// Function called by proj_assign_context() when a new context is assigned to +// an existing PJ object. Mostly to deal with objects being passed between +// threads. +static void reassign_context(PJ *P, PJ_CONTEXT *ctx) { + auto *Q = (struct defmodelData *)P->opaque; + if (Q->evaluatorIface.ctx != ctx) { + Q->evaluator->clearGridCache(); + Q->evaluatorIface.ctx = ctx; + } +} + +PJ *TRANSFORMATION(defmodel, 1) { + // Pass a dummy ellipsoid definition that will be overridden just afterwards + auto cart = proj_create(P->ctx, "+proj=cart +a=1"); + if (cart == nullptr) + return destructor(P, ENOMEM); + + /* inherit ellipsoid definition from P to Q->cart */ + pj_inherit_ellipsoid_def(P, cart); + + auto Q = new defmodelData(P->ctx, cart); + P->opaque = (void *)Q; + P->destructor = destructor; + P->reassign_context = reassign_context; + + const char *model = pj_param(P->ctx, P->params, "smodel").s; + if (!model) { + proj_log_error(P, "defmodel: +model= should be specified."); + return destructor(P, PJD_ERR_NO_ARGS); + } + + auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, model); + if (nullptr == file) { + proj_log_error(P, "defmodel: Cannot open %s", model); + return destructor(P, PJD_ERR_INVALID_ARG); + } + file->seek(0, SEEK_END); + unsigned long long size = file->tell(); + // Arbitrary threshold to avoid ingesting an arbitrarily large JSON file, + // that could be a denial of service risk. 10 MB should be sufficiently + // large for any valid use ! + if (size > 10 * 1024 * 1024) { + proj_log_error(P, "defmodel: File %s too large", model); + return destructor(P, PJD_ERR_INVALID_ARG); + } + file->seek(0); + std::string jsonStr; + jsonStr.resize(static_cast<size_t>(size)); + if (file->read(&jsonStr[0], jsonStr.size()) != jsonStr.size()) { + proj_log_error(P, "defmodel: Cannot read %s", model); + return destructor(P, PJD_ERR_INVALID_ARG); + } + + try { + Q->evaluator.reset(new Evaluator<Grid, GridSet, EvaluatorIface>( + MasterFile::parse(jsonStr), Q->evaluatorIface, P->a, P->b)); + } catch (const std::exception &e) { + proj_log_error(P, "defmodel: invalid model: %s", e.what()); + return destructor(P, PJD_ERR_INVALID_ARG); + } + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + + if (Q->evaluator->isGeographicCRS()) { + P->left = PJ_IO_UNITS_RADIANS; + P->right = PJ_IO_UNITS_RADIANS; + } else { + P->left = PJ_IO_UNITS_PROJECTED; + P->right = PJ_IO_UNITS_PROJECTED; + } + + return P; +} diff --git a/src/transformations/defmodel.hpp b/src/transformations/defmodel.hpp new file mode 100644 index 00000000..9b28a7d8 --- /dev/null +++ b/src/transformations/defmodel.hpp @@ -0,0 +1,630 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to deformation model + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +/** This file implements the gridded deformation model proposol of + * https://docs.google.com/document/d/1wiyrAmzqh8MZlzHSp3wf594Ob_M1LeFtDA5swuzvLZY + * It is written in a generic way, independent of the rest of PROJ + * infrastructure. + * + * Verbose debugging info can be turned on by setting the DEBUG_DEFMODEL macro + */ + +#ifndef DEFMODEL_HPP +#define DEFMODEL_HPP + +#ifdef PROJ_COMPILATION +#include "proj/internal/include_nlohmann_json.hpp" +#else +#include "nlohmann/json.hpp" +#endif + +#include <algorithm> +#include <cmath> +#include <exception> +#include <limits> +#include <memory> +#include <string> +#include <vector> + +#ifndef DEFORMATON_MODEL_NAMESPACE +#define DEFORMATON_MODEL_NAMESPACE DeformationModel +#endif + +#include "defmodel_exceptions.hpp" + +namespace DEFORMATON_MODEL_NAMESPACE { + +using json = nlohmann::json; + +// --------------------------------------------------------------------------- + +/** Spatial extent as a bounding box. */ +class SpatialExtent { + public: + /** Parse the provided object as an extent. + * + * @throws ParsingException + */ + static SpatialExtent parse(const json &j); + + double minx() const { return mMinx; } + double miny() const { return mMiny; } + double maxx() const { return mMaxx; } + double maxy() const { return mMaxy; } + + double minxNormalized(bool bIsGeographic) const { + return bIsGeographic ? mMinxRad : mMinx; + } + double minyNormalized(bool bIsGeographic) const { + return bIsGeographic ? mMinyRad : mMiny; + } + double maxxNormalized(bool bIsGeographic) const { + return bIsGeographic ? mMaxxRad : mMaxx; + } + double maxyNormalized(bool bIsGeographic) const { + return bIsGeographic ? mMaxyRad : mMaxy; + } + + protected: + friend class MasterFile; + friend class Component; + SpatialExtent() = default; + + private: + double mMinx = std::numeric_limits<double>::quiet_NaN(); + double mMiny = std::numeric_limits<double>::quiet_NaN(); + double mMaxx = std::numeric_limits<double>::quiet_NaN(); + double mMaxy = std::numeric_limits<double>::quiet_NaN(); + double mMinxRad = std::numeric_limits<double>::quiet_NaN(); + double mMinyRad = std::numeric_limits<double>::quiet_NaN(); + double mMaxxRad = std::numeric_limits<double>::quiet_NaN(); + double mMaxyRad = std::numeric_limits<double>::quiet_NaN(); +}; + +// --------------------------------------------------------------------------- + +/** Epoch */ +class Epoch { + public: + /** Constructor from a ISO 8601 date-time. May throw ParsingException */ + explicit Epoch(const std::string &dt = std::string()); + + /** Return ISO 8601 date-time */ + const std::string &toString() const { return mDt; } + + /** Return decimal year */ + double toDecimalYear() const; + + private: + std::string mDt{}; + double mDecimalYear = 0; +}; + +// --------------------------------------------------------------------------- + +/** Component of a deformation model. */ +class Component { + public: + /** Parse the provided object as a component. + * + * @throws ParsingException + */ + static Component parse(const json &j); + + /** Get a text description of the component. */ + const std::string &description() const { return mDescription; } + + /** Get the region within which the component is defined. Outside this + * region the component evaluates to 0. */ + const SpatialExtent &extent() const { return mSpatialExtent; } + + /** Get the displacement parameters defined by the component, one of + * "none", "horizontal", "vertical", and "3d". The "none" option allows + * for a component which defines uncertainty with different grids to those + * defining displacement. */ + const std::string &displacementType() const { return mDisplacementType; } + + /** Get the uncertainty parameters defined by the component, + * one of "none", "horizontal", "vertical", "3d". */ + const std::string &uncertaintyType() const { return mUncertaintyType; } + + /** Get the horizontal uncertainty to use if it is not defined explicitly + * in the spatial model. */ + double horizontalUncertainty() const { return mHorizontalUncertainty; } + + /** Get the vertical uncertainty to use if it is not defined explicitly in + * the spatial model. */ + double verticalUncertainty() const { return mVerticalUncertainty; } + + struct SpatialModel { + /** Specifies the type of the spatial model data file. Initially it + * is proposed that only "GeoTIFF" is supported. */ + std::string type{}; + + /** How values in model should be interpolated. This proposal will + * support "bilinear" and "geocentric_bilinear". */ + std::string interpolationMethod{}; + + /** Specifies location of the spatial model GeoTIFF file relative to + * the master JSON file. */ + std::string filename{}; + + /** A hex encoded MD5 checksum of the grid file that can be used to + * validate that it is the correct version of the file. */ + std::string md5Checksum{}; + }; + + /** Get the spatial model. */ + const SpatialModel &spatialModel() const { return mSpatialModel; } + + /** Generic type for a type functon */ + struct TimeFunction { + std::string type{}; + + virtual ~TimeFunction(); + + virtual double evaluateAt(double dt) const = 0; + + protected: + TimeFunction() = default; + }; + struct ConstantTimeFunction : public TimeFunction { + + virtual double evaluateAt(double dt) const override; + }; + struct VelocityTimeFunction : public TimeFunction { + /** Date/time at which the velocity function is zero. */ + Epoch referenceEpoch{}; + + virtual double evaluateAt(double dt) const override; + }; + + struct StepTimeFunction : public TimeFunction { + /** Epoch at which the step function transitions from 0 to 1. */ + Epoch stepEpoch{}; + + virtual double evaluateAt(double dt) const override; + }; + + struct ReverseStepTimeFunction : public TimeFunction { + /** Epoch at which the reverse step function transitions from 1. to 0 */ + Epoch stepEpoch{}; + + virtual double evaluateAt(double dt) const override; + }; + + struct PiecewiseTimeFunction : public TimeFunction { + /** One of "zero", "constant", and "linear", defines the behaviour of + * the function before the first defined epoch */ + std::string beforeFirst{}; + + /** One of "zero", "constant", and "linear", defines the behaviour of + * the function after the last defined epoch */ + std::string afterLast{}; + + struct EpochScaleFactorTuple { + /** Defines the date/time of the data point. */ + Epoch epoch{}; + + /** Function value at the above epoch */ + double scaleFactor = std::numeric_limits<double>::quiet_NaN(); + }; + + /** A sorted array data points each defined by two elements. + * The array is sorted in order of increasing epoch. + * Note: where the time function includes a step it is represented by + * two consecutive data points with the same epoch. The first defines + * the scale factor that applies before the epoch and the second the + * scale factor that applies after the epoch. */ + std::vector<EpochScaleFactorTuple> model{}; + + virtual double evaluateAt(double dt) const override; + }; + + struct ExponentialTimeFunction : public TimeFunction { + /** The date/time at which the exponential decay starts. */ + Epoch referenceEpoch{}; + + /** The date/time at which the exponential decay ends. */ + Epoch endEpoch{}; + + /** The relaxation constant in years. */ + double relaxationConstant = std::numeric_limits<double>::quiet_NaN(); + + /** The scale factor that applies before the reference epoch. */ + double beforeScaleFactor = std::numeric_limits<double>::quiet_NaN(); + + /** Initial scale factor. */ + double initialScaleFactor = std::numeric_limits<double>::quiet_NaN(); + + /** The scale factor the exponential function approaches. */ + double finalScaleFactor = std::numeric_limits<double>::quiet_NaN(); + + virtual double evaluateAt(double dt) const override; + }; + + /** Get the time function. */ + const TimeFunction *timeFunction() const { return mTimeFunction.get(); } + + private: + Component() = default; + + std::string mDescription{}; + SpatialExtent mSpatialExtent{}; + std::string mDisplacementType{}; + std::string mUncertaintyType{}; + double mHorizontalUncertainty = std::numeric_limits<double>::quiet_NaN(); + double mVerticalUncertainty = std::numeric_limits<double>::quiet_NaN(); + SpatialModel mSpatialModel{}; + std::unique_ptr<TimeFunction> mTimeFunction{}; +}; + +Component::TimeFunction::~TimeFunction() = default; + +// --------------------------------------------------------------------------- + +/** Master file of a deformation model. */ +class MasterFile { + public: + /** Parse the provided serialized JSON content and return an object. + * + * @throws ParsingException + */ + static std::unique_ptr<MasterFile> parse(const std::string &text); + + /** Get file type. Should always be "deformation_model_master_file" */ + const std::string &fileType() const { return mFileType; } + + /** Get the version of the format. At time of writing, only "1.0" is known + */ + const std::string &formatVersion() const { return mFormatVersion; } + + /** Get brief descriptive name of the deformation model. */ + const std::string &name() const { return mName; } + + /** Get a string identifying the version of the deformation model. + * The format for specifying version is defined by the agency + * responsible for the deformation model. */ + const std::string &version() const { return mVersion; } + + /** Get a string identifying the license of the file. + * e.g "Create Commons Attribution 4.0 International" */ + const std::string &license() const { return mLicense; } + + /** Get a text description of the model. Intended to be longer than name() + */ + const std::string &description() const { return mDescription; } + + /** Get a text description of the model. Intended to be longer than name() + */ + const std::string &publicationDate() const { return mPublicationDate; } + + /** Basic information on the agency responsible for the model. */ + struct Authority { + std::string name{}; + std::string url{}; + std::string address{}; + std::string email{}; + }; + + /** Get basic information on the agency responsible for the model. */ + const Authority &authority() const { return mAuthority; } + + /** Hyperlink related to the model. */ + struct Link { + /** URL holding the information */ + std::string href{}; + + /** Relationship to the dataset. e.g. "about", "source", "license", + * "metadata" */ + std::string rel{}; + + /** Mime type */ + std::string type{}; + + /** Description of the link */ + std::string title{}; + }; + + /** Get links to related information. */ + const std::vector<Link> links() const { return mLinks; } + + /** Get a string identifying the source CRS. That is the coordinate + * reference system to which the deformation model applies. Typically + * "EPSG:XXXX" */ + const std::string &sourceCRS() const { return mSourceCRS; } + + /** Get a string identifying the target CRS. That is, for a time + * dependent coordinate transformation, the coordinate reference + * system resulting from applying the deformation. + * Typically "EPSG:XXXX" */ + const std::string &targetCRS() const { return mTargetCRS; } + + /** Get a string identifying the definition CRS. That is, the + * coordinate reference system used to define the component spatial + * models. Typically "EPSG:XXXX" */ + const std::string &definitionCRS() const { return mDefinitionCRS; } + + /** Get the nominal reference epoch of the deformation model. Formated + * as a ISO-8601 date-time. This is not necessarily used to calculate + * the deformation model - each component defines its own time function. */ + const std::string &referenceEpoch() const { return mReferenceEpoch; } + + /** Get the epoch at which the uncertainties of the deformation model + * are calculated. Formated as a ISO-8601 date-time. */ + const std::string &uncertaintyReferenceEpoch() const { + return mUncertaintyReferenceEpoch; + } + + /** Unit of horizontal offsets. Only "metre" and "degree" are supported. */ + const std::string &horizontalOffsetUnit() const { + return mHorizontalOffsetUnit; + } + + /** Unit of vertical offsets. Only "metre" is supported. */ + const std::string &verticalOffsetUnit() const { + return mVerticalOffsetUnit; + } + + /** Type of horizontal uncertainty. e.g "circular 95% confidence limit" */ + const std::string &horizontalUncertaintyType() const { + return mHorizontalUncertaintyType; + } + + /** Unit of horizontal uncertainty. Only "metre" is supported. */ + const std::string &horizontalUncertaintyUnit() const { + return mHorizontalUncertaintyUnit; + } + + /** Type of vertical uncertainty. e.g "circular 95% confidence limit" */ + const std::string &verticalUncertaintyType() const { + return mVerticalUncertaintyType; + } + + /** Unit of vertical uncertainty. Only "metre" is supported. */ + const std::string &verticalUncertaintyUnit() const { + return mVerticalUncertaintyUnit; + } + + /** Defines how the horizontal offsets are applied to geographic + * coordinates. Only "addition" and "geocentric" are supported */ + const std::string &horizontalOffsetMethod() const { + return mHorizontalOffsetMethod; + } + + /** Get the region within which the deformation model is defined. + * It cannot be calculated outside this region */ + const SpatialExtent &extent() const { return mSpatialExtent; } + + /** Defines the range of times for which the model is valid, specified + * by a first and a last value. The deformation model is undefined for + * dates outside this range. */ + struct TimeExtent { + Epoch first{}; + Epoch last{}; + }; + + /** Get the range of times for which the model is valid. */ + const TimeExtent &timeExtent() const { return mTimeExtent; } + + /** Get an array of the components comprising the deformation model. */ + const std::vector<Component> &components() const { return mComponents; } + + private: + MasterFile() = default; + + std::string mFileType{}; + std::string mFormatVersion{}; + std::string mName{}; + std::string mVersion{}; + std::string mLicense{}; + std::string mDescription{}; + std::string mPublicationDate{}; + Authority mAuthority{}; + std::vector<Link> mLinks{}; + std::string mSourceCRS{}; + std::string mTargetCRS{}; + std::string mDefinitionCRS{}; + std::string mReferenceEpoch{}; + std::string mUncertaintyReferenceEpoch{}; + std::string mHorizontalOffsetUnit{}; + std::string mVerticalOffsetUnit{}; + std::string mHorizontalUncertaintyType{}; + std::string mHorizontalUncertaintyUnit{}; + std::string mVerticalUncertaintyType{}; + std::string mVerticalUncertaintyUnit{}; + std::string mHorizontalOffsetMethod{}; + SpatialExtent mSpatialExtent{}; + TimeExtent mTimeExtent{}; + std::vector<Component> mComponents{}; +}; + +// --------------------------------------------------------------------------- + +/** Prototype for a Grid used by GridSet. Intended to be implemented + * by user code */ +struct GridPrototype { + double minx = 0; + double miny = 0; + double resx = 0; + double resy = 0; + int width = 0; + int height = 0; + + bool getLonLatOffset(int /*ix*/, int /*iy*/, double & /*lonOffsetRadian*/, + double & /*latOffsetRadian*/) const { + throw UnimplementedException("getLonLatOffset unimplemented"); + } + + bool getZOffset(int /*ix*/, int /*iy*/, double & /*zOffset*/) const { + throw UnimplementedException("getZOffset unimplemented"); + } + + bool getEastingNorthingOffset(int /*ix*/, int /*iy*/, + double & /*eastingOffset*/, + double & /*northingOffset*/) const { + throw UnimplementedException("getEastingNorthingOffset unimplemented"); + } + + bool getLonLatZOffset(int /*ix*/, int /*iy*/, double & /*lonOffsetRadian*/, + double & /*latOffsetRadian*/, + double & /*zOffset*/) const { + throw UnimplementedException("getLonLatZOffset unimplemented"); +#if 0 + return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) && + getZOffset(ix, iy, zOffset); +#endif + } + + bool getEastingNorthingZOffset(int /*ix*/, int /*iy*/, + double & /*eastingOffset*/, + double & /*northingOffset*/, + double & /*zOffset*/) const { + throw UnimplementedException("getEastingNorthingOffset unimplemented"); +#if 0 + return getEastingNorthingOffset(ix, iy, eastingOffset, + northingOffset) && + getZOffset(ix, iy, zOffset); +#endif + } + +#ifdef DEBUG_DEFMODEL + std::string name() const { + throw UnimplementedException("name() unimplemented"); + } +#endif +}; + +// --------------------------------------------------------------------------- + +/** Prototype for a GridSet used by EvaluatorIface. Intended to be implemented + * by user code */ +template <class Grid = GridPrototype> struct GridSetPrototype { + // The return pointer should remain "stable" over time for a given grid + // of a GridSet. + const Grid *gridAt(double /*x */, double /* y */) { + throw UnimplementedException("gridAt unimplemented"); + } +}; + +// --------------------------------------------------------------------------- + +/** Prototype for a EvaluatorIface used by Evaluator. Intended to be implemented + * by user code */ +template <class Grid = GridPrototype, class GridSet = GridSetPrototype<>> +struct EvaluatorIfacePrototype { + + std::unique_ptr<GridSet> open(const std::string & /* filename*/) { + throw UnimplementedException("open unimplemented"); + } + + void geographicToGeocentric(double /* lam */, double /* phi */, + double /* height*/, double /* a */, + double /* b */, double /*es*/, double & /* X */, + double & /* Y */, double & /* Z */) { + throw UnimplementedException("geographicToGeocentric unimplemented"); + } + + void geocentricToGeographic(double /* X */, double /* Y */, double /* Z */, + double /* a */, double /* b */, double /*es*/, + double & /* lam */, double & /* phi */, + double & /* height*/) { + throw UnimplementedException("geocentricToGeographic unimplemented"); + } + + bool isGeographicCRS(const std::string & /* crsDef */) { + throw UnimplementedException("isGeographicCRS unimplemented"); + } + +#ifdef DEBUG_DEFMODEL + void log(const std::string & /* msg */) { + throw UnimplementedException("log unimplemented"); + } +#endif +}; + +// --------------------------------------------------------------------------- + +/** Internal class to offer caching services over a Component */ +template <class Grid, class GridSet> struct ComponentEx; + +// --------------------------------------------------------------------------- + +/** Class to evaluate the transformation of a coordinate */ +template <class Grid = GridPrototype, class GridSet = GridSetPrototype<>, + class EvaluatorIface = EvaluatorIfacePrototype<>> +class Evaluator { + public: + /** Constructor. May throw EvaluatorException */ + explicit Evaluator(std::unique_ptr<MasterFile> &&model, + EvaluatorIface &iface, double a, double b); + + /** Evaluate displacement of a position given by (x,y,z,t) and + * return it in (x_out,y_out_,z_out). + * For geographic CRS (only supported at that time), x must be a + * longitude, and y a latitude. + */ + bool forward(EvaluatorIface &iface, double x, double y, double z, double t, + double &x_out, double &y_out, double &z_out) { + return forward(iface, x, y, z, t, false, x_out, y_out, z_out); + } + + /** Apply inverse transformation. */ + bool inverse(EvaluatorIface &iface, double x, double y, double z, double t, + double &x_out, double &y_out, double &z_out); + + /** Clear grid cache */ + void clearGridCache(); + + /** Return whether the definition CRS is a geographic CRS */ + bool isGeographicCRS() const { return mIsGeographicCRS; } + + private: + std::unique_ptr<MasterFile> mModel; + const double mA; + const double mB; + const double mEs; + const bool mIsHorizontalUnitDegree; /* degree vs metre */ + const bool mIsAddition; /* addition vs geocentric */ + const bool mIsGeographicCRS; + + bool forward(EvaluatorIface &iface, double x, double y, double z, double t, + bool forInverseComputation, double &x_out, double &y_out, + double &z_out); + + std::vector<std::unique_ptr<ComponentEx<Grid, GridSet>>> mComponents{}; +}; + +// --------------------------------------------------------------------------- + +} // namespace DeformationModel + +// --------------------------------------------------------------------------- + +#include "defmodel_impl.hpp" + +#endif // DEFMODEL_HPP diff --git a/src/transformations/defmodel_exceptions.hpp b/src/transformations/defmodel_exceptions.hpp new file mode 100644 index 00000000..8addf60f --- /dev/null +++ b/src/transformations/defmodel_exceptions.hpp @@ -0,0 +1,81 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to deformation model + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#ifndef DEFORMATON_MODEL_NAMESPACE +#error "Should be included only by defmodel.hpp" +#endif + +#include <exception> + +namespace DEFORMATON_MODEL_NAMESPACE { + +// --------------------------------------------------------------------------- + +/** Parsing exception. */ +class ParsingException : public std::exception { + public: + ParsingException(const std::string &msg) : msg_(msg) {} + const char *what() const noexcept override; + + private: + std::string msg_; +}; + +const char *ParsingException::what() const noexcept { return msg_.c_str(); } + +// --------------------------------------------------------------------------- + +class UnimplementedException : public std::exception { + public: + UnimplementedException(const std::string &msg) : msg_(msg) {} + const char *what() const noexcept override; + + private: + std::string msg_; +}; + +const char *UnimplementedException::what() const noexcept { + return msg_.c_str(); +} + +// --------------------------------------------------------------------------- + +/** Evaluator exception. */ +class EvaluatorException : public std::exception { + public: + EvaluatorException(const std::string &msg) : msg_(msg) {} + const char *what() const noexcept override; + + private: + std::string msg_; +}; + +const char *EvaluatorException::what() const noexcept { return msg_.c_str(); } + +// --------------------------------------------------------------------------- + +} // namespace DEFORMATON_MODEL_NAMESPACE diff --git a/src/transformations/defmodel_impl.hpp b/src/transformations/defmodel_impl.hpp new file mode 100644 index 00000000..a15137d7 --- /dev/null +++ b/src/transformations/defmodel_impl.hpp @@ -0,0 +1,1265 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Functionality related to deformation model + * Author: Even Rouault, <even.rouault at spatialys.com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ + +#ifndef DEFORMATON_MODEL_NAMESPACE +#error "Should be included only by defmodel.hpp" +#endif + +namespace DEFORMATON_MODEL_NAMESPACE { + +// --------------------------------------------------------------------------- + +static const std::string STR_DEGREE("degree"); +static const std::string STR_METRE("metre"); + +static const std::string STR_ADDITION("addition"); +static const std::string STR_GEOCENTRIC("geocentric"); + +static const std::string STR_BILINEAR("bilinear"); +static const std::string STR_GEOCENTRIC_BILINEAR("geocentric_bilinear"); + +static const std::string STR_NONE("none"); +static const std::string STR_HORIZONTAL("horizontal"); +static const std::string STR_VERTICAL("vertical"); +static const std::string STR_3D("3d"); + +constexpr double DEFMODEL_PI = 3.14159265358979323846; +constexpr double DEG_TO_RAD_CONSTANT = 3.14159265358979323846 / 180.; +inline constexpr double DegToRad(double d) { return d * DEG_TO_RAD_CONSTANT; } + +// --------------------------------------------------------------------------- + +enum class DisplacementType { NONE, HORIZONTAL, VERTICAL, THREE_D }; + +// --------------------------------------------------------------------------- + +/** Internal class to offer caching services over a Grid */ +template <class Grid> struct GridEx { + const Grid *grid; + + bool smallResx; // lesser than one degree + + double sinhalfresx; + double coshalfresx; + + double sinresy; + double cosresy; + + int last_ix0 = -1; + int last_iy0 = -1; + double dX00 = 0; + double dY00 = 0; + double dZ00 = 0; + double dX01 = 0; + double dY01 = 0; + double dZ01 = 0; + double dX10 = 0; + double dY10 = 0; + double dZ10 = 0; + double dX11 = 0; + double dY11 = 0; + double dZ11 = 0; + double sinphi0 = 0; + double cosphi0 = 0; + double sinphi1 = 0; + double cosphi1 = 0; + + explicit GridEx(const Grid *gridIn) + : grid(gridIn), smallResx(grid->resx < DegToRad(1)), + sinhalfresx(sin(grid->resx / 2)), coshalfresx(cos(grid->resx / 2)), + sinresy(sin(grid->resy)), cosresy(cos(grid->resy)) {} + + // Return geocentric offset (dX, dY, dZ) relative to a point + // where x0 = -resx / 2 + inline void getBilinearGeocentric(int ix0, int iy0, double de00, + double dn00, double de01, double dn01, + double de10, double dn10, double de11, + double dn11, double m00, double m01, + double m10, double m11, double &dX, + double &dY, double &dZ) { + // If interpolating in the same cell as before, then we can skip + // the recomputation of dXij, dYij and dZij + if (ix0 != last_ix0 || iy0 != last_iy0) { + + last_ix0 = ix0; + if (iy0 != last_iy0) { + const double y0 = grid->miny + iy0 * grid->resy; + sinphi0 = sin(y0); + cosphi0 = cos(y0); + + // Use trigonometric formulas to avoid new calls to sin/cos + sinphi1 = + /* sin(y0+grid->resyRad) = */ sinphi0 * cosresy + + cosphi0 * sinresy; + cosphi1 = + /* cos(y0+grid->resyRad) = */ cosphi0 * cosresy - + sinphi0 * sinresy; + + last_iy0 = iy0; + } + + // Using "traditional" formulas to convert from easting, northing + // offsets to geocentric offsets + const double sinlam00 = -sinhalfresx; + const double coslam00 = coshalfresx; + const double sinphi00 = sinphi0; + const double cosphi00 = cosphi0; + const double dn00sinphi00 = dn00 * sinphi00; + dX00 = -de00 * sinlam00 - dn00sinphi00 * coslam00; + dY00 = de00 * coslam00 - dn00sinphi00 * sinlam00; + dZ00 = dn00 * cosphi00; + + const double sinlam01 = -sinhalfresx; + const double coslam01 = coshalfresx; + const double sinphi01 = sinphi1; + const double cosphi01 = cosphi1; + const double dn01sinphi01 = dn01 * sinphi01; + dX01 = -de01 * sinlam01 - dn01sinphi01 * coslam01; + dY01 = de01 * coslam01 - dn01sinphi01 * sinlam01; + dZ01 = dn01 * cosphi01; + + const double sinlam10 = sinhalfresx; + const double coslam10 = coshalfresx; + const double sinphi10 = sinphi0; + const double cosphi10 = cosphi0; + const double dn10sinphi10 = dn10 * sinphi10; + dX10 = -de10 * sinlam10 - dn10sinphi10 * coslam10; + dY10 = de10 * coslam10 - dn10sinphi10 * sinlam10; + dZ10 = dn10 * cosphi10; + + const double sinlam11 = sinhalfresx; + const double coslam11 = coshalfresx; + const double sinphi11 = sinphi1; + const double cosphi11 = cosphi1; + const double dn11sinphi11 = dn11 * sinphi11; + dX11 = -de11 * sinlam11 - dn11sinphi11 * coslam11; + dY11 = de11 * coslam11 - dn11sinphi11 * sinlam11; + dZ11 = dn11 * cosphi11; + } + + dX = m00 * dX00 + m01 * dX01 + m10 * dX10 + m11 * dX11; + dY = m00 * dY00 + m01 * dY01 + m10 * dY10 + m11 * dY11; + dZ = m00 * dZ00 + m01 * dZ01 + m10 * dZ10 + m11 * dZ11; + } +}; + +// --------------------------------------------------------------------------- + +/** Internal class to offer caching services over a Component */ +template <class Grid, class GridSet> struct ComponentEx { + const Component &component; + + const bool isBilinearInterpolation; /* bilinear vs geocentric_bilinear */ + + const DisplacementType displacementType; + + // Cache + std::unique_ptr<GridSet> gridSet{}; + std::map<const Grid *, GridEx<Grid>> mapGrids{}; + + private: + mutable double mCachedDt = 0; + mutable double mCachedValue = 0; + + static DisplacementType getDisplacementType(const std::string &s) { + if (s == STR_HORIZONTAL) + return DisplacementType::HORIZONTAL; + if (s == STR_VERTICAL) + return DisplacementType::VERTICAL; + if (s == STR_3D) + return DisplacementType::THREE_D; + return DisplacementType::NONE; + } + + public: + explicit ComponentEx(const Component &componentIn) + : component(componentIn), + isBilinearInterpolation( + componentIn.spatialModel().interpolationMethod == STR_BILINEAR), + displacementType(getDisplacementType(component.displacementType())) {} + + double evaluateAt(double dt) const { + if (dt == mCachedDt) + return mCachedValue; + mCachedDt = dt; + mCachedValue = component.timeFunction()->evaluateAt(dt); + return mCachedValue; + } + + void clearGridCache() { + gridSet.reset(); + mapGrids.clear(); + } +}; + +// --------------------------------------------------------------------------- + +/** Converts a ISO8601 date-time string, formatted as "YYYY-MM-DDTHH:MM:SSZ", + * into a decimal year. + * Leap years are taken into account, but not leap seconds. + */ +static double ISO8601ToDecimalYear(const std::string &dt) { + int year, month, day, hour, min, sec; + if (sscanf(dt.c_str(), "%04d-%02d-%02dT%02d:%02d:%02dZ", &year, &month, + &day, &hour, &min, &sec) != 6 || + year < 1582 || // Start of Gregorian calendar + month < 1 || month > 12 || day < 1 || day > 31 || hour < 0 || + hour >= 24 || min < 0 || min >= 60 || sec < 0 || sec >= 61) { + throw ParsingException("Wrong formatting / invalid date-time for " + + dt); + } + const bool isLeapYear = + (((year % 4) == 0 && (year % 100) != 0) || (year % 400) == 0); + // Given the intended use, we omit leap seconds... + const int month_table[2][12] = { + {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}, + {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}}; + int dayInYear = day - 1; + for (int m = 1; m < month; m++) { + dayInYear += month_table[isLeapYear ? 1 : 0][m - 1]; + } + if (day > month_table[isLeapYear ? 1 : 0][month - 1]) { + throw ParsingException("Wrong formatting / invalid date-time for " + + dt); + } + return year + + (dayInYear * 86400 + hour * 3600 + min * 60 + sec) / + (isLeapYear ? 86400. * 366 : 86400. * 365); +} + +// --------------------------------------------------------------------------- + +Epoch::Epoch(const std::string &dt) : mDt(dt) { + if (!dt.empty()) { + mDecimalYear = ISO8601ToDecimalYear(dt); + } +} + +// --------------------------------------------------------------------------- + +double Epoch::toDecimalYear() const { return mDecimalYear; } + +// --------------------------------------------------------------------------- + +static std::string getString(const json &j, const char *key, bool optional) { + if (!j.contains(key)) { + if (optional) { + return std::string(); + } + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json v = j[key]; + if (!v.is_string()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a string"); + } + return v.get<std::string>(); +} + +static std::string getReqString(const json &j, const char *key) { + return getString(j, key, false); +} + +static std::string getOptString(const json &j, const char *key) { + return getString(j, key, true); +} + +// --------------------------------------------------------------------------- + +static double getDouble(const json &j, const char *key, bool optional) { + if (!j.contains(key)) { + if (optional) { + return std::numeric_limits<double>::quiet_NaN(); + } + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json v = j[key]; + if (!v.is_number()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a number"); + } + return v.get<double>(); +} + +static double getReqDouble(const json &j, const char *key) { + return getDouble(j, key, false); +} +static double getOptDouble(const json &j, const char *key) { + return getDouble(j, key, true); +} + +// --------------------------------------------------------------------------- + +static json getObjectMember(const json &j, const char *key) { + if (!j.contains(key)) { + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json obj = j[key]; + if (!obj.is_object()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a object"); + } + return obj; +} + +// --------------------------------------------------------------------------- + +static json getArrayMember(const json &j, const char *key) { + if (!j.contains(key)) { + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json obj = j[key]; + if (!obj.is_array()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a array"); + } + return obj; +} + +// --------------------------------------------------------------------------- + +std::unique_ptr<MasterFile> MasterFile::parse(const std::string &text) { + std::unique_ptr<MasterFile> dmmf(new MasterFile()); + json j; + try { + j = json::parse(text); + } catch (const std::exception &e) { + throw ParsingException(e.what()); + } + if (!j.is_object()) { + throw ParsingException("Not an object"); + } + dmmf->mFileType = getReqString(j, "file_type"); + dmmf->mFormatVersion = getReqString(j, "format_version"); + dmmf->mName = getOptString(j, "name"); + dmmf->mVersion = getOptString(j, "version"); + dmmf->mLicense = getOptString(j, "license"); + dmmf->mDescription = getOptString(j, "description"); + dmmf->mPublicationDate = getOptString(j, "publication_date"); + + if (j.contains("authority")) { + const json jAuthority = j["authority"]; + if (!jAuthority.is_object()) { + throw ParsingException("authority is not a object"); + } + dmmf->mAuthority.name = getOptString(jAuthority, "name"); + dmmf->mAuthority.url = getOptString(jAuthority, "url"); + dmmf->mAuthority.address = getOptString(jAuthority, "address"); + dmmf->mAuthority.email = getOptString(jAuthority, "email"); + } + + if (j.contains("links")) { + const json jLinks = j["links"]; + if (!jLinks.is_array()) { + throw ParsingException("links is not an array"); + } + for (const json &jLink : jLinks) { + if (!jLink.is_object()) { + throw ParsingException("links[] item is not an object"); + } + Link link; + link.href = getOptString(jLink, "href"); + link.rel = getOptString(jLink, "rel"); + link.type = getOptString(jLink, "type"); + link.title = getOptString(jLink, "title"); + dmmf->mLinks.emplace_back(std::move(link)); + } + } + dmmf->mSourceCRS = getReqString(j, "source_crs"); + dmmf->mTargetCRS = getReqString(j, "target_crs"); + dmmf->mDefinitionCRS = getReqString(j, "definition_crs"); + if (dmmf->mSourceCRS != dmmf->mDefinitionCRS) { + throw ParsingException( + "source_crs != definition_crs not currently supported"); + } + dmmf->mReferenceEpoch = getOptString(j, "reference_epoch"); + dmmf->mUncertaintyReferenceEpoch = + getOptString(j, "uncertainty_reference_epoch"); + dmmf->mHorizontalOffsetUnit = getOptString(j, "horizontal_offset_unit"); + if (!dmmf->mHorizontalOffsetUnit.empty() && + dmmf->mHorizontalOffsetUnit != STR_METRE && + dmmf->mHorizontalOffsetUnit != STR_DEGREE) { + throw ParsingException("Unsupported value for horizontal_offset_unit"); + } + dmmf->mVerticalOffsetUnit = getOptString(j, "vertical_offset_unit"); + if (!dmmf->mVerticalOffsetUnit.empty() && + dmmf->mVerticalOffsetUnit != STR_METRE) { + throw ParsingException("Unsupported value for vertical_offset_unit"); + } + dmmf->mHorizontalUncertaintyType = + getOptString(j, "horizontal_uncertainty_type"); + dmmf->mHorizontalUncertaintyUnit = + getOptString(j, "horizontal_uncertainty_unit"); + dmmf->mVerticalUncertaintyType = + getOptString(j, "vertical_uncertainty_type"); + dmmf->mVerticalUncertaintyUnit = + getOptString(j, "vertical_uncertainty_unit"); + dmmf->mHorizontalOffsetMethod = getOptString(j, "horizontal_offset_method"); + if (!dmmf->mHorizontalOffsetMethod.empty() && + dmmf->mHorizontalOffsetMethod != STR_ADDITION && + dmmf->mHorizontalOffsetMethod != STR_GEOCENTRIC) { + throw ParsingException( + "Unsupported value for horizontal_offset_method"); + } + dmmf->mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent")); + + const json jTimeExtent = getObjectMember(j, "time_extent"); + dmmf->mTimeExtent.first = Epoch(getReqString(jTimeExtent, "first")); + dmmf->mTimeExtent.last = Epoch(getReqString(jTimeExtent, "last")); + + const json jComponents = getArrayMember(j, "components"); + for (const json &jComponent : jComponents) { + dmmf->mComponents.emplace_back(Component::parse(jComponent)); + const auto &comp = dmmf->mComponents.back(); + if (comp.displacementType() == STR_HORIZONTAL || + comp.displacementType() == STR_3D) { + if (dmmf->mHorizontalOffsetUnit.empty()) { + throw ParsingException("horizontal_offset_unit should be " + "defined as there is a component with " + "displacement_type = horizontal/3d"); + } + if (dmmf->mHorizontalOffsetMethod.empty()) { + throw ParsingException("horizontal_offset_method should be " + "defined as there is a component with " + "displacement_type = horizontal/3d"); + } + } + if (comp.displacementType() == STR_VERTICAL || + comp.displacementType() == STR_3D) { + if (dmmf->mVerticalOffsetUnit.empty()) { + throw ParsingException("vertical_offset_unit should be defined " + "as there is a component with " + "displacement_type = vertical/3d"); + } + } + if (dmmf->mHorizontalOffsetUnit == STR_DEGREE && + comp.spatialModel().interpolationMethod != STR_BILINEAR) { + throw ParsingException("horizontal_offset_unit = degree can only " + "be used with interpolation_method = " + "bilinear"); + } + } + + if (dmmf->mHorizontalOffsetUnit == STR_DEGREE && + dmmf->mHorizontalOffsetMethod != STR_ADDITION) { + throw ParsingException("horizontal_offset_unit = degree can only be " + "used with horizontal_offset_method = addition"); + } + + return dmmf; +} + +// --------------------------------------------------------------------------- + +SpatialExtent SpatialExtent::parse(const json &j) { + SpatialExtent ex; + + const std::string type = getReqString(j, "type"); + if (type != "bbox") { + throw ParsingException("unsupported type of extent"); + } + + const json jParameter = getObjectMember(j, "parameters"); + const json jBbox = getArrayMember(jParameter, "bbox"); + if (jBbox.size() != 4) { + throw ParsingException("bbox is not an array of 4 numeric elements"); + } + for (int i = 0; i < 4; i++) { + if (!jBbox[i].is_number()) { + throw ParsingException( + "bbox is not an array of 4 numeric elements"); + } + } + ex.mMinx = jBbox[0].get<double>(); + ex.mMiny = jBbox[1].get<double>(); + ex.mMaxx = jBbox[2].get<double>(); + ex.mMaxy = jBbox[3].get<double>(); + + ex.mMinxRad = DegToRad(ex.mMinx); + ex.mMinyRad = DegToRad(ex.mMiny); + ex.mMaxxRad = DegToRad(ex.mMaxx); + ex.mMaxyRad = DegToRad(ex.mMaxy); + + return ex; +} + +// --------------------------------------------------------------------------- + +Component Component::parse(const json &j) { + Component comp; + if (!j.is_object()) { + throw ParsingException("component is not an object"); + } + comp.mDescription = getOptString(j, "description"); + comp.mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent")); + comp.mDisplacementType = getReqString(j, "displacement_type"); + if (comp.mDisplacementType != STR_NONE && + comp.mDisplacementType != STR_HORIZONTAL && + comp.mDisplacementType != STR_VERTICAL && + comp.mDisplacementType != STR_3D) { + throw ParsingException("Unsupported value for displacement_type"); + } + comp.mUncertaintyType = getReqString(j, "uncertainty_type"); + comp.mHorizontalUncertainty = getOptDouble(j, "horizontal_uncertainty"); + comp.mVerticalUncertainty = getOptDouble(j, "vertical_uncertainty"); + + const json jSpatialModel = getObjectMember(j, "spatial_model"); + comp.mSpatialModel.type = getReqString(jSpatialModel, "type"); + comp.mSpatialModel.interpolationMethod = + getReqString(jSpatialModel, "interpolation_method"); + if (comp.mSpatialModel.interpolationMethod != STR_BILINEAR && + comp.mSpatialModel.interpolationMethod != STR_GEOCENTRIC_BILINEAR) { + throw ParsingException("Unsupported value for interpolation_method"); + } + comp.mSpatialModel.filename = getReqString(jSpatialModel, "filename"); + comp.mSpatialModel.md5Checksum = + getOptString(jSpatialModel, "md5_checksum"); + + const json jTimeFunction = getObjectMember(j, "time_function"); + const std::string timeFunctionType = getReqString(jTimeFunction, "type"); + const json jParameters = timeFunctionType == "constant" + ? json() + : getObjectMember(jTimeFunction, "parameters"); + + if (timeFunctionType == "constant") { + std::unique_ptr<ConstantTimeFunction> tf(new ConstantTimeFunction()); + tf->type = timeFunctionType; + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "velocity") { + std::unique_ptr<VelocityTimeFunction> tf(new VelocityTimeFunction()); + tf->type = timeFunctionType; + tf->referenceEpoch = + Epoch(getReqString(jParameters, "reference_epoch")); + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "step") { + std::unique_ptr<StepTimeFunction> tf(new StepTimeFunction()); + tf->type = timeFunctionType; + tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch")); + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "reverse_step") { + std::unique_ptr<ReverseStepTimeFunction> tf( + new ReverseStepTimeFunction()); + tf->type = timeFunctionType; + tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch")); + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "piecewise") { + std::unique_ptr<PiecewiseTimeFunction> tf(new PiecewiseTimeFunction()); + tf->type = timeFunctionType; + tf->beforeFirst = getReqString(jParameters, "before_first"); + if (tf->beforeFirst != "zero" && tf->beforeFirst != "constant" && + tf->beforeFirst != "linear") { + throw ParsingException("Unsupported value for before_first"); + } + tf->afterLast = getReqString(jParameters, "after_last"); + if (tf->afterLast != "zero" && tf->afterLast != "constant" && + tf->afterLast != "linear") { + throw ParsingException("Unsupported value for afterLast"); + } + const json jModel = getArrayMember(jParameters, "model"); + for (const json &jModelElt : jModel) { + if (!jModelElt.is_object()) { + throw ParsingException("model[] element is not an object"); + } + PiecewiseTimeFunction::EpochScaleFactorTuple tuple; + tuple.epoch = Epoch(getReqString(jModelElt, "epoch")); + tuple.scaleFactor = getReqDouble(jModelElt, "scale_factor"); + tf->model.emplace_back(std::move(tuple)); + } + + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "exponential") { + std::unique_ptr<ExponentialTimeFunction> tf( + new ExponentialTimeFunction()); + tf->type = timeFunctionType; + tf->referenceEpoch = + Epoch(getReqString(jParameters, "reference_epoch")); + tf->endEpoch = Epoch(getOptString(jParameters, "end_epoch")); + tf->relaxationConstant = + getReqDouble(jParameters, "relaxation_constant"); + if (tf->relaxationConstant <= 0.0) { + throw ParsingException("Invalid value for relaxation_constant"); + } + tf->beforeScaleFactor = + getReqDouble(jParameters, "before_scale_factor"); + tf->initialScaleFactor = + getReqDouble(jParameters, "initial_scale_factor"); + tf->finalScaleFactor = getReqDouble(jParameters, "final_scale_factor"); + comp.mTimeFunction = std::move(tf); + } else { + throw ParsingException("Unsupported type of time function: " + + timeFunctionType); + } + + return comp; +} + +// --------------------------------------------------------------------------- + +double Component::ConstantTimeFunction::evaluateAt(double) const { return 1.0; } + +// --------------------------------------------------------------------------- + +double Component::VelocityTimeFunction::evaluateAt(double dt) const { + return dt - referenceEpoch.toDecimalYear(); +} + +// --------------------------------------------------------------------------- + +double Component::StepTimeFunction::evaluateAt(double dt) const { + if (dt < stepEpoch.toDecimalYear()) + return 0.0; + return 1.0; +} + +// --------------------------------------------------------------------------- + +double Component::ReverseStepTimeFunction::evaluateAt(double dt) const { + if (dt < stepEpoch.toDecimalYear()) + return -1.0; + return 0.0; +} + +// --------------------------------------------------------------------------- + +double Component::PiecewiseTimeFunction::evaluateAt(double dt) const { + if (model.empty()) { + return 0.0; + } + + const double dt1 = model[0].epoch.toDecimalYear(); + if (dt < dt1) { + if (beforeFirst == "zero") + return 0.0; + if (beforeFirst == "constant" || model.size() == 1) + return model[0].scaleFactor; + + // linear + const double f1 = model[0].scaleFactor; + const double dt2 = model[1].epoch.toDecimalYear(); + const double f2 = model[1].scaleFactor; + if (dt1 == dt2) + return f1; + return (f1 * (dt2 - dt) + f2 * (dt - dt1)) / (dt2 - dt1); + } + for (size_t i = 1; i < model.size(); i++) { + const double dtip1 = model[i].epoch.toDecimalYear(); + if (dt < dtip1) { + const double dti = model[i - 1].epoch.toDecimalYear(); + const double fip1 = model[i].scaleFactor; + const double fi = model[i - 1].scaleFactor; + return (fi * (dtip1 - dt) + fip1 * (dt - dti)) / (dtip1 - dti); + } + } + if (afterLast == "zero") { + return 0.0; + } + if (afterLast == "constant" || model.size() == 1) + return model.back().scaleFactor; + + // linear + const double dtnm1 = model[model.size() - 2].epoch.toDecimalYear(); + const double fnm1 = model[model.size() - 2].scaleFactor; + const double dtn = model.back().epoch.toDecimalYear(); + const double fn = model.back().scaleFactor; + if (dtnm1 == dtn) + return fn; + return (fnm1 * (dtn - dt) + fn * (dt - dtnm1)) / (dtn - dtnm1); +} + +// --------------------------------------------------------------------------- + +double Component::ExponentialTimeFunction::evaluateAt(double dt) const { + const double t0 = referenceEpoch.toDecimalYear(); + if (dt < t0) + return beforeScaleFactor; + if (!endEpoch.toString().empty()) { + dt = std::min(dt, endEpoch.toDecimalYear()); + } + return initialScaleFactor + + (finalScaleFactor - initialScaleFactor) * + (1.0 - std::exp(-(dt - t0) / relaxationConstant)); +} + +// --------------------------------------------------------------------------- + +inline void DeltaEastingNorthingToLongLat(double cosphi, double de, double dn, + double a, double b, double es, + double &dlam, double &dphi) { + const double oneMinuX = es * (1 - cosphi * cosphi); + const double X = 1 - oneMinuX; + const double sqrtX = sqrt(X); +#if 0 + // With es of Earth, absolute/relative error is at most 2e-8 + // const double sqrtX = 1 - 0.5 * oneMinuX * (1 + 0.25 * oneMinuX); +#endif + dlam = de * sqrtX / (a * cosphi); + dphi = dn * a * sqrtX * X / (b * b); +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +Evaluator<Grid, GridSet, EvaluatorIface>::Evaluator( + std::unique_ptr<MasterFile> &&model, EvaluatorIface &iface, double a, + double b) + : mModel(std::move(model)), mA(a), mB(b), mEs(1 - (b * b) / (a * a)), + mIsHorizontalUnitDegree(mModel->horizontalOffsetUnit() == STR_DEGREE), + mIsAddition(mModel->horizontalOffsetMethod() == STR_ADDITION), + mIsGeographicCRS(iface.isGeographicCRS(mModel->definitionCRS())) { + if (!mIsGeographicCRS && mIsHorizontalUnitDegree) { + throw EvaluatorException( + "definition_crs = projected CRS and " + "horizontal_offset_unit = degree are incompatible"); + } + if (!mIsGeographicCRS && !mIsAddition) { + throw EvaluatorException( + "definition_crs = projected CRS and " + "horizontal_offset_method = geocentric are incompatible"); + } + mComponents.reserve(mModel->components().size()); + for (const auto &comp : mModel->components()) { + mComponents.emplace_back(std::unique_ptr<ComponentEx<Grid, GridSet>>( + new ComponentEx<Grid, GridSet>(comp))); + if (!mIsGeographicCRS && !mComponents.back()->isBilinearInterpolation) { + throw EvaluatorException( + "definition_crs = projected CRS and " + "interpolation_method = geocentric_bilinear are incompatible"); + } + } +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +void Evaluator<Grid, GridSet, EvaluatorIface>::clearGridCache() { + for (auto &comp : mComponents) { + comp->clearGridCache(); + } +} + +// --------------------------------------------------------------------------- + +#ifdef DEBUG_DEFMODEL + +static std::string shortName(const Component &comp) { + const auto &desc = comp.description(); + return desc.substr(0, desc.find('\n')) + " (" + + comp.spatialModel().filename + ")"; +} + +static std::string toString(double val) { + char buffer[32]; + snprintf(buffer, sizeof(buffer), "%.9g", val); + return buffer; +} + +#endif + +// --------------------------------------------------------------------------- + +static bool bboxCheck(double &x, double &y, bool forInverseComputation, + const double minx, const double miny, const double maxx, + const double maxy, const double EPS, + const double extraMarginForInverse) { + if (x < minx - EPS || x > maxx + EPS || y < miny - EPS || y > maxy + EPS) { + if (!forInverseComputation) { + return false; + } + // In case of iterative computation for inverse, allow to be a + // slightly bit outside of the grid and clamp to the edges + bool xOk = false; + if (x >= minx - EPS && x <= maxx + EPS) { + xOk = true; + } else if (x > minx - extraMarginForInverse && x < minx) { + x = minx; + xOk = true; + } else if (x < maxx + extraMarginForInverse && x > maxx) { + x = maxx; + xOk = true; + } + + bool yOk = false; + if (y >= miny - EPS && y <= maxy + EPS) { + yOk = true; + } else if (y > miny - extraMarginForInverse && y < miny) { + y = miny; + yOk = true; + } else if (y < maxy + extraMarginForInverse && y > maxy) { + y = maxy; + yOk = true; + } + + return xOk && yOk; + } + return true; +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +bool Evaluator<Grid, GridSet, EvaluatorIface>::forward( + EvaluatorIface &iface, double x, double y, double z, double t, + bool forInverseComputation, double &x_out, double &y_out, double &z_out) + +{ + x_out = x; + y_out = y; + z_out = z; + + const double EPS = mIsGeographicCRS ? 1e-10 : 1e-5; + + // Check against global model spatial extent, potentially wrapping + // longitude to match + { + const auto &extent = mModel->extent(); + const double minx = extent.minxNormalized(mIsGeographicCRS); + const double maxx = extent.maxxNormalized(mIsGeographicCRS); + if (mIsGeographicCRS) { + while (x < minx - EPS) { + x += 2.0 * DEFMODEL_PI; + } + while (x > maxx + EPS) { + x -= 2.0 * DEFMODEL_PI; + } + } + const double miny = extent.minyNormalized(mIsGeographicCRS); + const double maxy = extent.maxyNormalized(mIsGeographicCRS); + const double extraMarginForInverse = + mIsGeographicCRS ? DegToRad(0.1) : 10000; + if (!bboxCheck(x, y, forInverseComputation, minx, miny, maxx, maxy, EPS, + extraMarginForInverse)) { +#ifdef DEBUG_DEFMODEL + iface.log("Calculation point " + toString(x) + "," + toString(y) + + " is outside the extents of the deformation model"); +#endif + return false; + } + } + + // Check against global model temporal extent + { + const auto &timeExtent = mModel->timeExtent(); + if (t < timeExtent.first.toDecimalYear() || + t > timeExtent.last.toDecimalYear()) { +#ifdef DEBUG_DEFMODEL + iface.log("Calculation epoch " + toString(t) + + " is not valid for the deformation model"); +#endif + return false; + } + } + + // For mIsHorizontalUnitDegree + double dlam = 0; + double dphi = 0; + + // For !mIsHorizontalUnitDegree + double de = 0; + double dn = 0; + + double dz = 0; + + bool sincosphiInitialized = false; + double sinphi = 0; + double cosphi = 0; + + for (auto &compEx : mComponents) { + const auto &comp = compEx->component; + if (compEx->displacementType == DisplacementType::NONE) { + continue; + } + const auto &extent = comp.extent(); + double xForGrid = x; + double yForGrid = y; + const double minx = extent.minxNormalized(mIsGeographicCRS); + const double maxx = extent.maxxNormalized(mIsGeographicCRS); + const double miny = extent.minyNormalized(mIsGeographicCRS); + const double maxy = extent.maxyNormalized(mIsGeographicCRS); + const double extraMarginForInverse = 0; + if (!bboxCheck(xForGrid, yForGrid, forInverseComputation, minx, miny, + maxx, maxy, EPS, extraMarginForInverse)) { +#ifdef DEBUG_DEFMODEL + iface.log( + "Skipping component " + shortName(comp) + + " due to point being outside of its declared spatial extent."); +#endif + continue; + } + xForGrid = std::max(xForGrid, minx); + yForGrid = std::max(yForGrid, miny); + xForGrid = std::min(xForGrid, maxx); + yForGrid = std::min(yForGrid, maxy); + const auto tfactor = compEx->evaluateAt(t); + if (tfactor == 0.0) { +#ifdef DEBUG_DEFMODEL + iface.log("Skipping component " + shortName(comp) + + " due to time function evaluating to 0."); +#endif + continue; + } + +#ifdef DEBUG_DEFMODEL + iface.log("Entering component " + shortName(comp) + + " with time function evaluating to " + toString(tfactor) + + "."); +#endif + + if (compEx->gridSet == nullptr) { + compEx->gridSet = iface.open(comp.spatialModel().filename); + if (compEx->gridSet == nullptr) { + return false; + } + } + const Grid *grid = compEx->gridSet->gridAt(xForGrid, yForGrid); + if (grid == nullptr) { +#ifdef DEBUG_DEFMODEL + iface.log("Skipping component " + shortName(comp) + + " due to no grid found for this point in the grid set."); +#endif + continue; + } + if (grid->width < 2 || grid->height < 2) { + return false; + } + const double ix_d = (xForGrid - grid->minx) / grid->resx; + const double iy_d = (yForGrid - grid->miny) / grid->resy; + if (ix_d < -EPS || iy_d < -EPS || ix_d + 1 >= grid->width + EPS || + iy_d + 1 >= grid->height + EPS) { +#ifdef DEBUG_DEFMODEL + iface.log("Skipping component " + shortName(comp) + + " due to point being outside of actual spatial extent of " + "grid " + + grid->name() + "."); +#endif + continue; + } + const int ix0 = std::min(static_cast<int>(ix_d), grid->width - 2); + const int iy0 = std::min(static_cast<int>(iy_d), grid->height - 2); + const int ix1 = ix0 + 1; + const int iy1 = iy0 + 1; + const double frct_x = ix_d - ix0; + const double frct_y = iy_d - iy0; + const double one_minus_frct_x = 1. - frct_x; + const double one_minus_frct_y = 1. - frct_y; + const double m00 = one_minus_frct_x * one_minus_frct_y; + const double m10 = frct_x * one_minus_frct_y; + const double m01 = one_minus_frct_x * frct_y; + const double m11 = frct_x * frct_y; + + if (compEx->displacementType == DisplacementType::VERTICAL) { + double dz00 = 0; + double dz01 = 0; + double dz10 = 0; + double dz11 = 0; + if (!grid->getZOffset(ix0, iy0, dz00) || + !grid->getZOffset(ix1, iy0, dz10) || + !grid->getZOffset(ix0, iy1, dz01) || + !grid->getZOffset(ix1, iy1, dz11)) { + return false; + } + const double dzInterp = + dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + + toString(dzInterp) + "."); +#endif + dz += tfactor * dzInterp; + } else if (mIsHorizontalUnitDegree) { + double dx00 = 0; + double dy00 = 0; + double dx01 = 0; + double dy01 = 0; + double dx10 = 0; + double dy10 = 0; + double dx11 = 0; + double dy11 = 0; + if (compEx->displacementType == DisplacementType::HORIZONTAL) { + if (!grid->getLonLatOffset(ix0, iy0, dx00, dy00) || + !grid->getLonLatOffset(ix1, iy0, dx10, dy10) || + !grid->getLonLatOffset(ix0, iy1, dx01, dy01) || + !grid->getLonLatOffset(ix1, iy1, dx11, dy11)) { + return false; + } + } else /* if (compEx->displacementType == DisplacementType::THREE_D) */ { + double dz00 = 0; + double dz01 = 0; + double dz10 = 0; + double dz11 = 0; + if (!grid->getLonLatZOffset(ix0, iy0, dx00, dy00, dz00) || + !grid->getLonLatZOffset(ix1, iy0, dx10, dy10, dz10) || + !grid->getLonLatZOffset(ix0, iy1, dx01, dy01, dz01) || + !grid->getLonLatZOffset(ix1, iy1, dx11, dy11, dz11)) { + return false; + } + const double dzInterp = + dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + + toString(dzInterp) + "."); +#endif + dz += tfactor * dzInterp; + } + const double dlamInterp = + dx00 * m00 + dx01 * m01 + dx10 * m10 + dx11 * m11; + const double dphiInterp = + dy00 * m00 + dy01 * m01 + dy10 * m10 + dy11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dlamInterp = " + toString(tfactor) + " * " + + toString(dlamInterp) + "."); + iface.log("tfactor * dphiInterp = " + toString(tfactor) + " * " + + toString(dphiInterp) + "."); +#endif + dlam += tfactor * dlamInterp; + dphi += tfactor * dphiInterp; + } else /* horizontal unit is metre */ { + double de00 = 0; + double dn00 = 0; + double de01 = 0; + double dn01 = 0; + double de10 = 0; + double dn10 = 0; + double de11 = 0; + double dn11 = 0; + if (compEx->displacementType == DisplacementType::HORIZONTAL) { + if (!grid->getEastingNorthingOffset(ix0, iy0, de00, dn00) || + !grid->getEastingNorthingOffset(ix1, iy0, de10, dn10) || + !grid->getEastingNorthingOffset(ix0, iy1, de01, dn01) || + !grid->getEastingNorthingOffset(ix1, iy1, de11, dn11)) { + return false; + } + } else /* if (compEx->displacementType == DisplacementType::THREE_D) */ { + double dz00 = 0; + double dz01 = 0; + double dz10 = 0; + double dz11 = 0; + if (!grid->getEastingNorthingZOffset(ix0, iy0, de00, dn00, + dz00) || + !grid->getEastingNorthingZOffset(ix1, iy0, de10, dn10, + dz10) || + !grid->getEastingNorthingZOffset(ix0, iy1, de01, dn01, + dz01) || + !grid->getEastingNorthingZOffset(ix1, iy1, de11, dn11, + dz11)) { + return false; + } + const double dzInterp = + dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + + toString(dzInterp) + "."); +#endif + dz += tfactor * dzInterp; + } + if (compEx->isBilinearInterpolation) { + const double deInterp = + de00 * m00 + de01 * m01 + de10 * m10 + de11 * m11; + const double dnInterp = + dn00 * m00 + dn01 * m01 + dn10 * m10 + dn11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * deInterp = " + toString(tfactor) + " * " + + toString(deInterp) + "."); + iface.log("tfactor * dnInterp = " + toString(tfactor) + " * " + + toString(dnInterp) + "."); +#endif + de += tfactor * deInterp; + dn += tfactor * dnInterp; + } else /* geocentric_bilinear */ { + double dX; + double dY; + double dZ; + + auto iter = compEx->mapGrids.find(grid); + if (iter == compEx->mapGrids.end()) { + GridEx<Grid> gridWithCache(grid); + iter = compEx->mapGrids + .insert(std::pair<const Grid *, GridEx<Grid>>( + grid, std::move(gridWithCache))) + .first; + } + GridEx<Grid> &gridwithCacheRef = iter->second; + + gridwithCacheRef.getBilinearGeocentric( + ix0, iy0, de00, dn00, de01, dn01, de10, dn10, de11, dn11, + m00, m01, m10, m11, dX, dY, dZ); + if (!sincosphiInitialized) { + sincosphiInitialized = true; + sinphi = sin(y); + cosphi = cos(y); + } + const double lam_rel_to_cell_center = + (frct_x - 0.5) * grid->resx; + // Use small-angle approximation of sin/cos when reasonable + // Max abs/rel error on cos is 3.9e-9 and on sin 1.3e-11 + const double sinlam = + gridwithCacheRef.smallResx + ? lam_rel_to_cell_center * + (1 - + (1. / 6) * (lam_rel_to_cell_center * + lam_rel_to_cell_center)) + : sin(lam_rel_to_cell_center); + const double coslam = gridwithCacheRef.smallResx + ? (1 - + 0.5 * (lam_rel_to_cell_center * + lam_rel_to_cell_center)) + : cos(lam_rel_to_cell_center); + + // Convert back from geocentric deltas to easting, northing + // deltas + const double deInterp = -dX * sinlam + dY * coslam; + const double dnInterp = + (-dX * coslam - dY * sinlam) * sinphi + dZ * cosphi; +#ifdef DEBUG_DEFMODEL + iface.log("After geocentric_bilinear interpolation: tfactor * " + "deInterp = " + + toString(tfactor) + " * " + toString(deInterp) + "."); + iface.log("After geocentric_bilinear interpolation: tfactor * " + "dnInterp = " + + toString(tfactor) + " * " + toString(dnInterp) + "."); +#endif + de += tfactor * deInterp; + dn += tfactor * dnInterp; + } + } + } + + // Apply shifts depending on horizontal_offset_unit and + // horizontal_offset_method + if (mIsHorizontalUnitDegree) { + x_out += dlam; + y_out += dphi; + } else { +#ifdef DEBUG_DEFMODEL + iface.log("Total sum of de: " + toString(de)); + iface.log("Total sum of dn: " + toString(dn)); +#endif + if (mIsAddition && !mIsGeographicCRS) { + x_out += de; + y_out += dn; + } else if (mIsAddition) { + // Simple way of adding the offset + if (!sincosphiInitialized) { + cosphi = cos(y); + } + DeltaEastingNorthingToLongLat(cosphi, de, dn, mA, mB, mEs, dlam, + dphi); +#ifdef DEBUG_DEFMODEL + iface.log("Result dlam: " + toString(dlam)); + iface.log("Result dphi: " + toString(dphi)); +#endif + x_out += dlam; + y_out += dphi; + } else { + // Geocentric way of adding the offset + if (!sincosphiInitialized) { + sinphi = sin(y); + cosphi = cos(y); + } + const double sinlam = sin(x); + const double coslam = cos(x); + const double dnsinphi = dn * sinphi; + const double dX = -de * sinlam - dnsinphi * coslam; + const double dY = de * coslam - dnsinphi * sinlam; + const double dZ = dn * cosphi; + + double X; + double Y; + double Z; + iface.geographicToGeocentric(x, y, 0, mA, mB, mEs, X, Y, Z); +#ifdef DEBUG_DEFMODEL + iface.log("Geocentric coordinate before: " + toString(X) + "," + + toString(Y) + "," + toString(Z)); + iface.log("Geocentric shift: " + toString(dX) + "," + toString(dY) + + "," + toString(dZ)); +#endif + X += dX; + Y += dY; + Z += dZ; +#ifdef DEBUG_DEFMODEL + iface.log("Geocentric coordinate after: " + toString(X) + "," + + toString(Y) + "," + toString(Z)); +#endif + + double h_out_ignored; + iface.geocentricToGeographic(X, Y, Z, mA, mB, mEs, x_out, y_out, + h_out_ignored); + } + } +#ifdef DEBUG_DEFMODEL + iface.log("Total sum of dz: " + toString(dz)); +#endif + z_out += dz; + + return true; +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +bool Evaluator<Grid, GridSet, EvaluatorIface>::inverse( + EvaluatorIface &iface, double x, double y, double z, double t, + double &x_out, double &y_out, double &z_out) + +{ + x_out = x; + y_out = y; + z_out = z; + constexpr double EPS_HORIZ = 1e-12; + constexpr double EPS_VERT = 1e-3; + constexpr bool forInverseComputation = true; + for (int i = 0; i < 10; i++) { +#ifdef DEBUG_DEFMODEL + iface.log("Iteration " + std::to_string(i) + ": before foward: x=" + + toString(x_out) + ", y=" + toString(y_out)); +#endif + double x_new; + double y_new; + double z_new; + if (!forward(iface, x_out, y_out, z_out, t, forInverseComputation, + x_new, y_new, z_new)) { + return false; + } +#ifdef DEBUG_DEFMODEL + iface.log("After foward: x=" + toString(x_new) + ", y=" + + toString(y_new)); +#endif + const double dx = x_new - x; + const double dy = y_new - y; + const double dz = z_new - z; + x_out -= dx; + y_out -= dy; + z_out -= dz; + if (std::max(std::fabs(dx), std::fabs(dy)) < EPS_HORIZ && + std::fabs(dz) < EPS_VERT) { + return true; + } + } + return false; +} + +// --------------------------------------------------------------------------- + +} // namespace DEFORMATON_MODEL_NAMESPACE diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp index 8aee50c9..8ce02bee 100644 --- a/src/transformations/deformation.cpp +++ b/src/transformations/deformation.cpp @@ -124,7 +124,7 @@ static bool get_grid_values(PJ* P, } bool must_retry = false; - if( !pj_bilinear_interpolation_three_samples(grid, lp, + if( !pj_bilinear_interpolation_three_samples(P->ctx, grid, lp, sampleE, sampleN, sampleU, vx, vy, vz, must_retry) ) diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp index e1c76518..e37e874d 100644 --- a/src/transformations/xyzgridshift.cpp +++ b/src/transformations/xyzgridshift.cpp @@ -105,7 +105,7 @@ static bool get_grid_values(PJ* P, } bool must_retry = false; - if( !pj_bilinear_interpolation_three_samples(grid, lp, + if( !pj_bilinear_interpolation_three_samples(P->ctx, grid, lp, sampleX, sampleY, sampleZ, dx, dy, dz, must_retry) ) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index dcb4c26d..9fc29103 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -16,6 +16,7 @@ proj_add_gie_test("adams_ws1" "gie/adams_ws1.gie") 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("defmodel" "gie/defmodel.gie") # GIGS tests. Uncommented tests are expected to fail due to issues with # various projections. Should be investigated further and fixed. diff --git a/test/gie/Makefile.am b/test/gie/Makefile.am index 1637a38e..e6bb9582 100644 --- a/test/gie/Makefile.am +++ b/test/gie/Makefile.am @@ -15,8 +15,8 @@ EXTRA_DIST = 4D-API_cs2cs-style.gie \ adams_ws1.gie \ adams_ws2.gie \ guyou.gie \ - peirce_q.gie - + peirce_q.gie \ + defmodel.gie PROJ_LIB ?= ../../data/for_tests @@ -65,4 +65,7 @@ guyou: guyou.gie peirce_q: peirce_q.gie PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $< -check-local: 4D-API-cs2cs-style GDA axisswap builtins deformation ellipsoid more_builtins unitconvert DHDN_ETRS89 geotiff_grids adams_hemi adams_ws1 adams_ws2 guyou peirce_q +defmodel: defmodel.gie + PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $< + +check-local: 4D-API-cs2cs-style GDA axisswap builtins deformation ellipsoid more_builtins unitconvert DHDN_ETRS89 geotiff_grids adams_hemi adams_ws1 adams_ws2 guyou peirce_q defmodel diff --git a/test/gie/defmodel.gie b/test/gie/defmodel.gie new file mode 100644 index 00000000..94dd5784 --- /dev/null +++ b/test/gie/defmodel.gie @@ -0,0 +1,124 @@ + +------------------------------------------------------------------------------- +=============================================================================== +Test +proj=defmodel +=============================================================================== + +<gie-strict> + +# Missing +model +operation +proj=defmodel +expect failure errno no_args + +# +model doesn't point to an existing file +operation +proj=defmodel +model=i_do_not_exist +expect failure errno invalid_arg + +# Not a JSON file +operation +proj=defmodel +model=proj.ini +expect failure errno invalid_arg + +# Horizontal deformation with horizontal unit = degree +operation +proj=defmodel +model=tests/simple_model_degree_horizontal.json +tolerance 0.1 mm +accept 2 49 30 2020 +expect 3 51 30 2020 +roundtrip 1 + +# 3D deformation with horizontal unit = degree +operation +proj=defmodel +model=tests/simple_model_degree_3d.json +tolerance 0.1 mm +accept 2 49 30 2020 +expect 3 51 33 2020 +roundtrip 1 + +# Horizontal deformation with horizontal unit = metre +operation +proj=pipeline +step +inv +proj=merc +step +proj=defmodel +model=tests/simple_model_metre_horizontal.json +step +proj=merc +tolerance 0.1 mm +accept 10 20 30 2020 +expect 11 22 30 2020 +roundtrip 1 + +# 3D deformation with horizontal unit = metre +operation +proj=pipeline +step +inv +proj=merc +step +proj=defmodel +model=tests/simple_model_metre_3d.json +step +proj=merc +tolerance 0.1 mm +accept 10 20 30 2020 +expect 11 22 33 2020 +roundtrip 1 + +# 3D deformation with horizontal unit = metre and a projeced grid +operation +proj=pipeline +step +proj=defmodel +model=tests/simple_model_projected.json +tolerance 0.1 mm + +accept 1500200.0 5400400.0 30 2020 +expect 1500200.588 5400399.722 30.6084 2020 +roundtrip 1 + +# South-west corner +accept 1500000.0 5400000.0 30 2020 +expect 1500000.4 5399999.8 30.84 2020 +roundtrip 1 + +# South-east corner +accept 1501000.0 5400000.0 30 2020 +expect 1501000.5 5399999.75 30.75 2020 +roundtrip 1 + +# North-west corner +accept 1500000.0 5401000.0 30 2020 +expect 1500000.8 5400999.6 30.36 2020 +roundtrip 1 + +# North-east corner +accept 1501000.0 5401000.0 30 2020 +expect 1501001.0 5400999.7 30 2020 +roundtrip 1 + +# Test geocentric addition of components +operation +proj=pipeline +step +inv +proj=merc +step +proj=defmodel +model=tests/simple_model_metre_3d_geocentric.json +step +proj=merc +tolerance 0.1 mm +accept 10 20 30 2020 +expect 11 22 33 2020 +roundtrip 1 + +# Vertical deformation with vertical unit = metre +operation +proj=defmodel +model=tests/simple_model_metre_vertical.json +tolerance 0.1 mm +accept 2 49 30 2020 +expect 2 49 33 2020 +roundtrip 1 + +# Adjust for 360 degree longitude offsets +operation +proj=defmodel +model=tests/simple_model_metre_vertical.json +tolerance 0.1 mm + +accept 362 49 30 2020 +expect 2 49 33 2020 + +operation +proj=defmodel +model=tests/simple_model_wrap_east.json + +accept 165.9 -37.3 10 2020 +expect 165.9 -37.3 10.4525 2020 + +operation +proj=defmodel +model=tests/simple_model_wrap_west.json + +accept 165.9 -37.3 10 2020 +expect 165.9 -37.3 10.4525 2020 + +# Test geocentric bilinear interpolation method +operation +proj=defmodel +model=tests/simple_model_polar.json +tolerance 0.1 mm + +accept 20 -90 15 2020 +expect 27.4743245365 -89.9999747721 18.0000 2020 + +accept 120 -90 15 2020 +expect 27.4737934098 -89.9999747718 18.0000 2020 + +accept 235 -89.5 15 2020 +expect -124.9986638571 -89.5000223708 17.3750 2020 + +accept 45 -89.5 15 2020 +expect 44.9991295392 -89.4999759438 18.5469 2020 + +</gie-strict> diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 62932aae..04c6acb6 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -170,3 +170,14 @@ target_link_libraries(test_network add_test(NAME test_network COMMAND test_network) set_property(TEST test_network PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) + + +add_executable(test_defmodel + main.cpp + test_defmodel.cpp) +target_link_libraries(test_defmodel + GTest::gtest + ${PROJ_LIBRARIES}) +add_test(NAME test_defmodel COMMAND test_defmodel) +set_property(TEST test_defmodel + PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am index 11a6473c..2dd0dbff 100644 --- a/test/unit/Makefile.am +++ b/test/unit/Makefile.am @@ -18,6 +18,7 @@ noinst_PROGRAMS += test_cpp_api noinst_PROGRAMS += gie_self_tests noinst_PROGRAMS += include_proj_h_from_c noinst_PROGRAMS += test_network +noinst_PROGRAMS += test_defmodel pj_transform_test_SOURCES = pj_transform_test.cpp main.cpp pj_transform_test_LDADD = ../../src/libproj.la @GTEST_LIBS@ @@ -70,4 +71,12 @@ test_network_LDADD = ../../src/libproj.la @GTEST_LIBS@ @SQLITE3_LIBS@ @CURL_LIBS test_network-check: test_network PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) PROJ_SOURCE_DATA=$(PROJ_LIB) ./test_network -check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check proj_angular_io_test-check proj_context_test-check test_cpp_api-check gie_self_tests-check test_network-check +test_defmodel_SOURCES = test_defmodel.cpp main.cpp +test_defmodel_LDADD = ../../src/libproj.la @GTEST_LIBS@ + +test_defmodel-check: test_defmodel + PROJ_LIB=$(PROJ_LIB) ./test_defmodel + +check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check \ + proj_angular_io_test-check proj_context_test-check test_cpp_api-check \ + gie_self_tests-check test_network-check test_defmodel-check diff --git a/test/unit/test_defmodel.cpp b/test/unit/test_defmodel.cpp new file mode 100644 index 00000000..b3228fd0 --- /dev/null +++ b/test/unit/test_defmodel.cpp @@ -0,0 +1,1515 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Test deformation model + * Author: Even Rouault <even dot rouault at spatialys dot com> + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault <even dot rouault at spatialys dot com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "gtest_include.h" + +// to be able to use internal::toString +#define FROM_PROJ_CPP +#include "proj/internal/internal.hpp" + +#include "proj.h" +#include "proj_internal.h" + +// Silence C4702 (unreachable code) due to some dummy implementation of the +// interfaces of defmodel.hpp +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable : 4702) +#endif + +#define PROJ_COMPILATION +#define DEFORMATON_MODEL_NAMESPACE TestDeformationModel +#include "transformations/defmodel.hpp" + +using namespace DEFORMATON_MODEL_NAMESPACE; + +namespace { + +constexpr double modelMinX = 158; +constexpr double modelMinY = -58; +constexpr double modelMaxX = 194; +constexpr double modelMaxY = -25; + +static json getMinValidContent() { + json j; + j["file_type"] = "GeoTIFF"; + j["format_version"] = "1.0"; + j["source_crs"] = "EPSG:4959"; + j["target_crs"] = "EPSG:7907"; + j["definition_crs"] = "EPSG:4959"; + j["extent"]["type"] = "bbox"; + j["extent"] + ["parameters"] = {{"bbox", {modelMinX, modelMinY, modelMaxX, modelMaxY}}}; + j["time_extent"]["first"] = "1900-01-01T00:00:00Z"; + j["time_extent"]["last"] = "2050-01-01T00:00:00Z"; + j["components"] = json::array(); + + return j; +} + +// --------------------------------------------------------------------------- + +constexpr int IDX_CONSTANT = 0; +constexpr int IDX_VELOCITY = 1; +constexpr int IDX_STEP = 2; +constexpr int IDX_REVERSE_STEP = 3; +constexpr int IDX_PIECEWISE = 4; +constexpr int IDX_EXPONENTIAL = 5; + +static json getFullValidContent() { + json j(getMinValidContent()); + j["name"] = "name"; + j["version"] = "version"; + j["publication_date"] = "2018-07-01T00:00:00Z"; + j["license"] = "license"; + j["description"] = "description"; + j["authority"]["name"] = "authority_name"; + j["authority"]["url"] = "authority_url"; + j["authority"]["address"] = "authority_address"; + j["authority"]["email"] = "authority_email"; + j["links"] = {{{"href", "href"}, + {"rel", "rel"}, + {"type", "type"}, + {"title", "title"}}}; + j["reference_epoch"] = "2000-01-01T00:00:00Z"; + j["uncertainty_reference_epoch"] = "2018-12-15T00:00:00Z"; + j["horizontal_offset_method"] = "addition"; + j["horizontal_offset_unit"] = "metre"; + j["vertical_offset_unit"] = "metre"; + j["horizontal_uncertainty_type"] = "circular 95% confidence limit"; + j["horizontal_uncertainty_unit"] = "metre"; + j["vertical_uncertainty_type"] = "95% confidence limit"; + j["vertical_uncertainty_unit"] = "metre"; + j["components"] = {{ + {"description", "description"}, + {"displacement_type", "horizontal"}, + {"uncertainty_type", "none"}, + {"horizontal_uncertainty", 0.01}, + {"vertical_uncertainty", 0.02}, + {"extent", + {{"type", "bbox"}, + {"parameters", + {{"bbox", {modelMinX, modelMinY, modelMaxX, modelMaxY}}}}}}, + {"spatial_model", + { + {"type", "GeoTIFF"}, + {"interpolation_method", "bilinear"}, + {"filename", "nzgd2000-ndm-grid02.tif"}, + {"md5_checksum", "49fce8ab267be2c8d00d43683060a032"}, + }}, + {"time_function", + { + {"type", "constant"}, {"parameters", json::object()}, + }}, + }}; + + j["components"].push_back(j["components"][0]); + j["components"][IDX_VELOCITY]["time_function"] = { + {"type", "velocity"}, + {"parameters", {{"reference_epoch", "2000-01-01T00:00:00Z"}}}, + }; + + j["components"].push_back(j["components"][0]); + j["components"][IDX_STEP]["time_function"] = { + {"type", "step"}, + {"parameters", {{"step_epoch", "2000-01-01T00:00:00Z"}}}, + }; + + j["components"].push_back(j["components"][0]); + j["components"][IDX_REVERSE_STEP]["time_function"] = { + {"type", "reverse_step"}, + {"parameters", {{"step_epoch", "2000-01-01T00:00:00Z"}}}, + }; + + j["components"].push_back(j["components"][0]); + j["components"][IDX_PIECEWISE]["time_function"] = { + {"type", "piecewise"}, + {"parameters", + {{"before_first", "zero"}, + {"after_last", "constant"}, + {"model", + {{{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}}, + {{"epoch", "2017-01-01T00:00:00Z"}, {"scale_factor", 1.0}}, + {{"epoch", "2017-01-01T00:00:00Z"}, {"scale_factor", 2.0}}, + {{"epoch", "2018-01-01T00:00:00Z"}, {"scale_factor", 1.0}}}}}}}; + + j["components"].push_back(j["components"][0]); + j["components"][IDX_EXPONENTIAL]["time_function"] = { + {"type", "exponential"}, + {"parameters", + { + {"reference_epoch", "2000-01-01T00:00:00Z"}, + {"end_epoch", "2001-01-01T00:00:00Z"}, + {"relaxation_constant", 2.0}, + {"before_scale_factor", 0.0}, + {"initial_scale_factor", 1.0}, + {"final_scale_factor", 3.0}, + }}, + }; + + return j; +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, basic) { + EXPECT_THROW(MasterFile::parse("foo"), ParsingException); + EXPECT_THROW(MasterFile::parse("null"), ParsingException); + EXPECT_THROW(MasterFile::parse("{}"), ParsingException); + + const auto jMinValid(getMinValidContent()); + { + auto mf = MasterFile::parse(jMinValid.dump()); + EXPECT_EQ(mf->fileType(), "GeoTIFF"); + EXPECT_EQ(mf->formatVersion(), "1.0"); + EXPECT_EQ(mf->sourceCRS(), "EPSG:4959"); + EXPECT_EQ(mf->targetCRS(), "EPSG:7907"); + EXPECT_EQ(mf->definitionCRS(), "EPSG:4959"); + EXPECT_EQ(mf->extent().minx(), modelMinX); + EXPECT_EQ(mf->extent().miny(), modelMinY); + EXPECT_EQ(mf->extent().maxx(), modelMaxX); + EXPECT_EQ(mf->extent().maxy(), modelMaxY); + EXPECT_EQ(mf->timeExtent().first.toString(), "1900-01-01T00:00:00Z"); + EXPECT_EQ(mf->timeExtent().last.toString(), "2050-01-01T00:00:00Z"); + } + + // Check that removing one of each required key causes an exception + for (const auto &kv : jMinValid.items()) { + json jcopy(jMinValid); + jcopy.erase(kv.key()); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["definition_crs"] = "EPSG:4326"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["file_type"] = 1; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"].erase("type"); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"].erase("parameters"); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"]["parameters"].clear(); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"]["parameters"].erase("bbox"); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"]["parameters"]["bbox"] = "foo"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"]["parameters"]["bbox"] = {0, 1, 2}; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["extent"]["parameters"]["bbox"] = {0, 1, 2, "foo"}; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["time_extent"] = "foo"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["time_extent"].erase("first"); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jMinValid); + jcopy["time_extent"].erase("last"); + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, full) { + const auto jFullValid(getFullValidContent()); + auto mf = MasterFile::parse(jFullValid.dump()); + EXPECT_EQ(mf->name(), "name"); + EXPECT_EQ(mf->version(), "version"); + EXPECT_EQ(mf->publicationDate(), "2018-07-01T00:00:00Z"); + EXPECT_EQ(mf->license(), "license"); + EXPECT_EQ(mf->description(), "description"); + EXPECT_EQ(mf->authority().name, "authority_name"); + EXPECT_EQ(mf->authority().url, "authority_url"); + EXPECT_EQ(mf->authority().address, "authority_address"); + EXPECT_EQ(mf->authority().email, "authority_email"); + EXPECT_EQ(mf->links().size(), 1U); + EXPECT_EQ(mf->links()[0].href, "href"); + EXPECT_EQ(mf->links()[0].rel, "rel"); + EXPECT_EQ(mf->links()[0].type, "type"); + EXPECT_EQ(mf->links()[0].title, "title"); + EXPECT_EQ(mf->referenceEpoch(), "2000-01-01T00:00:00Z"); + EXPECT_EQ(mf->uncertaintyReferenceEpoch(), "2018-12-15T00:00:00Z"); + EXPECT_EQ(mf->horizontalOffsetUnit(), "metre"); + EXPECT_EQ(mf->verticalOffsetUnit(), "metre"); + EXPECT_EQ(mf->horizontalUncertaintyType(), "circular 95% confidence limit"); + EXPECT_EQ(mf->horizontalUncertaintyUnit(), "metre"); + EXPECT_EQ(mf->verticalUncertaintyType(), "95% confidence limit"); + EXPECT_EQ(mf->verticalUncertaintyUnit(), "metre"); + EXPECT_EQ(mf->horizontalOffsetMethod(), "addition"); + ASSERT_EQ(mf->components().size(), 6U); + { + const auto &comp = mf->components()[IDX_CONSTANT]; + EXPECT_EQ(comp.description(), "description"); + EXPECT_EQ(comp.displacementType(), "horizontal"); + EXPECT_EQ(comp.uncertaintyType(), "none"); + EXPECT_EQ(comp.horizontalUncertainty(), 0.01); + EXPECT_EQ(comp.verticalUncertainty(), 0.02); + EXPECT_EQ(comp.extent().minx(), modelMinX); + EXPECT_EQ(comp.extent().miny(), modelMinY); + EXPECT_EQ(comp.extent().maxx(), modelMaxX); + EXPECT_EQ(comp.extent().maxy(), modelMaxY); + EXPECT_EQ(comp.spatialModel().type, "GeoTIFF"); + EXPECT_EQ(comp.spatialModel().interpolationMethod, "bilinear"); + EXPECT_EQ(comp.spatialModel().filename, "nzgd2000-ndm-grid02.tif"); + EXPECT_EQ(comp.spatialModel().md5Checksum, + "49fce8ab267be2c8d00d43683060a032"); + ASSERT_NE(comp.timeFunction(), nullptr); + ASSERT_EQ(comp.timeFunction()->type, "constant"); + } + { + const auto &comp = mf->components()[IDX_VELOCITY]; + ASSERT_NE(comp.timeFunction(), nullptr); + ASSERT_EQ(comp.timeFunction()->type, "velocity"); + const auto velocity = + static_cast<const Component::VelocityTimeFunction *>( + comp.timeFunction()); + EXPECT_EQ(velocity->referenceEpoch.toString(), "2000-01-01T00:00:00Z"); + } + { + const auto &comp = mf->components()[IDX_STEP]; + ASSERT_NE(comp.timeFunction(), nullptr); + ASSERT_EQ(comp.timeFunction()->type, "step"); + const auto step = static_cast<const Component::StepTimeFunction *>( + comp.timeFunction()); + EXPECT_EQ(step->stepEpoch.toString(), "2000-01-01T00:00:00Z"); + } + { + const auto &comp = mf->components()[IDX_REVERSE_STEP]; + ASSERT_NE(comp.timeFunction(), nullptr); + ASSERT_EQ(comp.timeFunction()->type, "reverse_step"); + const auto step = + static_cast<const Component::ReverseStepTimeFunction *>( + comp.timeFunction()); + EXPECT_EQ(step->stepEpoch.toString(), "2000-01-01T00:00:00Z"); + } + { + const auto &comp = mf->components()[IDX_PIECEWISE]; + ASSERT_NE(comp.timeFunction(), nullptr); + ASSERT_EQ(comp.timeFunction()->type, "piecewise"); + const auto piecewise = + static_cast<const Component::PiecewiseTimeFunction *>( + comp.timeFunction()); + EXPECT_EQ(piecewise->beforeFirst, "zero"); + EXPECT_EQ(piecewise->afterLast, "constant"); + EXPECT_EQ(piecewise->model.size(), 4U); + EXPECT_EQ(piecewise->model[0].epoch.toString(), "2016-01-01T00:00:00Z"); + EXPECT_EQ(piecewise->model[0].scaleFactor, 0.5); + } + { + const auto &comp = mf->components()[IDX_EXPONENTIAL]; + ASSERT_NE(comp.timeFunction(), nullptr); + ASSERT_EQ(comp.timeFunction()->type, "exponential"); + const auto exponential = + static_cast<const Component::ExponentialTimeFunction *>( + comp.timeFunction()); + EXPECT_EQ(exponential->referenceEpoch.toString(), + "2000-01-01T00:00:00Z"); + EXPECT_EQ(exponential->endEpoch.toString(), "2001-01-01T00:00:00Z"); + EXPECT_EQ(exponential->relaxationConstant, 2.0); + EXPECT_EQ(exponential->beforeScaleFactor, 0.0); + EXPECT_EQ(exponential->initialScaleFactor, 1.0); + EXPECT_EQ(exponential->finalScaleFactor, 3.0); + } +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, error_cases) { + const auto jFullValid(getFullValidContent()); + { + json jcopy(jFullValid); + jcopy["horizontal_offset_method"] = "unsupported"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["horizontal_offset_unit"] = "unsupported"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["vertical_offset_unit"] = "unsupported"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_CONSTANT]["spatial_model"] + ["interpolation_method"] = "unsupported"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_CONSTANT]["displacement_type"] = "unsupported"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["model"] = "foo"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["before_first"] = "illegal"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["after_last"] = "illegal"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + { + json jcopy(jFullValid); + jcopy["components"][0]["time_function"]["type"] = "unknown"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + // Unsupported combination + { + json jcopy(jFullValid); + jcopy["horizontal_offset_method"] = "geocentric"; + EXPECT_NO_THROW(MasterFile::parse(jcopy.dump())); + } + { + json jcopy(jFullValid); + jcopy["horizontal_offset_unit"] = "degree"; + EXPECT_NO_THROW(MasterFile::parse(jcopy.dump())); + } + { + json jcopy(jFullValid); + jcopy["horizontal_offset_method"] = "geocentric"; + jcopy["horizontal_offset_unit"] = "degree"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } + + // Unsupported combination + { + json jcopy(jFullValid); + jcopy["components"][IDX_VELOCITY]["spatial_model"] + ["interpolation_method"] = "geocentric_bilinear"; + EXPECT_NO_THROW(MasterFile::parse(jcopy.dump())); + } + { + json jcopy(jFullValid); + jcopy["horizontal_offset_unit"] = "degree"; + EXPECT_NO_THROW(MasterFile::parse(jcopy.dump())); + } + { + json jcopy(jFullValid); + jcopy["components"][IDX_VELOCITY]["spatial_model"] + ["interpolation_method"] = "geocentric_bilinear"; + jcopy["horizontal_offset_unit"] = "degree"; + EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException); + } +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, ISO8601ToDecimalYear) { + EXPECT_EQ(ISO8601ToDecimalYear("2000-01-01T00:00:00Z"), 2000.0); + EXPECT_EQ(ISO8601ToDecimalYear("2000-02-29T12:00:00Z"), + 2000.0 + ((31 + 28) * 86400. + 12 * 3600) / (366 * 86400)); + EXPECT_EQ(ISO8601ToDecimalYear("2000-12-31T23:59:59Z"), + 2000.0 + (366 * 86400 - 1.) / (366 * 86400)); + EXPECT_EQ(ISO8601ToDecimalYear("2001-01-01T00:00:00Z"), 2001.0); + EXPECT_EQ(ISO8601ToDecimalYear("2001-12-31T23:59:59Z"), + 2001.0 + (365 * 86400 - 1.) / (365 * 86400)); + EXPECT_THROW(ISO8601ToDecimalYear(""), ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("0000-01-01T00:00:00Z"), + ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("2001-02-29T00:00:00Z"), + ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("2000-13-01T00:00:00Z"), + ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("2000-01-32T00:00:00Z"), + ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T24:00:00Z"), + ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T00:60:00Z"), + ParsingException); + EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T00:00:61Z"), + ParsingException); +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluate_constant) { + const auto jFullValid(getFullValidContent()); + const auto mf = MasterFile::parse(jFullValid.dump()); + const auto &comp = mf->components()[IDX_CONSTANT]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.0), 1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.0), 1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.0), 1.0); +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluate_velocity) { + const auto jFullValid(getFullValidContent()); + const auto mf = MasterFile::parse(jFullValid.dump()); + const auto &comp = mf->components()[IDX_VELOCITY]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.0), -1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.0), 0.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.0), 1.0); +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluate_step) { + const auto jFullValid(getFullValidContent()); + const auto mf = MasterFile::parse(jFullValid.dump()); + const auto &comp = mf->components()[IDX_STEP]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), 0.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.01), 1.0); +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluate_reverse_step) { + const auto jFullValid(getFullValidContent()); + const auto mf = MasterFile::parse(jFullValid.dump()); + const auto &comp = mf->components()[IDX_REVERSE_STEP]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), -1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 0.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.01), 0.0); +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluate_piecewise) { + const auto jFullValid(getFullValidContent()); + { + const auto mf = MasterFile::parse(jFullValid.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.99), 0.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.00), 0.5); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 0.75); + EXPECT_NEAR(comp.timeFunction()->evaluateAt(2017 - 1e-9), 1.0, 1e-9); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2017.0), 2.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2017.5), 1.5); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.0), 1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2019.0), 1.0); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["before_first"] = "zero"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.0); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["before_first"] = "constant"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["before_first"] = "linear"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.25); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["after_last"] = "zero"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 0.0); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["after_last"] = "constant"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 1.0); + } + + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["after_last"] = "linear"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 0.5); + } + + // No epoch + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["model"] + .clear(); + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.0); + } + + // Just one epoch + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["model"] = { + {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}}}; + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["before_first"] = "linear"; + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["after_last"] = "linear"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 0.5); + } + + // Two identical epochs + { + json jcopy(jFullValid); + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["model"] = { + {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}}, + {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 1.0}}}; + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["before_first"] = "linear"; + jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"] + ["after_last"] = "linear"; + const auto mf = MasterFile::parse(jcopy.dump()); + const auto &comp = mf->components()[IDX_PIECEWISE]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 1.0); + } +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluate_exponential) { + const auto jFullValid(getFullValidContent()); + const auto mf = MasterFile::parse(jFullValid.dump()); + const auto &comp = mf->components()[IDX_EXPONENTIAL]; + EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), 0.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 1.0); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.50), + 1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2000.50 - 2000.00) / 2.0))); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.00), + 1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2001.00 - 2000.00) / 2.0))); + EXPECT_EQ(comp.timeFunction()->evaluateAt(2002.00), + 1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2001.00 - 2000.00) / 2.0))); +} + +// --------------------------------------------------------------------------- + +inline double RadToDeg(double d) { return d / DEG_TO_RAD_CONSTANT; } + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluator_horizontal_unit_degree) { + json j(getMinValidContent()); + j["horizontal_offset_method"] = "addition"; + j["horizontal_offset_unit"] = "degree"; + constexpr double tFactor = 0.5; + constexpr double gridMinX = 160; + constexpr double gridMinY = -50; + constexpr double gridMaxX = 190; + constexpr double gridMaxY = -30; + j["components"] = { + {{"displacement_type", "horizontal"}, + {"uncertainty_type", "none"}, + {"extent", + {{"type", "bbox"}, + {"parameters", + {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}}}}, + {"spatial_model", + { + {"type", "GeoTIFF"}, + {"interpolation_method", "bilinear"}, + {"filename", "bla.tif"}, + }}, + {"time_function", + {{"type", "piecewise"}, + {"parameters", + {{"before_first", "zero"}, + {"after_last", "zero"}, + {"model", + {{{"epoch", "2010-01-01T00:00:00Z"}, {"scale_factor", tFactor}}, + {{"epoch", "2020-01-01T00:00:00Z"}, + {"scale_factor", tFactor}}}}}}}}}}; + + constexpr int iQueriedX = 1; + constexpr int iQueriedY = 3; + constexpr double lonOffsetQueriedX = 0.01; + constexpr double lonOffsetQueriedXp1 = 0.02; + constexpr double latOffsetQueriedY = 0.03; + constexpr double latOffsetQueriedYp1 = 0.04; + constexpr double zOffsetQueriedXY = 10.; + constexpr double zOffsetQueriedXp1Y = 11.; + constexpr double zOffsetQueriedXYp1 = 11.; + constexpr double zOffsetQueriedXp1Yp1 = 12.; + constexpr double gridResX = 2; + constexpr double gridResY = 0.5; + + struct Grid : public GridPrototype { + bool getLonLatOffset(int ix, int iy, double &lonOffsetRadian, + double &latOffsetRadian) const { + if (ix == iQueriedX) { + lonOffsetRadian = DegToRad(lonOffsetQueriedX); + } else if (ix == iQueriedX + 1) { + lonOffsetRadian = DegToRad(lonOffsetQueriedXp1); + } else { + return false; + } + if (iy == iQueriedY) { + latOffsetRadian = DegToRad(latOffsetQueriedY); + } else if (iy == iQueriedY + 1) { + latOffsetRadian = DegToRad(latOffsetQueriedYp1); + } else { + return false; + } + return true; + } + + bool getZOffset(int ix, int iy, double &zOffset) const { + if (ix == iQueriedX && iy == iQueriedY) { + zOffset = zOffsetQueriedXY; + } else if (ix == iQueriedX + 1 && iy == iQueriedY) { + zOffset = zOffsetQueriedXp1Y; + } else if (ix == iQueriedX && iy == iQueriedY + 1) { + zOffset = zOffsetQueriedXYp1; + } else if (ix == iQueriedX + 1 && iy == iQueriedY + 1) { + zOffset = zOffsetQueriedXp1Yp1; + } else { + return false; + } + return true; + } + + bool getLonLatZOffset(int ix, int iy, double &lonOffsetRadian, + double &latOffsetRadian, double &zOffset) const { + return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) && + getZOffset(ix, iy, zOffset); + } + +#ifdef DEBUG_DEFMODEL + std::string name() const { return std::string(); } +#endif + }; + + struct GridSet : public GridSetPrototype<Grid> { + + Grid grid{}; + + GridSet() { + grid.minx = DegToRad(gridMinX); + grid.miny = DegToRad(gridMinY); + grid.resx = DegToRad(gridResX); + grid.resy = DegToRad(gridResY); + grid.width = + 1 + static_cast<int>(0.5 + (gridMaxX - gridMinX) / gridResX); + grid.height = + 1 + static_cast<int>(0.5 + (gridMaxY - gridMinY) / gridResY); + } + + const Grid *gridAt(double /*x */, double /* y */) { return &grid; } + }; + + struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> { + std::unique_ptr<GridSet> open(const std::string &filename) { + if (filename != "bla.tif") + return nullptr; + return std::unique_ptr<GridSet>(new GridSet()); + } + + bool isGeographicCRS(const std::string & /* crsDef */) { return true; } + +#ifdef DEBUG_DEFMODEL + void log(const std::string & /* msg */) {} +#endif + }; + + EvaluatorIface iface; + + Evaluator<Grid, GridSet, EvaluatorIface> eval(MasterFile::parse(j.dump()), + iface, 1, 1); + double newLon; + double newLat; + double newZ; + constexpr double tValid = 2018; + constexpr double EPS = 1e-9; + constexpr double zVal = 100; + + // Query on exact grid intersection + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon + tFactor * lonOffsetQueriedX, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat + tFactor * latOffsetQueriedY, EPS); + EXPECT_EQ(newZ, zVal); + } + + // Query between grid points + { + constexpr double alphaX = 0.25; + constexpr double alphaY = 0.125; + const double lon = gridMinX + iQueriedX * gridResX + alphaX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR( + RadToDeg(newLon), + lon + + tFactor * (lonOffsetQueriedX + + alphaX * (lonOffsetQueriedXp1 - lonOffsetQueriedX)), + EPS); + EXPECT_NEAR( + RadToDeg(newLat), + lat + + tFactor * (latOffsetQueriedY + + alphaY * (latOffsetQueriedYp1 - latOffsetQueriedY)), + EPS); + EXPECT_EQ(newZ, zVal); + } + + // Longitude < model min + { + const double lon = modelMinX - 1e-1; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + } + + // Longitude > model max + { + const double lon = modelMaxX + 1e-1; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + } + + // Latitude < model min + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = modelMinY - 1e-1; + EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + } + + // Latitude > model max + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = modelMaxY + 1e-1; + EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + } + + // Before timeExtent.first + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + 1000, newLon, newLat, newZ)); + } + + // After timeExtent.last + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + 3000, newLon, newLat, newZ)); + } + + // Longitude < grid min + { + const double lon = gridMinX - 1e-1; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat, EPS); + EXPECT_EQ(newZ, zVal); + } + + // Longitude > grid max + { + const double lon = gridMaxX + 1e-1; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat, EPS); + EXPECT_EQ(newZ, zVal); + } + + // Latitude < grid min + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = gridMinY - 1e-1; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat, EPS); + EXPECT_EQ(newZ, zVal); + } + + // Latitude > grid max + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = gridMaxY + 1e-1; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat, EPS); + EXPECT_EQ(newZ, zVal); + } + + // Time function values to zero + { + const double lon = gridMinX + iQueriedX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + 2000, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat, EPS); + EXPECT_EQ(newZ, zVal); + } + + // Test vertical + j["components"][0]["displacement_type"] = "vertical"; + j["vertical_offset_unit"] = "metre"; + Evaluator<Grid, GridSet, EvaluatorIface> evalVertical( + MasterFile::parse(j.dump()), iface, 1, 1); + { + constexpr double alphaX = 0.25; + constexpr double alphaY = 0.125; + const double lon = gridMinX + iQueriedX * gridResX + alphaX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY; + EXPECT_TRUE(evalVertical.forward(iface, DegToRad(lon), DegToRad(lat), + zVal, tValid, newLon, newLat, newZ)); + EXPECT_NEAR(RadToDeg(newLon), lon, EPS); + EXPECT_NEAR(RadToDeg(newLat), lat, EPS); + + const double zBottom = + zOffsetQueriedXY + alphaX * (zOffsetQueriedXp1Y - zOffsetQueriedXY); + const double zTop = + zOffsetQueriedXYp1 + + alphaX * (zOffsetQueriedXp1Yp1 - zOffsetQueriedXYp1); + EXPECT_NEAR( + newZ, zVal + tFactor * (zBottom + alphaY * (zTop - zBottom)), EPS); + } + + // Test 3d + j["components"][0]["displacement_type"] = "3d"; + j["vertical_offset_unit"] = "metre"; + Evaluator<Grid, GridSet, EvaluatorIface> eval3d(MasterFile::parse(j.dump()), + iface, 1, 1); + { + constexpr double alphaX = 0.25; + constexpr double alphaY = 0.125; + const double lon = gridMinX + iQueriedX * gridResX + alphaX * gridResX; + const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY; + EXPECT_TRUE(eval3d.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + EXPECT_NEAR( + RadToDeg(newLon), + lon + + tFactor * (lonOffsetQueriedX + + alphaX * (lonOffsetQueriedXp1 - lonOffsetQueriedX)), + EPS); + EXPECT_NEAR( + RadToDeg(newLat), + lat + + tFactor * (latOffsetQueriedY + + alphaY * (latOffsetQueriedYp1 - latOffsetQueriedY)), + EPS); + + const double zBottom = + zOffsetQueriedXY + alphaX * (zOffsetQueriedXp1Y - zOffsetQueriedXY); + const double zTop = + zOffsetQueriedXYp1 + + alphaX * (zOffsetQueriedXp1Yp1 - zOffsetQueriedXYp1); + EXPECT_NEAR( + newZ, zVal + tFactor * (zBottom + alphaY * (zTop - zBottom)), EPS); + } +} + +// --------------------------------------------------------------------------- + +inline void DeltaLongLatToEastingNorthing(double phi, double dlam, double dphi, + double a, double b, double &de, + double &dn) { + const double sinphi = sin(phi); + const double cosphi = cos(phi); + const double a2 = a * a; + const double b2 = b * b; + const double X = a2 * (cosphi * cosphi) + b2 * (sinphi * sinphi); + const double sqrtX = sqrt(X); + de = dlam * (a2 * cosphi) / sqrtX; + dn = dphi * a2 * b2 / (sqrtX * X); +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluator_horizontal_unit_metre) { + + json j(getMinValidContent()); + j["horizontal_offset_method"] = "addition"; + j["horizontal_offset_unit"] = "metre"; + j["vertical_offset_unit"] = "metre"; + constexpr double tFactor = 0.5; + constexpr double gridMinX = 165.8; + constexpr double gridMinY = -37.5; + constexpr double gridMaxX = 166.2; + constexpr double gridMaxY = -37.2; + constexpr double gridResX = gridMaxX - gridMinX; + constexpr double gridResY = gridMaxY - gridMinY; + constexpr int extraPointX = 1; + constexpr int extraPointY = 1; + + j["components"] = { + {{"displacement_type", "horizontal"}, + {"uncertainty_type", "none"}, + {"extent", + {{"type", "bbox"}, + {"parameters", + {{"bbox", + {gridMinX - extraPointX * gridResX, + gridMinY - extraPointY * gridResY, gridMaxX, gridMaxY}}}}}}, + {"spatial_model", + { + {"type", "GeoTIFF"}, + {"interpolation_method", "XXXXXXX"}, + {"filename", "bla.tif"}, + }}, + {"time_function", + {{"type", "piecewise"}, + {"parameters", + {{"before_first", "zero"}, + {"after_last", "zero"}, + {"model", + {{{"epoch", "2010-01-01T00:00:00Z"}, {"scale_factor", tFactor}}, + {{"epoch", "2020-01-01T00:00:00Z"}, + {"scale_factor", tFactor}}}}}}}}}}; + + struct Grid : public GridPrototype { + bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset, + double &northingOffset) const { + ix -= extraPointX; + iy -= extraPointY; + if (ix == -1) + ix = 0; + if (iy == -1) + iy = 0; + if (ix == 0 && iy == 0) { + eastingOffset = 0.4f; + northingOffset = -0.2f; + } else if (ix == 1 && iy == 0) { + eastingOffset = 0.5f; + northingOffset = -0.25f; + } else if (ix == 0 && iy == 1) { + eastingOffset = 0.8f; + northingOffset = -0.4f; + } else if (ix == 1 && iy == 1) { + eastingOffset = 1.f; + northingOffset = -0.3f; + } else { + return false; + } + return true; + } + + bool getZOffset(int ix, int iy, double &zOffset) const { + ix -= extraPointX; + iy -= extraPointY; + if (ix == -1) + ix = 0; + if (iy == -1) + iy = 0; + if (ix == 0 && iy == 0) { + zOffset = 0.84f; + } else if (ix == 1 && iy == 0) { + zOffset = 0.75f; + } else if (ix == 0 && iy == 1) { + zOffset = 0.36f; + } else if (ix == 1 && iy == 1) { + zOffset = 0.f; + } else { + return false; + } + return true; + } + + bool getEastingNorthingZOffset(int ix, int iy, double &eastingOffset, + double &northingOffset, + double &zOffset) const { + return getEastingNorthingOffset(ix, iy, eastingOffset, + northingOffset) && + getZOffset(ix, iy, zOffset); + } + +#ifdef DEBUG_DEFMODEL + std::string name() const { return std::string(); } +#endif + }; + + struct GridSet : public GridSetPrototype<Grid> { + + Grid grid{}; + + GridSet() { + grid.minx = DegToRad(gridMinX - extraPointX * gridResX); + grid.miny = DegToRad(gridMinY - extraPointY * gridResY); + grid.resx = DegToRad(gridResX); + grid.resy = DegToRad(gridResY); + grid.width = 2 + extraPointX; + grid.height = 2 + extraPointY; + } + + const Grid *gridAt(double /*x */, double /* y */) { return &grid; } + }; + + struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> { + std::unique_ptr<GridSet> open(const std::string &filename) { + if (filename != "bla.tif") + return nullptr; + return std::unique_ptr<GridSet>(new GridSet()); + } + + bool isGeographicCRS(const std::string & /* crsDef */) { return true; } + +#ifdef DEBUG_DEFMODEL + void log(const std::string & /* msg */) {} +#endif + + void geographicToGeocentric(double lam, double phi, double height, + double a, double /*b*/, double es, + double &X, double &Y, double &Z) { + PJ_CONTEXT *ctx = proj_context_create(); + PJ *cart = proj_create( + ctx, + ("+proj=cart +a=" + osgeo::proj::internal::toString(a, 18) + + " +es=" + osgeo::proj::internal::toString(es, 18)) + .c_str()); + PJ_LPZ lpz; + lpz.lam = lam; + lpz.phi = phi; + lpz.z = height; + PJ_XYZ xyz = cart->fwd3d(lpz, cart); + X = xyz.x; + Y = xyz.y; + Z = xyz.z; + proj_destroy(cart); + proj_context_destroy(ctx); + } + + void geocentricToGeographic(double X, double Y, double Z, double a, + double /*b*/, double es, double &lam, + double &phi, double &height) { + PJ_CONTEXT *ctx = proj_context_create(); + PJ *cart = proj_create( + ctx, + ("+proj=cart +a=" + osgeo::proj::internal::toString(a, 18) + + " +es=" + osgeo::proj::internal::toString(es, 18)) + .c_str()); + PJ_XYZ xyz; + xyz.x = X; + xyz.y = Y; + xyz.z = Z; + PJ_LPZ lpz = cart->inv3d(xyz, cart); + lam = lpz.lam; + phi = lpz.phi; + height = lpz.z; + proj_destroy(cart); + proj_context_destroy(ctx); + } + }; + EvaluatorIface iface; + + constexpr double a = 6378137; + constexpr double b = 6356752.314140; + constexpr double tValid = 2018; + constexpr double zVal = 100; + + const struct { + double lon; + double lat; + double expected_de; + double expected_dn; + double expected_dz; + const char *displacement_type; + const char *interpolation_method; + } testPoints[] = { + {gridMinX - extraPointX * gridResX - 1e-11, + gridMinY - extraPointY * gridResY - 1e-11, 0.4, -0.2, 0, "horizontal", + "bilinear"}, + {gridMinX, gridMinY, 0.4, -0.2, 0, "horizontal", "bilinear"}, + {gridMaxX, gridMinY, 0.5, -0.25, 0, "horizontal", "bilinear"}, + {gridMinX, gridMaxY, 0.8, -0.4, 0, "horizontal", "bilinear"}, + {gridMaxX, gridMaxY, 1, -0.3, 0, "horizontal", "bilinear"}, + {gridMaxX + 1e-11, gridMaxY + 1e-11, 1, -0.3, 0, "horizontal", + "bilinear"}, + {165.9, -37.3, 0.70833334, -0.32083334, 0, "horizontal", "bilinear"}, + {165.9, -37.3, 0.70833334, -0.32083334, 0.4525, "3d", "bilinear"}, + + {gridMinX, gridMinY, 0.4, -0.2, 0, "horizontal", "geocentric_bilinear"}, + {gridMaxX, gridMinY, 0.5, -0.25, 0, "horizontal", + "geocentric_bilinear"}, + {gridMinX, gridMaxY, 0.8, -0.4, 0, "horizontal", "geocentric_bilinear"}, + {gridMaxX, gridMaxY, 1, -0.3, 0, "horizontal", "geocentric_bilinear"}, + {165.9, -37.3, 0.7083692044608846, -0.3209642339711405, 0, "horizontal", + "geocentric_bilinear"}, + {165.9, -37.3, 0.7083692044608846, -0.3209642339711405, 0.4525, "3d", + "geocentric_bilinear"}, + }; + + for (const auto &testPoint : testPoints) { + j["components"][0]["displacement_type"] = testPoint.displacement_type; + j["components"][0]["spatial_model"]["interpolation_method"] = + testPoint.interpolation_method; + Evaluator<Grid, GridSet, EvaluatorIface> eval( + MasterFile::parse(j.dump()), iface, a, b); + + const double lon = testPoint.lon; + const double lat = testPoint.lat; + double newLon; + double newLat; + double newZ; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)) + << lon << " " << lat << " " << testPoint.displacement_type + << testPoint.interpolation_method; + EXPECT_NEAR(newZ - zVal, tFactor * testPoint.expected_dz, 1e-8) + << lon << " " << lat << " " << testPoint.displacement_type + << testPoint.interpolation_method; + + double de; + double dn; + DeltaLongLatToEastingNorthing(DegToRad(lat), newLon - DegToRad(lon), + newLat - DegToRad(lat), a, b, de, dn); + EXPECT_NEAR(de, tFactor * testPoint.expected_de, 1e-8) + << lon << " " << lat << " " << testPoint.displacement_type + << testPoint.interpolation_method; + EXPECT_NEAR(dn, tFactor * testPoint.expected_dn, 1e-8) + << lon << " " << lat << " " << testPoint.displacement_type + << testPoint.interpolation_method; + + if (lon == gridMinX && lat == gridMinY) { + // Redo the exact same test, to test caching + double newLon2; + double newLat2; + double newZ2; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon2, newLat2, newZ2)); + EXPECT_EQ(newLon2, newLon); + EXPECT_EQ(newLat2, newLat); + EXPECT_EQ(newZ2, newZ); + + // Shift in longitude + EXPECT_TRUE(eval.forward(iface, DegToRad(lon - gridResX / 2), + DegToRad(lat), zVal, tValid, newLon2, + newLat2, newZ2)); + + // Redo test at original position + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon2, newLat2, newZ2)); + EXPECT_EQ(newLon2, newLon); + EXPECT_EQ(newLat2, newLat); + EXPECT_EQ(newZ2, newZ); + + // Shift in latitude + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), + DegToRad(lat - gridResY / 2), zVal, tValid, + newLon2, newLat2, newZ2)); + + // Redo test at original position + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon2, newLat2, newZ2)); + EXPECT_EQ(newLon2, newLon); + EXPECT_EQ(newLat2, newLat); + EXPECT_EQ(newZ2, newZ); + } + } + + // Test inverse() + { + j["horizontal_offset_method"] = "addition"; + j["components"][0]["displacement_type"] = "3d"; + j["components"][0]["spatial_model"]["interpolation_method"] = + "bilinear"; + Evaluator<Grid, GridSet, EvaluatorIface> eval( + MasterFile::parse(j.dump()), iface, a, b); + + const double lon = 165.9; + const double lat = -37.3; + double newLon; + double newLat; + double newZ; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + + double invLon; + double invLat; + double invZ; + EXPECT_TRUE(eval.inverse(iface, newLon, newLat, newZ, tValid, invLon, + invLat, invZ)); + EXPECT_NEAR(RadToDeg(invLon), lon, 1e-10); + EXPECT_NEAR(RadToDeg(invLat), lat, 1e-10); + EXPECT_NEAR(invZ, zVal, 1e-4); + } + + // Test horizontal_offset_method = geocentric + { + j["horizontal_offset_method"] = "geocentric"; + j["components"][0]["displacement_type"] = "3d"; + j["components"][0]["spatial_model"]["interpolation_method"] = + "bilinear"; + Evaluator<Grid, GridSet, EvaluatorIface> eval( + MasterFile::parse(j.dump()), iface, a, b); + + const double lon = gridMinX; + const double lat = gridMinY; + double newLon; + double newLat; + double newZ; + EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal, + tValid, newLon, newLat, newZ)); + + double de; + double dn; + constexpr double expected_de = 0.40000000948081327; + constexpr double expected_dn = -0.19999999810542682; + DeltaLongLatToEastingNorthing(DegToRad(lat), newLon - DegToRad(lon), + newLat - DegToRad(lat), a, b, de, dn); + EXPECT_NEAR(de, tFactor * expected_de, 1e-10); + EXPECT_NEAR(dn, tFactor * expected_dn, 1e-10); + EXPECT_NEAR(newZ - zVal, tFactor * 0.84, 1e-4); + } +} + +// --------------------------------------------------------------------------- + +TEST(defmodel, evaluator_projected_crs) { + + json j(getMinValidContent()); + j["horizontal_offset_method"] = "addition"; + j["horizontal_offset_unit"] = "metre"; + j["vertical_offset_unit"] = "metre"; + constexpr double gridMinX = 10000; + constexpr double gridMinY = 20000; + constexpr double gridMaxX = 30000; + constexpr double gridMaxY = 40000; + constexpr double gridResX = gridMaxX - gridMinX; + constexpr double gridResY = gridMaxY - gridMinY; + + j["extent"] + ["parameters"] = {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}; + j["components"] = { + {{"displacement_type", "horizontal"}, + {"uncertainty_type", "none"}, + {"extent", + {{"type", "bbox"}, + {"parameters", + {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}}}}, + {"spatial_model", + { + {"type", "GeoTIFF"}, + {"interpolation_method", "bilinear"}, + {"filename", "bla.tif"}, + }}, + {"time_function", {{"type", "constant"}}}}}; + + struct Grid : public GridPrototype { + bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset, + double &northingOffset) const { + if (ix == 0 && iy == 0) { + eastingOffset = 0.4; + northingOffset = -0.2; + } else if (ix == 1 && iy == 0) { + eastingOffset = 0.5; + northingOffset = -0.25; + } else if (ix == 0 && iy == 1) { + eastingOffset = 0.8; + northingOffset = -0.4; + } else if (ix == 1 && iy == 1) { + eastingOffset = 1.; + northingOffset = -0.3; + } else { + return false; + } + return true; + } + +#ifdef DEBUG_DEFMODEL + std::string name() const { return std::string(); } +#endif + }; + + struct GridSet : public GridSetPrototype<Grid> { + + Grid grid{}; + + GridSet() { + grid.minx = gridMinX; + grid.miny = gridMinY; + grid.resx = gridResX; + grid.resy = gridResY; + grid.width = 2; + grid.height = 2; + } + + const Grid *gridAt(double /*x */, double /* y */) { return &grid; } + }; + + struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> { + std::unique_ptr<GridSet> open(const std::string &filename) { + if (filename != "bla.tif") + return nullptr; + return std::unique_ptr<GridSet>(new GridSet()); + } + + bool isGeographicCRS(const std::string & /* crsDef */) { return false; } + +#ifdef DEBUG_DEFMODEL + void log(const std::string & /* msg */) {} +#endif + }; + EvaluatorIface iface; + + constexpr double a = 6378137; + constexpr double b = 6356752.314140; + constexpr double tValid = 2018; + constexpr double zVal = 100; + + Evaluator<Grid, GridSet, EvaluatorIface> eval(MasterFile::parse(j.dump()), + iface, a, b); + + double newX; + double newY; + double newZ; + EXPECT_TRUE(eval.forward(iface, gridMinX, gridMinY, zVal, tValid, newX, + newY, newZ)); + EXPECT_NEAR(newX - gridMinX, 0.4, 1e-8); + EXPECT_NEAR(newY - gridMinY, -0.2, 1e-8); + EXPECT_NEAR(newZ - zVal, 0, 1e-8); + + { + json jcopy(j); + jcopy["horizontal_offset_unit"] = "degree"; + EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>( + MasterFile::parse(jcopy.dump()), iface, a, b)), + EvaluatorException); + } + + { + json jcopy(j); + jcopy["horizontal_offset_method"] = "geocentric"; + EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>( + MasterFile::parse(jcopy.dump()), iface, a, b)), + EvaluatorException); + } + + { + json jcopy(j); + jcopy["components"][0]["spatial_model"]["interpolation_method"] = + "geocentric_bilinear"; + EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>( + MasterFile::parse(jcopy.dump()), iface, a, b)), + EvaluatorException); + } +} + +} // namespace + +#ifdef _MSC_VER +#pragma warning(pop) +#endif diff --git a/test/unit/test_grids.cpp b/test/unit/test_grids.cpp index 5240949e..e9841acd 100644 --- a/test/unit/test_grids.cpp +++ b/test/unit/test_grids.cpp @@ -66,7 +66,7 @@ TEST_F(GridTest, VerticalShiftGridSet_null) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 3); EXPECT_EQ(grid->height(), 3); - EXPECT_EQ(grid->extentAndRes().westLon, -M_PI); + EXPECT_EQ(grid->extentAndRes().west, -M_PI); EXPECT_TRUE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out = -1.0f; @@ -103,7 +103,7 @@ TEST_F(GridTest, HorizontalShiftGridSet_null) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 3); EXPECT_EQ(grid->height(), 3); - EXPECT_EQ(grid->extentAndRes().westLon, -M_PI); + EXPECT_EQ(grid->extentAndRes().west, -M_PI); EXPECT_TRUE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out1 = -1.0f; @@ -131,7 +131,7 @@ TEST_F(GridTest, HorizontalShiftGridSet_gtiff) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 4); EXPECT_EQ(grid->height(), 4); - EXPECT_EQ(grid->extentAndRes().westLon, 4.0 / 180 * M_PI); + EXPECT_EQ(grid->extentAndRes().west, 4.0 / 180 * M_PI); EXPECT_FALSE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out1 = -1.0f; @@ -154,7 +154,7 @@ TEST_F(GridTest, GenericShiftGridSet_null) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 3); EXPECT_EQ(grid->height(), 3); - EXPECT_EQ(grid->extentAndRes().westLon, -M_PI); + EXPECT_EQ(grid->extentAndRes().west, -M_PI); EXPECT_TRUE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out = -1.0f; @@ -180,7 +180,8 @@ TEST_F(GridTest, GenericShiftGridSet_gtiff) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 5); EXPECT_EQ(grid->height(), 5); - EXPECT_EQ(grid->extentAndRes().westLon, 21.0 / 180 * M_PI); + EXPECT_EQ(grid->extentAndRes().isGeographic, true); + EXPECT_EQ(grid->extentAndRes().west, 21.0 / 180 * M_PI); EXPECT_FALSE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out = -1.0f; @@ -224,4 +225,24 @@ TEST_F(GridTest, EXPECT_EQ(grid->height(), 8); } +// --------------------------------------------------------------------------- + +TEST_F(GridTest, GenericShiftGridSet_gtiff_projected) { + auto gridSet = NS_PROJ::GenericShiftGridSet::open( + m_ctxt, "tests/test_3d_grid_projected.tif"); + ASSERT_NE(gridSet, nullptr); + ASSERT_EQ(gridSet->gridAt(-1000, -1000), nullptr); + auto grid = gridSet->gridAt(1500300.0, 5400300.0); + ASSERT_NE(grid, nullptr); + EXPECT_EQ(grid->width(), 2); + EXPECT_EQ(grid->height(), 2); + EXPECT_EQ(grid->extentAndRes().isGeographic, false); + EXPECT_EQ(grid->extentAndRes().west, 1500000.0); + EXPECT_EQ(grid->extentAndRes().east, 1501000.0); + EXPECT_EQ(grid->extentAndRes().south, 5400000.0); + EXPECT_EQ(grid->extentAndRes().north, 5401000.0); + EXPECT_EQ(grid->extentAndRes().resX, 1000); + EXPECT_EQ(grid->extentAndRes().resY, 1000); +} + } // namespace |