import { PJD_3PARAM, PJD_7PARAM, PJD_GRIDSHIFT, PJD_NODATUM, R2D, SRS_WGS84_ESQUARED, SRS_WGS84_SEMIMAJOR, SRS_WGS84_SEMIMINOR } from './constants/values'; import {geodeticToGeocentric, geocentricToGeodetic, geocentricToWgs84, geocentricFromWgs84, compareDatums} from './datumUtils'; import adjust_lon from "./common/adjust_lon"; function checkParams(type) { return (type === PJD_3PARAM || type === PJD_7PARAM); } export default function(source, dest, point) { // Short cut if the datums are identical. if (compareDatums(source, dest)) { return point; // in this case, zero is sucess, // whereas cs_compare_datums returns 1 to indicate TRUE // confusing, should fix this } // Explicitly skip datum transform by setting 'datum=none' as parameter for either source or dest if (source.datum_type === PJD_NODATUM || dest.datum_type === PJD_NODATUM) { return 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.datum_type === PJD_GRIDSHIFT) { var gridShiftCode = applyGridShift(source, false, point); if (gridShiftCode !== 0) { return undefined; } source_a = SRS_WGS84_SEMIMAJOR; source_es = SRS_WGS84_ESQUARED; } var dest_a = dest.a; var dest_b = dest.b; var dest_es = dest.es; if (dest.datum_type === PJD_GRIDSHIFT) { dest_a = SRS_WGS84_SEMIMAJOR; dest_b = SRS_WGS84_SEMIMINOR; dest_es = SRS_WGS84_ESQUARED; } // Do we need to go through geocentric coordinates? if (source_es === dest_es && source_a === dest_a && !checkParams(source.datum_type) && !checkParams(dest.datum_type)) { return point; } // Convert to geocentric coordinates. point = geodeticToGeocentric(point, source_es, source_a); // Convert between datums if (checkParams(source.datum_type)) { point = geocentricToWgs84(point, source.datum_type, source.datum_params); } if (checkParams(dest.datum_type)) { point = geocentricFromWgs84(point, dest.datum_type, dest.datum_params); } point = geocentricToGeodetic(point, dest_es, dest_a, dest_b); if (dest.datum_type === PJD_GRIDSHIFT) { var destGridShiftResult = applyGridShift(dest, true, point); if (destGridShiftResult !== 0) { return undefined; } } return point; } export function applyGridShift(source, inverse, point) { if (source.grids === null || source.grids.length === 0) { console.log('Grid shift grids not found'); return -1; } var input = {x: -point.x, y: point.y}; var output = {x: Number.NaN, y: Number.NaN}; var onlyMandatoryGrids = false; var attemptedGrids = []; for (var i = 0; i < source.grids.length; i++) { var grid = source.grids[i]; attemptedGrids.push(grid.name); if (grid.isNull) { output = input; break; } onlyMandatoryGrids = grid.mandatory; if (grid.grid === null) { if (grid.mandatory) { console.log("Unable to find mandatory grid '" + grid.name + "'"); return -1; } continue; } var subgrid = grid.grid.subgrids[0]; // skip tables that don't match our point at all var epsilon = (Math.abs(subgrid.del[1]) + Math.abs(subgrid.del[0])) / 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 (!isNaN(output.x)) { break; } } if (isNaN(output.x)) { console.log("Failed to find a grid shift table for location '"+ -input.x * R2D + " " + input.y * R2D + " tried: '" + attemptedGrids + "'"); return -1; } point.x = -output.x; point.y = output.y; return 0; } function applySubgridShift(pin, inverse, ct) { var val = {x: Number.NaN, y: Number.NaN}; if (isNaN(pin.x)) { return val; } var tb = {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 = nadInterpolate(tb, ct); if (inverse) { if (isNaN(t.x)) { return val; } t.x = tb.x - t.x; t.y = tb.y - t.y; var i = 9, tol = 1e-12; var dif, del; do { del = nadInterpolate(t, ct); if (isNaN(del.x)) { console.log("Inverse grid shift iteration failed, presumably at grid edge. Using first approximation."); break; } dif = {x: tb.x - (del.x + t.x), y: tb.y - (del.y + t.y)}; t.x += dif.x; t.y += dif.y; } while (i-- && Math.abs(dif.x) > tol && Math.abs(dif.y) > tol); if (i < 0) { console.log("Inverse grid shift iterator failed to converge."); return val; } val.x = adjust_lon(t.x + ct.ll[0]); val.y = t.y + ct.ll[1]; } else { if (!isNaN(t.x)) { val.x = pin.x + t.x; val.y = pin.y + t.y; } } return val; } function nadInterpolate(pin, ct) { var t = {x: pin.x / ct.del[0], y: pin.y / ct.del[1]}; var indx = {x: Math.floor(t.x), y: Math.floor(t.y)}; var frct = {x: t.x - 1.0 * indx.x, y: t.y - 1.0 * indx.y}; var val= {x: Number.NaN, y: Number.NaN}; var 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; var f00 = {x: ct.cvs[inx][0], y: ct.cvs[inx][1]}; inx++; var f10= {x: ct.cvs[inx][0], y: ct.cvs[inx][1]}; inx += ct.lim[0]; var f11 = {x: ct.cvs[inx][0], y: ct.cvs[inx][1]}; inx--; var f01 = {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; }