diff options
Diffstat (limited to 'app/utm.js')
| -rw-r--r-- | app/utm.js | 240 |
1 files changed, 240 insertions, 0 deletions
diff --git a/app/utm.js b/app/utm.js new file mode 100644 index 0000000..1a8790f --- /dev/null +++ b/app/utm.js @@ -0,0 +1,240 @@ +/** + * UTM (Universal Transverse Mercator) coordinate conversion utilities + * Converts WGS84 coordinates (latitude/longitude) to UTM coordinates + */ +export class UTMConverter { + // WGS84 ellipsoid parameters + static WGS84_A = 6378137.0; // semi-major axis in meters + static WGS84_F = 1 / 298.257223563; // flattening + static WGS84_E = 0.081819191; // eccentricity + static WGS84_ESQ = 0.00669438; // eccentricity squared + + // UTM parameters + static UTM_SCALE_FACTOR = 0.9996; + static UTM_FALSE_EASTING = 500000; // in meters + static UTM_FALSE_NORTHING_N = 0; // for northern hemisphere + static UGM_FALSE_NORTHING_S = 10000000; // for southern hemisphere + + /** + * Convert WGS84 coordinates to UTM + * @param {number} longitude - Longitude in decimal degrees + * @param {number} latitude - Latitude in decimal degrees + * @returns {Object} UTM coordinates with zone information + */ + static toUTM(longitude, latitude) { + // Validate input coordinates + if (longitude < -180 || longitude > 180) { + throw new Error(`Longitude ${longitude} is out of valid range [-180, 180]`); + } + if (latitude < -90 || latitude > 90) { + throw new Error(`Latitude ${latitude} is out of valid range [-90, 90]`); + } + + const zone = UTMConverter.getUTMZone(longitude, latitude); + const hemisphere = latitude >= 0 ? "N" : "S"; + + // Convert to radians + const latRad = (latitude * Math.PI) / 180; + const lonRad = (longitude * Math.PI) / 180; + + // Central meridian for the zone (in radians) + const lonOrigin = (((zone - 1) * 6 - 180 + 3) * Math.PI) / 180; + const dLon = lonRad - lonOrigin; + + // Ellipsoid constants + const e = UTMConverter.WGS84_E; + const e2 = UTMConverter.WGS84_ESQ; + const e4 = e2 * e2; + const e6 = e4 * e2; + + const N = UTMConverter.WGS84_A / Math.sqrt(1 - e2 * Math.sin(latRad) * Math.sin(latRad)); + const T = Math.tan(latRad) * Math.tan(latRad); + const C = (e2 * Math.cos(latRad) * Math.cos(latRad)) / (1 - e2); + const A = Math.cos(latRad) * dLon; + + const M = + UTMConverter.WGS84_A * + ((1 - e2 / 4 - (3 * e4) / 64 - (5 * e6) / 256) * latRad - + ((3 * e2) / 8 + (3 * e4) / 32 + (45 * e6) / 1024) * Math.sin(2 * latRad) + + ((15 * e4) / 256 + (45 * e6) / 1024) * Math.sin(4 * latRad) - + ((35 * e6) / 3072) * Math.sin(6 * latRad)); + + // Calculate easting + const easting = + UTMConverter.UTM_SCALE_FACTOR * + N * + (A + + ((1 - T + C) * A * A * A) / 6 + + ((5 - 18 * T + T * T + 72 * C - 58 * e2) * A * A * A * A * A) / 120) + + UTMConverter.UTM_FALSE_EASTING; + + // Calculate northing + const northing = + UTMConverter.UTM_SCALE_FACTOR * + (M + + N * + Math.tan(latRad) * + ((A * A) / 2 + + ((5 - T + 9 * C + 4 * C * C) * A * A * A * A) / 24 + + ((61 - 58 * T + T * T + 600 * C - 330 * e2) * A * A * A * A * A * A) / 720)); + + // Adjust northing for southern hemisphere + const adjustedNorthing = + hemisphere === "S" ? northing + UTMConverter.UGM_FALSE_NORTHING_S : northing; + + return { + datum: "WGS84", + easting: Math.round(easting * 100) / 100, // Round to 2 decimal places + hemisphere: hemisphere, + northing: Math.round(adjustedNorthing * 100) / 100, + zone: zone, + }; + } + + /** + * Determine UTM zone from longitude and latitude + * @param {number} longitude - Longitude in decimal degrees + * @param {number} latitude - Latitude in decimal degrees + * @returns {number} UTM zone number + */ + static getUTMZone(longitude, latitude) { + // Special zones for Norway and Svalbard + if (latitude >= 56.0 && latitude < 64.0 && longitude >= 3.0 && longitude < 12.0) { + return 32; + } + + // Special zones for Svalbard + if (latitude >= 72.0 && latitude < 84.0) { + if (longitude >= 0.0 && longitude < 9.0) return 31; + if (longitude >= 9.0 && longitude < 21.0) return 33; + if (longitude >= 21.0 && longitude < 33.0) return 35; + if (longitude >= 33.0 && longitude < 42.0) return 37; + } + + // Standard zone calculation + let zone = Math.floor((longitude + 180) / 6) + 1; + + // Handle zones 31-37 for northern Europe special cases + if (latitude >= 72.0 && latitude < 84.0 && longitude >= 0.0 && longitude < 42.0) { + if (longitude < 9.0) zone = 31; + else if (longitude < 21.0) zone = 33; + else if (longitude < 33.0) zone = 35; + else zone = 37; + } + + return Math.min(Math.max(zone, 1), 60); + } + + /** + * Convert UTM coordinates back to WGS84 + * @param {number} easting - UTM easting in meters + * @param {number} northing - UTM northing in meters + * @param {number} zone - UTM zone number + * @param {string} hemisphere - 'N' for northern, 'S' for southern + * @returns {Object} WGS84 coordinates {longitude, latitude} + */ + static fromUTM(easting, northing, zone, hemisphere = "N") { + if (zone < 1 || zone > 60) { + throw new Error(`UTM zone ${zone} is out of valid range [1, 60]`); + } + + // Adjust northing for southern hemisphere + const adjustedNorthing = + hemisphere === "S" ? northing - UTMConverter.UGM_FALSE_NORTHING_S : northing; + + const e = UTMConverter.WGS84_E; + const e2 = UTMConverter.WGS84_ESQ; + const e4 = e2 * e2; + const e6 = e4 * e2; + + // Remove false easting and scale factor + const x = (easting - UTMConverter.UTM_FALSE_EASTING) / UTMConverter.UTM_SCALE_FACTOR; + const y = adjustedNorthing / UTMConverter.UTM_SCALE_FACTOR; + + // Footprint latitude + const M = y; + const mu = M / (UTMConverter.WGS84_A * (1 - e2 / 4 - (3 * e4) / 64 - (5 * e6) / 256)); + + const e1 = (1 - Math.sqrt(1 - e2)) / (1 + Math.sqrt(1 - e2)); + const fp = + mu + + ((3 * e1) / 2 - (27 * e1 * e1 * e1) / 32) * Math.sin(2 * mu) + + ((21 * e1 * e1) / 16 - (55 * e1 * e1 * e1 * e1) / 32) * Math.sin(4 * mu) + + ((151 * e1 * e1 * e1) / 96) * Math.sin(6 * mu); + + // Calculate latitude and longitude + const C1 = (e2 * Math.cos(fp) * Math.cos(fp)) / (1 - e2); + const T1 = Math.tan(fp) * Math.tan(fp); + const N1 = UTMConverter.WGS84_A / Math.sqrt(1 - e2 * Math.sin(fp) * Math.sin(fp)); + const R1 = (UTMConverter.WGS84_A * (1 - e2)) / (1 - e2 * Math.sin(fp) * Math.sin(fp)) ** 1.5; + const D = x / N1; + + const lat = + fp - + ((N1 * Math.tan(fp)) / R1) * + ((D * D) / 2 - + ((5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * e2) * D * D * D * D) / 24 + + ((61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * e2 - 3 * C1 * C1) * + D * + D * + D * + D * + D * + D) / + 720); + + const lonOrigin = (((zone - 1) * 6 - 180 + 3) * Math.PI) / 180; + const lon = + lonOrigin + + (D - + ((1 + 2 * T1 + C1) * D * D * D) / 6 + + ((5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * e2 + 24 * T1 * T1) * D * D * D * D * D) / 120) / + Math.cos(fp); + + return { + latitude: (lat * 180) / Math.PI, + longitude: (lon * 180) / Math.PI, + }; + } + + /** + * Convert Point geometry to UTM coordinates + * @param {Point} point - Point geometry + * @returns {Object} UTM coordinates + */ + static pointToUTM(point) { + return UTMConverter.toUTM(point.lng, point.lat); + } + + /** + * Convert multiple coordinates to UTM + * @param {Coordinate[]} coordinates - Array of [longitude, latitude] pairs + * @returns {Object[]} Array of UTM coordinates + */ + static coordinatesToUTM(coordinates) { + return coordinates.map((coord) => UTMConverter.toUTM(coord[0], coord[1])); + } +} + +// Convenience function for direct coordinate conversion +/** + * Convert WGS84 to UTM coordinates + * @param {number} lng - Longitude + * @param {number} lat - Latitude + * @returns {Object} UTM coordinates {easting, northing, zone, hemisphere} + */ +export function wgs84ToUTM(lng, lat) { + return UTMConverter.toUTM(lng, lat); +} + +/** + * Convert UTM to WGS84 coordinates + * @param {number} easting - UTM easting + * @param {number} northing - UTM northing + * @param {number} zone - UTM zone + * @param {string} hemisphere - 'N' or 'S' + * @returns {Object} WGS84 coordinates {longitude, latitude} + */ +export function utmToWGS84(easting, northing, zone, hemisphere = "N") { + return UTMConverter.fromUTM(easting, northing, zone, hemisphere); +} |