diff --git a/.vscode/launch.json b/.vscode/launch.json index 0cc545d..e9a66f2 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -1,19 +1,29 @@ { "version": "0.2.0", - "configurations": [{ - "name": "Dart-proj4", - "program": "example/proj4dart_example.dart", - "request": "launch", - "type": "dart" - }, { - "name": "Dart-ogc-wkt", - "program": "example/proj4dart_ogc_wkt_example.dart", - "request": "launch", - "type": "dart" - }, { - "name": "Dart-esri-wkt", - "program": "example/proj4dart_esri_wkt_example.dart", - "request": "launch", - "type": "dart" - }] + "configurations": [ + { + "name": "Dart-proj4", + "program": "example/proj4dart_example.dart", + "request": "launch", + "type": "dart" + }, + { + "name": "Dart-ogc-wkt", + "program": "example/proj4dart_ogc_wkt_example.dart", + "request": "launch", + "type": "dart" + }, + { + "name": "Dart-esri-wkt", + "program": "example/proj4dart_esri_wkt_example.dart", + "request": "launch", + "type": "dart" + }, + { + "name": "Dart-shift-grid", + "program": "example/proj4dart_shift_grid_example.dart", + "request": "launch", + "type": "dart" + } + ] } \ No newline at end of file diff --git a/README.md b/README.md index d2ebb9e..709753c 100644 --- a/README.md +++ b/README.md @@ -99,6 +99,43 @@ var projection = Projection.parse(def); For full example visit [example/proj4dart_esri_wkt_example.dart](example/proj4dart_esri_wkt_example.dart) +#### Grid Based Datum Adjustments + +To use `+nadgrids=` in a proj definition, first read your NTv2 `.gsb` file into an `Uint8List`, then pass to `Projection.nadgrid`. E.g: + +Dart: +```dart +final bytes = await File( + 'assets/my_grid.gsb', +).readAsBytes(); + +Projection.nadgrid('key', bytes); +``` + +Flutter: +```dart +import 'package:flutter/services.dart' show rootBundle; +final bytes = (await rootBundle.load(fileName)).buffer.asUint8List(); + +Projection.nadgrid('key', bytes); +``` + +then use the given key in your definition, e.g. `+nadgrids=@key,null`. See [Proj4 General Parameters](http://proj.maptools.org/gen_parms.html). + +Nadgrid example: + +```dart +import 'package:flutter/services.dart' show rootBundle; + +final bytes = (await rootBundle.load('assets/nzgd2kgrid0005.gsb')).buffer.asUint8List(); +Projection.nadgrid('nzgd2kgrid0005.gsb', bytes); + +var def = '+proj=longlat +datum=nzgd49 +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 +nadgrids=nzgd2kgrid0005.gsb +no_defs'; +var namedProjection = Projection.add('EPSG:4272', def); +``` + +For full example visit [example/proj4dart_shift_grid_example.dart](example/proj4dart_shift_grid_example.dart) + ### Transform between Projections ```dart diff --git a/example/proj4dart_shift_grid_example.dart b/example/proj4dart_shift_grid_example.dart new file mode 100644 index 0000000..551970d --- /dev/null +++ b/example/proj4dart_shift_grid_example.dart @@ -0,0 +1,41 @@ +import 'dart:io'; + +import 'package:proj4dart/proj4dart.dart'; + +void main() async { + // read bytes from filesystem + // In case of Flutter assets: + // import 'package:flutter/services.dart' show rootBundle; + // final bytes = (await rootBundle.load(fileName)).buffer.asUint8List(); + final bytes = await File( + 'test/test_resources/ntv2_0_downsampled.gsb', + ).readAsBytes(); + + // load nadgrid + Projection.nadgrid('ntv2', bytes); + + // Define Point + var pointSrc = Point(x: -44.382211538462, y: 40.3768); + + // Define some projection which has '+nadgrids=' + var projSrc = Projection.add( + 'ntv2_from', + '+proj=longlat +ellps=clrk66 +nadgrids=@ignorable,ntv2,null', + ); + ; + + // Use built-in projection + var projDst = Projection.WGS84; + + // Forward transform (projected crs -> lonlat) + var pointForward = projSrc.transform(projDst, pointSrc); + print( + 'FORWARD: Transform point ${pointSrc.toArray()} from ntv2_from to WGS84: ${pointForward.toArray()}'); + // FORWARD: Transform point [-44.382211538462, 40.3768] from ntv2_from to WGS84: [-44.38074905319326, 40.37745745991217] + + // Inverse transform (lonlat -> projected crs) + var pointInverse = projDst.transform(projSrc, pointForward); + print( + 'INVERSE: Transform point ${pointForward.toArray()} from WGS84 to ntv2_from: ${pointInverse.toArray()}'); + // INVERSE: Transform point [-44.38074905319326, 40.37745745991217] from WGS84 to ntv2_from: [-44.38074905319326, 40.37745745991217] +} diff --git a/lib/src/classes/datum.dart b/lib/src/classes/datum.dart index 62eb9f2..a011eb2 100644 --- a/lib/src/classes/datum.dart +++ b/lib/src/classes/datum.dart @@ -1,3 +1,4 @@ +import 'package:proj4dart/src/classes/nadgrid.dart'; import 'package:proj4dart/src/constants/values.dart' as consts; class Datum { @@ -7,13 +8,15 @@ class Datum { final double b; final double es; final double ep2; + final List? grids; Datum(String? datumCode, List? datum_params, double a, double b, - double es, double ep2) + double es, double ep2, List? nadgrids) : a = a, b = b, es = es, - ep2 = ep2 { + ep2 = ep2, + grids = nadgrids { if (datumCode == null || datumCode == 'none') { datumType = consts.PJD_NODATUM; } else { @@ -37,5 +40,9 @@ class Datum { } } } + + if (nadgrids != null) { + datumType = consts.PJD_GRIDSHIFT; + } } } diff --git a/lib/src/classes/nadgrid.dart b/lib/src/classes/nadgrid.dart new file mode 100644 index 0000000..264395b --- /dev/null +++ b/lib/src/classes/nadgrid.dart @@ -0,0 +1,268 @@ +import 'dart:typed_data'; +import 'dart:math' as math; + +import 'package:proj4dart/src/globals/nadgrid_store.dart'; + +/// Resources for details of NTv2 file formats: +/// - https://web.archive.org/web/20140127204822if_/http://www.mgs.gov.on.ca:80/stdprodconsume/groups/content/@mgs/@iandit/documents/resourcelist/stel02_047447.pdf +/// - http://mimaka.com/help/gs/html/004_NTV2%20Data%20Format.htm +class NadgridParam { + final String name; + final bool mandatory; + final Nadgird? grid; + final bool isNull; + + NadgridParam({ + required this.name, + required this.mandatory, + this.grid, + required this.isNull, + }); +} + +class GridHeader { + final int nFields; + final int nSubgridFields; + final int nSubgrids; + final String shiftType; + final double fromSemiMajorAxis; + final double fromSemiMinorAxis; + final double toSemiMajorAxis; + final double toSemiMinorAxis; + + GridHeader({ + required this.nFields, + required this.nSubgridFields, + required this.nSubgrids, + required this.shiftType, + required this.fromSemiMajorAxis, + required this.fromSemiMinorAxis, + required this.toSemiMajorAxis, + required this.toSemiMinorAxis, + }); +} + +class SubGridHeader { + final String name; + final String parent; + final double lowerLatitude; + final double upperLatitude; + final double lowerLongitude; + final double upperLongitude; + final double latitudeInterval; + final double longitudeInterval; + final int gridNodeCount; + + SubGridHeader({ + required this.name, + required this.parent, + required this.lowerLatitude, + required this.upperLatitude, + required this.lowerLongitude, + required this.upperLongitude, + required this.latitudeInterval, + required this.longitudeInterval, + required this.gridNodeCount, + }); +} + +class GridNode { + final double latitudeShift; + final double longitudeShift; + final double latitudeAccuracy; + final double longitudeAccuracy; + + GridNode({ + required this.latitudeShift, + required this.longitudeShift, + required this.latitudeAccuracy, + required this.longitudeAccuracy, + }); +} + +class SubGrid { + final List ll; + final List del; + final List lim; + final int count; + final List> cvs; + + SubGrid({ + required this.ll, + required this.del, + required this.lim, + required this.count, + required this.cvs, + }); +} + +class Nadgird { + final GridHeader header; + final List subgrids; + + Nadgird({ + required this.header, + required this.subgrids, + }); + + /// Load a binary NTv2 file (.gsb) to a key that can be used in a proj string like +nadgrids=. Pass the NTv2 file + /// as an Uint8List. + factory Nadgird.factory(String key, Uint8List data) { + var view = ByteData.view(data.buffer); + var endian = detectEndian(view); + var header = readHeader(view, endian); + if (header.nSubgrids > 1) { + print( + 'Only single NTv2 subgrids are currently supported, subsequent sub grids are ignored'); + } + var subgrids = readSubgrids(view, header, endian); + var nadgrid = Nadgird(header: header, subgrids: subgrids); + return NadgridStore().register(key, nadgrid); + } + + /// Given a proj4 value for nadgrids, return an array of loaded grids + static List? getNadgrids(String? nadgrids) { + // Format details: http://proj.maptools.org/gen_parms.html + if (nadgrids == null) { + return null; + } + var grids = nadgrids.split(','); + return grids.map(parseNadgridString).toList(); + } + + static NadgridParam parseNadgridString(String value) { + // if (value.isEmpty) { + // return null; + // } + var optional = value[0] == '@'; + if (optional) { + value = value.substring(1); + } + if (value == 'null') { + return NadgridParam( + name: 'null', + mandatory: !optional, + grid: null, + isNull: true, + ); + } + return NadgridParam( + name: value, + mandatory: !optional, + grid: NadgridStore().get(value), + isNull: false, + ); + } + + static double secondsToRadians(double seconds) { + return (seconds / 3600) * math.pi / 180; + } + + static Endian detectEndian(ByteData view) { + return view.getUint8(8) == 11 ? Endian.little : Endian.big; + } + + static GridHeader readHeader(ByteData view, Endian endian) { + return GridHeader( + nFields: view.getInt32(8, endian), + nSubgridFields: view.getInt32(24, endian), + nSubgrids: view.getInt32(40, endian), + shiftType: decodeString(view, 56, 56 + 8).trim(), + fromSemiMajorAxis: view.getFloat64(120, endian), + fromSemiMinorAxis: view.getFloat64(136, endian), + toSemiMajorAxis: view.getFloat64(152, endian), + toSemiMinorAxis: view.getFloat64(168, endian), + ); + } + + static String decodeString(ByteData view, int start, int end) { + return String.fromCharCodes( + Iterable.generate(end - start, (x) => view.getUint8(start + x)), + ); + } + + static List readSubgrids( + ByteData view, GridHeader header, Endian endian) { + var gridOffset = 176; + var grids = []; + for (var i = 0; i < header.nSubgrids; i++) { + var subHeader = readGridHeader(view, gridOffset, endian); + var nodes = readGridNodes(view, gridOffset, subHeader, endian); + var lngColumnCount = (1 + + (subHeader.upperLongitude - subHeader.lowerLongitude) / + subHeader.longitudeInterval) + .round(); + var latColumnCount = (1 + + (subHeader.upperLatitude - subHeader.lowerLatitude) / + subHeader.latitudeInterval) + .round(); + + // Proj4 operates on radians whereas the coordinates are in seconds in the grid + grids.add( + SubGrid( + ll: [ + secondsToRadians(subHeader.lowerLongitude), + secondsToRadians(subHeader.lowerLatitude) + ], + del: [ + secondsToRadians(subHeader.longitudeInterval), + secondsToRadians(subHeader.latitudeInterval) + ], + lim: [lngColumnCount, latColumnCount], + count: subHeader.gridNodeCount, + cvs: mapNodes(nodes), + ), + ); + } + return grids; + } + + static List> mapNodes(List nodes) { + return nodes + .map( + (r) => [ + secondsToRadians(r.longitudeShift), + secondsToRadians(r.latitudeShift) + ], + ) + .toList(); + } + + static SubGridHeader readGridHeader( + ByteData view, int offset, Endian endian) { + return SubGridHeader( + name: decodeString(view, offset + 8, offset + 16).trim(), + parent: decodeString(view, offset + 24, offset + 24 + 8).trim(), + lowerLatitude: view.getFloat64(offset + 72, endian), + upperLatitude: view.getFloat64(offset + 88, endian), + lowerLongitude: view.getFloat64(offset + 104, endian), + upperLongitude: view.getFloat64(offset + 120, endian), + latitudeInterval: view.getFloat64(offset + 136, endian), + longitudeInterval: view.getFloat64(offset + 152, endian), + gridNodeCount: view.getInt32(offset + 168, endian), + ); + } + + static List readGridNodes( + ByteData view, int offset, SubGridHeader gridHeader, Endian endian) { + var nodesOffset = offset + 176; + var gridRecordLength = 16; + var gridShiftRecords = []; + + for (var i = 0; i < gridHeader.gridNodeCount; i++) { + var record = GridNode( + latitudeShift: + view.getFloat32(nodesOffset + i * gridRecordLength, endian), + longitudeShift: + view.getFloat32(nodesOffset + i * gridRecordLength + 4, endian), + latitudeAccuracy: + view.getFloat32(nodesOffset + i * gridRecordLength + 8, endian), + longitudeAccuracy: + view.getFloat32(nodesOffset + i * gridRecordLength + 12, endian), + ); + gridShiftRecords.add(record); + } + + return gridShiftRecords; + } +} diff --git a/lib/src/classes/proj_params.dart b/lib/src/classes/proj_params.dart index e91ef7c..1ef7cdf 100644 --- a/lib/src/classes/proj_params.dart +++ b/lib/src/classes/proj_params.dart @@ -1,4 +1,5 @@ import 'package:proj4dart/src/classes/datum.dart'; +import 'package:proj4dart/src/classes/nadgrid.dart'; import 'package:proj4dart/src/common/derive_constants.dart' as dc; import 'package:proj4dart/src/constants/datums.dart' as datums; import 'package:proj4dart/src/constants/prime_meridians.dart' as consts_pm; @@ -222,6 +223,7 @@ class ProjParams { var sphere = dc.sphere(a, b, rf, ellps!, map['sphere'] as bool?); var ecc = dc.eccentricity(sphere['a'] as double, sphere['b'] as double, sphere['rf'] as double?, R_A); + var parsedNadgrids = Nadgird.getNadgrids(nadgrids); map['a'] = sphere['a']; map['b'] = sphere['b']; map['rf'] = sphere['rf']; @@ -230,7 +232,8 @@ class ProjParams { map['e'] = ecc['e']; map['ep2'] = ecc['ep2']; if (datum == null) { - map['datum'] = Datum(datumCode, datum_params, a!, b!, es!, ep2!); + map['datum'] = + Datum(datumCode, datum_params, a!, b!, es!, ep2!, parsedNadgrids); } } diff --git a/lib/src/classes/projection.dart b/lib/src/classes/projection.dart index a9b3b40..87d7db2 100644 --- a/lib/src/classes/projection.dart +++ b/lib/src/classes/projection.dart @@ -1,5 +1,8 @@ +import 'dart:typed_data'; + import 'package:proj4dart/proj4dart.dart'; import 'package:proj4dart/src/classes/datum.dart'; +import 'package:proj4dart/src/classes/nadgrid.dart'; import 'package:proj4dart/src/classes/point.dart'; import 'package:proj4dart/src/classes/proj_params.dart'; import 'package:proj4dart/src/common/datum_transform.dart' as dt; @@ -105,6 +108,10 @@ abstract class Projection { return initializer(params); } + static void nadgrid(String key, Uint8List data) { + Nadgird.factory(key, data); + } + static final _mercatorCodes = ['3857', '900913', '3785', '102113']; /// Checks whether it is EPSG:3857 or one of its versions diff --git a/lib/src/common/datum_transform.dart b/lib/src/common/datum_transform.dart index 5e72ffa..4abc84a 100644 --- a/lib/src/common/datum_transform.dart +++ b/lib/src/common/datum_transform.dart @@ -1,7 +1,10 @@ import 'package:proj4dart/src/classes/datum.dart'; +import 'package:proj4dart/src/classes/nadgrid.dart'; import 'package:proj4dart/src/classes/point.dart'; import 'package:proj4dart/src/common/datum_utils.dart' as datum_utils; import 'package:proj4dart/src/constants/values.dart' as consts; +import 'dart:math' as math; +import 'package:proj4dart/src/common/utils.dart' as utils; bool checkParams(int type) { return (type == consts.PJD_3PARAM || type == consts.PJD_7PARAM); @@ -22,14 +25,33 @@ Point transform(Datum source, Datum dest, Point point) { } // If this datum requires grid shifts, then apply it to geodetic coordinates. + var source_a = source.a; + var source_es = source.es; + if (source.datumType == consts.PJD_GRIDSHIFT) { + applyGridShift(source, false, point); + source_a = consts.SRS_WGS84_SEMIMAJOR; + source_es = consts.SRS_WGS84_ESQUARED; + } + + var dest_a = dest.a; + var dest_b = dest.b; + var dest_es = dest.es; + if (dest.datumType == consts.PJD_GRIDSHIFT) { + dest_a = consts.SRS_WGS84_SEMIMAJOR; + dest_b = consts.SRS_WGS84_SEMIMINOR; + dest_es = consts.SRS_WGS84_ESQUARED; + } // Do we need to go through geocentric coordinates? - // if (source.es === dest.es && source.a === dest.a && !checkParams(source.datumType) && !checkParams(dest.datumType)) { - // return point; - // } + if (source_es == dest_es && + source_a == dest_a && + !checkParams(source.datumType) && + !checkParams(dest.datumType)) { + return point; + } // Convert to geocentric coordinates. - point = datum_utils.geodeticToGeocentric(point, source.es, source.a); + point = datum_utils.geodeticToGeocentric(point, source_es, source_a); // Convert between datums if (checkParams(source.datumType)) { point = datum_utils.geocentricToWgs84( @@ -39,5 +61,130 @@ Point transform(Datum source, Datum dest, Point point) { point = datum_utils.geocentricFromWgs84( point, dest.datumType, dest.datumParams); } - return datum_utils.geocentricToGeodetic(point, dest.es, dest.a, dest.b); + point = datum_utils.geocentricToGeodetic(point, dest_es, dest_a, dest_b); + + if (dest.datumType == consts.PJD_GRIDSHIFT) { + applyGridShift(dest, true, point); + } + + return point; +} + +void applyGridShift(Datum source, bool inverse, Point point) { + if (source.grids == null || source.grids!.isEmpty) { + throw Exception('Grid shift grids not found'); + } + var input = Point(x: -point.x, y: point.y); + var output = Point(x: double.nan, y: double.nan); + var onlyMandatoryGrids = false; + var attemptedGrids = []; + for (var i = 0; i < source.grids!.length; i++) { + var grid = source.grids![i]; + attemptedGrids.add(grid.name); + if (grid.isNull) { + output = input; + break; + } + onlyMandatoryGrids = grid.mandatory; + if (grid.grid == null) { + if (grid.mandatory) { + throw Exception("Unable to find mandatory grid '${grid.name}'"); + } + continue; + } + var subgrid = grid.grid!.subgrids[0]; + // skip tables that don't match our point at all + var epsilon = ((subgrid.del[1]).abs() + (subgrid.del[0]).abs()) / 10000.0; + var minX = subgrid.ll[0] - epsilon; + var minY = subgrid.ll[1] - epsilon; + var maxX = subgrid.ll[0] + (subgrid.lim[0] - 1) * subgrid.del[0] + epsilon; + var maxY = subgrid.ll[1] + (subgrid.lim[1] - 1) * subgrid.del[1] + epsilon; + if (minY > input.y || minX > input.x || maxY < input.y || maxX < input.x) { + continue; + } + output = applySubgridShift(input, inverse, subgrid); + if (!(output.x).isNaN) { + break; + } + } + if ((output.x).isNaN) { + throw Exception( + "Failed to find a grid shift table for location '${-input.x * consts.R2D} ${input.y * consts.R2D} tried: $attemptedGrids'", + ); + } + point.x = -output.x; + point.y = output.y; +} + +Point applySubgridShift(Point pin, bool inverse, SubGrid ct) { + var val = Point(x: double.nan, y: double.nan); + if (pin.x.isNaN) { + return val; + } + var tb = Point(x: pin.x, y: pin.y); + tb.x -= ct.ll[0]; + tb.y -= ct.ll[1]; + tb.x = utils.adjust_lon(tb.x - math.pi) + math.pi; + var t = nadInterpolate(tb, ct); + if (inverse) { + if (t.x.isNaN) { + return val; + } + t.x = tb.x - t.x; + t.y = tb.y - t.y; + var i = 9, tol = 1e-12; + Point dif, del; + do { + del = nadInterpolate(t, ct); + if (del.x.isNaN) { + print( + 'Inverse grid shift iteration failed, presumably at grid edge. Using first approximation.'); + break; + } + dif = Point(x: tb.x - (del.x + t.x), y: tb.y - (del.y + t.y)); + t.x += dif.x; + t.y += dif.y; + } while (i-- != 0 && dif.x.abs() > tol && dif.y.abs() > tol); + if (i < 0) { + print('Inverse grid shift iterator failed to converge.'); + return val; + } + val.x = utils.adjust_lon(t.x + ct.ll[0]); + val.y = t.y + ct.ll[1]; + } else { + if (!t.x.isNaN) { + val.x = pin.x + t.x; + val.y = pin.y + t.y; + } + } + return val; +} + +Point nadInterpolate(Point pin, SubGrid ct) { + var t = Point(x: pin.x / ct.del[0], y: pin.y / ct.del[1]); + var indx = Point(x: t.x.floorToDouble(), y: t.y.floorToDouble()); + var frct = Point(x: t.x - 1.0 * indx.x, y: t.y - 1.0 * indx.y); + var val = Point(x: double.nan, y: double.nan); + int inx; + if (indx.x < 0 || indx.x >= ct.lim[0]) { + return val; + } + if (indx.y < 0 || indx.y >= ct.lim[1]) { + return val; + } + inx = ((indx.y * ct.lim[0]) + indx.x).toInt(); + var f00 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); + inx++; + var f10 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); + inx += ct.lim[0]; + var f11 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); + inx--; + var f01 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); + var m11 = frct.x * frct.y, + m10 = frct.x * (1.0 - frct.y), + m00 = (1.0 - frct.x) * (1.0 - frct.y), + m01 = (1.0 - frct.x) * frct.y; + val.x = (m00 * f00.x + m10 * f10.x + m01 * f01.x + m11 * f11.x); + val.y = (m00 * f00.y + m10 * f10.y + m01 * f01.y + m11 * f11.y); + return val; } diff --git a/lib/src/common/utils.dart b/lib/src/common/utils.dart index 91e7d81..4702cc7 100644 --- a/lib/src/common/utils.dart +++ b/lib/src/common/utils.dart @@ -178,33 +178,6 @@ double imlfn(double ml, double e0, double e1, double e2, double e3) { return double.nan; } -Point inverseNadCvt(Point t, Point val, tb, ct) { - if ((t.x).isNaN) { - return val; - } - t.x = tb.x + t.x; - t.y = tb.y - t.y; - var i = 9; - var tol = 1e-12; - Point dif; - Point del; - do { - del = nad_intr(t, ct); - if ((del.x).isNaN) { - break; - } - dif = Point(x: t.x - del.x - tb.x, y: t.y + del.y - tb.y); - t.x -= dif.x; - t.y -= dif.y; - } while (i-- != 0 && (dif.x).abs() > tol && (dif.y).abs() > tol); - if (i < 0) { - return val; - } - val.x = adjust_lon(t.x + ct.ll[0]); - val.y = t.y + ct.ll[1]; - return val; -} - double invlatiso(double eccent, double ts) { var phi = fL(1, ts); var Iphi = 0.0; @@ -286,90 +259,6 @@ double msfnz(double eccent, double sinphi, double cosphi) { return cosphi / (math.sqrt(1 - con * con)); } -Point nad_cvt(Point pin, bool inverse, ct) { - var val = Point(x: double.nan, y: double.nan); - if (pin.x.isNaN) { - return val; - } - var tb = Point(x: pin.x, y: pin.y); - tb.x -= ct.ll[0]; - tb.y -= ct.ll[1]; - tb.x = adjust_lon(tb.x - math.pi) + math.pi; - var t = nad_intr(tb, ct); - if (inverse) { - return inverseNadCvt(t, val, tb, ct); - } else { - if (!t.x.isNaN) { - val.x = pin.x - t.x; - val.y = pin.y + t.y; - } - } - return val; -} - -Point nad_intr(pin, ct) { - // force computation by decreasing by 1e-7 to be as closed as possible - // from computation under C:C++ by leveraging rounding problems ... - var t = Point(x: (pin.x - 1e-7) / ct.del[0], y: (pin.y - 1e-7) / ct.del[1]); - var indx = Point(x: (t.x).floorToDouble(), y: (t.y).floorToDouble()); - Point frct = Point(x: t.x - 1 * indx.x, y: t.y - 1 * indx.y); - var val = Point(x: double.nan, y: double.nan); - var temp = nadInterBreakout(indx, frct, 'x', 0, ct); - if (temp) { - indx = temp[0]; - frct = temp[1]; - } else { - return val; - } - temp = nadInterBreakout(indx, frct, 'y', 1, ct); - if (temp != false) { - indx = temp[0]; - frct = temp[1]; - } else { - return val; - } - var inx = (indx.y * ct.lim[0]) + indx.x; - var f00 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); - inx++; - var f10 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); - inx += ct.lim[0]; - var f11 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); - inx--; - var f01 = Point(x: ct.cvs[inx][0], y: ct.cvs[inx][1]); - var m11 = frct.x * frct.y, - m10 = frct.x * (1 - frct.y), - m00 = (1 - frct.x) * (1 - frct.y), - m01 = (1 - frct.x) * frct.y; - val.x = (m00 * f00.x + m10 * f10.x + m01 * f01.x + m11 * f11.x); - val.y = (m00 * f00.y + m10 * f10.y + m01 * f01.y + m11 * f11.y); - return val; -} - -dynamic nadInterBreakout(indx, frct, String letter, int number, ct) { - var inx; - if (indx[letter] < 0) { - if (!(indx[letter] == -1 && frct[letter] > 0.99999999999)) { - return false; - } - indx[letter]++; - frct[letter] = 0; - } else { - inx = indx[letter] + 1; - if (inx >= ct.lim[number]) { - if (!(inx == ct.lim[number] && frct[letter] < 1e-11)) { - return false; - } - if (letter == 'x') { - indx[letter]--; - } else { - indx[letter]++; - } - frct[letter] = 1; - } - } - return [indx, frct]; -} - double phi2z(double eccent, double ts) { var eccnth = 0.5 * eccent; double con, dphi; @@ -504,7 +393,8 @@ void checkSanity(Point point) { Point adjust_axis(Projection crs, bool denorm, Point point) { var xin = point.x, yin = point.y, zin = point.z ?? 0.0; var v, t, i; - var pointString = ''' + var pointString = + ''' { "x": ${point.x}, "y": ${point.y}, @@ -512,7 +402,8 @@ Point adjust_axis(Projection crs, bool denorm, Point point) { } '''; var pointObj = jsonDecode(pointString); - var outString = ''' + var outString = + ''' { "x": null, "y": null, diff --git a/lib/src/constants/values.dart b/lib/src/constants/values.dart index 3a32f95..49792d1 100644 --- a/lib/src/constants/values.dart +++ b/lib/src/constants/values.dart @@ -1,27 +1,32 @@ import 'dart:math' as math; -final PJD_3PARAM = 1; -final PJD_7PARAM = 2; -final PJD_WGS84 = 4; // WGS84 or equivalent -final PJD_NODATUM = 5; // WGS84 or equivalent -final SEC_TO_RAD = 4.84813681109535993589914102357e-6; -final HALF_PI = math.pi / 2; +const PJD_3PARAM = 1; +const PJD_7PARAM = 2; +const PJD_GRIDSHIFT = 3; +const PJD_WGS84 = 4; // WGS84 or equivalent +const PJD_NODATUM = 5; // WGS84 or equivalent +const SRS_WGS84_SEMIMAJOR = 6378137.0; // only used in grid shift transforms +const SRS_WGS84_SEMIMINOR = 6356752.314; // only used in grid shift transforms +const SRS_WGS84_ESQUARED = + 0.0066943799901413165; // only used in grid shift transforms +const SEC_TO_RAD = 4.84813681109535993589914102357e-6; +const HALF_PI = math.pi / 2; // ellipoid pj_set_ell.c -final SIXTH = 0.1666666666666666667; +const SIXTH = 0.1666666666666666667; // 1/6 -final RA4 = 0.04722222222222222222; +const RA4 = 0.04722222222222222222; // 17/360 -final RA6 = 0.02215608465608465608; -final EPSLN = 1.0e-10; +const RA6 = 0.02215608465608465608; +const EPSLN = 1.0e-10; // you'd think you could use Number.EPSILON above but that makes // Mollweide get into an infinate loop. -final D2R = 0.01745329251994329577; -final R2D = 57.29577951308232088; -final FORTPI = math.pi / 4; -final TWO_PI = math.pi * 2; +const D2R = 0.01745329251994329577; +const R2D = 57.29577951308232088; +const FORTPI = math.pi / 4; +const TWO_PI = math.pi * 2; // SPI is slightly greater than Math.PI, so values that exceed the -180..180 // degree range by a tiny amount don't get wrapped. This prevents points that // have drifted from their original location along the 180th meridian (due to // floating point error) from changing their sign. -final SPI = 3.14159265359; +const SPI = 3.14159265359; diff --git a/lib/src/globals/nadgrid_store.dart b/lib/src/globals/nadgrid_store.dart new file mode 100644 index 0000000..d789032 --- /dev/null +++ b/lib/src/globals/nadgrid_store.dart @@ -0,0 +1,34 @@ +import 'package:meta/meta.dart'; +import 'package:proj4dart/src/classes/nadgrid.dart'; + +/// Global class for storing Nadgrids +class NadgridStore { + final Map _nadgridCache = {}; + + bool get isEmpty => _nadgridCache.isEmpty; + + /// Clear Nadgird Cache only for testing purpose + @visibleForTesting + void clearNadgirdCache() { + _nadgridCache.clear(); + } + + static final NadgridStore _nadgridStore = NadgridStore._internal(); + + factory NadgridStore() { + return _nadgridStore; + } + + /// Private constructor + NadgridStore._internal(); + + /// Get Nadgird from the store + Nadgird? get(String gridName) { + return _nadgridCache[gridName]; + } + + /// Registers Nadgird to the store + Nadgird register(String gridName, Nadgird nadgrid) { + return _nadgridCache[gridName] = nadgrid; + } +} diff --git a/test/data/all_proj4_defs.dart b/test/data/all_proj4_defs.dart index 9405e6e..c98154c 100644 --- a/test/data/all_proj4_defs.dart +++ b/test/data/all_proj4_defs.dart @@ -60,6 +60,39 @@ const blackList = { 'ESRI:102497': 'undefined', 'ESRI:102498': 'undefined', 'ESRI:102590': 'undefined', + 'EPSG:4272': 'Failed to find a grid shift table for location', + 'EPSG:27205': 'Failed to find a grid shift table for location', + 'EPSG:27206': 'Failed to find a grid shift table for location', + 'EPSG:27207': 'Failed to find a grid shift table for location', + 'EPSG:27208': 'Failed to find a grid shift table for location', + 'EPSG:27209': 'Failed to find a grid shift table for location', + 'EPSG:27210': 'Failed to find a grid shift table for location', + 'EPSG:27211': 'Failed to find a grid shift table for location', + 'EPSG:27212': 'Failed to find a grid shift table for location', + 'EPSG:27214': 'Failed to find a grid shift table for location', + 'EPSG:27215': 'Failed to find a grid shift table for location', + 'EPSG:27216': 'Failed to find a grid shift table for location', + 'EPSG:27217': 'Failed to find a grid shift table for location', + 'EPSG:27218': 'Failed to find a grid shift table for location', + 'EPSG:27219': 'Failed to find a grid shift table for location', + 'EPSG:27220': 'Failed to find a grid shift table for location', + 'EPSG:27221': 'Failed to find a grid shift table for location', + 'EPSG:27222': 'Failed to find a grid shift table for location', + 'EPSG:27223': 'Failed to find a grid shift table for location', + 'EPSG:27224': 'Failed to find a grid shift table for location', + 'EPSG:27225': 'Failed to find a grid shift table for location', + 'EPSG:27226': 'Failed to find a grid shift table for location', + 'EPSG:27227': 'Failed to find a grid shift table for location', + 'EPSG:27228': 'Failed to find a grid shift table for location', + 'EPSG:27229': 'Failed to find a grid shift table for location', + 'EPSG:27230': 'Failed to find a grid shift table for location', + 'EPSG:27231': 'Failed to find a grid shift table for location', + 'EPSG:27232': 'Failed to find a grid shift table for location', + 'EPSG:27258': 'Failed to find a grid shift table for location', + 'EPSG:27259': 'Failed to find a grid shift table for location', + 'EPSG:27260': 'Failed to find a grid shift table for location', + 'EPSG:27291': 'Failed to find a grid shift table for location', + 'EPSG:27292': 'Failed to find a grid shift table for location', }; /// proj4 definitions based on PostGIS 3.0.1 (8500 definitions) diff --git a/test/proj4dart_test.dart b/test/proj4dart_test.dart index e8ee98f..7f39391 100644 --- a/test/proj4dart_test.dart +++ b/test/proj4dart_test.dart @@ -1,4 +1,5 @@ import 'dart:collection'; +import 'dart:io'; import 'package:proj4dart/proj4dart.dart'; import 'package:proj4dart/src/globals/projection_store.dart'; @@ -16,6 +17,101 @@ import './results/all_proj4_esri_wkt_results.dart' as all_proj4_esri_results; import 'classes/close_to_helper.dart'; void main() { + setUpAll(() async { + final bytes = await File( + 'test/test_resources/nzgd2kgrid0005.gsb', + ).readAsBytes(); + + Projection.nadgrid('nzgd2kgrid0005.gsb', bytes); + }); + + group('Nadgrid tests', () { + setUpAll(() async { + final bytes = await File( + 'test/test_resources/ntv2_0_downsampled.gsb', + ).readAsBytes(); + + Projection.nadgrid('ntv2', bytes); + }); + + setUp(() => ProjectionStore().clearProjectionCache()); + + test('Create all projections via proj4 def strings and find all of them', + () async { + var tests = >[ + [ + -44.382211538462, + 40.3768, + -44.380749, + 40.377457 + ], // just inside the lower limit + [ + -87.617788, + 59.623262, + -87.617659, + 59.623441 + ], // just inside the upper limit + [-44.5, 40.5, -44.498553, 40.500632], // inside the first square + [ + -60, + 50, + -59.999192, + 50.000058 + ], // a general point towards the middle of the grid + [0, 0, 0, 0] // fall back to null + ]; + + final from = Projection.add('ntv2_from', + '+proj=longlat +ellps=clrk66 +nadgrids=@ignorable,ntv2,null'); + final to = + Projection.add('ntv2_to', '+proj=longlat +datum=WGS84 +no_defs'); + var converter = ProjectionTuple(fromProj: from, toProj: to); + + tests.forEach((test) { + var fromLng = test[0]; + var fromLat = test[1]; + var toLng = test[2]; + var toLat = test[3]; + + var actual = converter.forward(Point(x: fromLng, y: fromLat)); + + expect( + actual.x, + closeTo(toLng, 0.000001), + ); + + expect( + actual.y, + closeTo(toLat, 0.000001), + ); + }); + + var inverseTests = >[ + [-44.5, 40.5, -44.498553, 40.500632], + [-60, 50, -59.999192, 50.000058] + ]; + + inverseTests.forEach((test) { + var fromLng = test[0]; + var fromLat = test[1]; + var toLng = test[2]; + var toLat = test[3]; + + var actual = converter.inverse(Point(x: toLng, y: toLat)); + + expect( + actual.x, + closeTo(fromLng, 0.000001), + ); + + expect( + actual.y, + closeTo(fromLat, 0.000001), + ); + }); + }); + }); + group('Bulk tests', () { setUp(() => ProjectionStore().clearProjectionCache()); diff --git a/test/test_resources/ntv2_0_downsampled.gsb b/test/test_resources/ntv2_0_downsampled.gsb new file mode 100644 index 0000000..06b8111 Binary files /dev/null and b/test/test_resources/ntv2_0_downsampled.gsb differ diff --git a/test/test_resources/nzgd2kgrid0005.gsb b/test/test_resources/nzgd2kgrid0005.gsb new file mode 100644 index 0000000..840bd00 Binary files /dev/null and b/test/test_resources/nzgd2kgrid0005.gsb differ