Module:Georèferencement
La documentation pour ce module peut être créée à Module:Georèferencement/doc
local Tools = require('Module:Outils')
local Coordinates = require('Module:Coordinates')
local wd = require("Module:Wikidata")
local p = {}
local function dms2dec(value) -- convertit les coordonnées DMS en coordonnées décimales
if tonumber(value) then
return value
end
local dms_object = Coordinates._parsedmsstring(value)
return Coordinates._dms2dec(dms_object)
end
function p.setCoord(args) -- retourne les coordonnées dms ou décimales en degrés ; valeur wikidata par défaut
if args.latitude and args.longitude then
local latitude = dms2dec(args.latitude)
local longitude = dms2dec(args.longitude)
return latitude, longitude
end
local claim = wd.getClaims({property = 'P625'})
if claim then
local coords = wd.formatSnak(claim[1].mainsnak)
return coords.latitude, coords.longitude
else
return 0, 0
end
end
function p.setRadCoord(args) -- retourne les coordonnées dms ou décimales en radian ; valeur wikidata par défaut
local latitude, longitude = p.setCoord(args)
return math.rad(latitude), math.rad(longitude)
end
function p.footer(date) -- ajoute la date de consultation
if date ~= nil and date ~= ' 'and date ~= '' then
return ' <small>(viu lo ' .. date .. ')</small>'
end
return ''
end
function p.noCoord(args) -- test l'existance de coordonnées entrées ou de coordonnées wikidata
local latitude, longitude = p.setCoord(args)
if latitude == 0 and longitude == 0 then
return 'Coordonâs manquentes.'
end
end
local function ellipsoide(nom, unit)
--[[
valeurs des paramètres de l'ellipsoide
*** Paramètres ***
- nom : nom de l'ellipsoide
- unit : unité ('km')
*** Retour ***
liste params = {a, b, f, e} avec a : demi grand axe ; b : demi petit axe : f : aplatissement, e : première excentricité
]]
if nom == 'WGS84' or nom == 'WGS 84' then
params = {6378137.0, 6356752.314245179497563967, 1/298.257223563}
elseif nom == 'GRS80' or nom == 'GRS 80' then
params = {6378137.0, 6356752.314140355847852106, 1/298.257222101}
elseif nom == 'Airy1830' or nom =='Airy 1830' then
params = {6377563.396, 6356256.909, 1/299.3249646}
elseif nom == 'Bessel' then
params = {6377397.155, 6356078.963, 1/299.1528153513233}
end
if unit == 'km' then
params[1] = params[1]/1000
params[2] = params[2]/1000
end
params[4] = math.sqrt((params[1]^2-params[2]^2)/params[1]^2)
return params
end
local function geodesicDatum(datum, unit)
--[[
valeurs des paramètres du système géodésique
*** Paramètres ***
- datum : système géodésique
- unit : unité ('km')
Pour plus d'informations, voir : https://fr.wikipedia.org/wiki/Syst%C3%A8me_g%C3%A9od%C3%A9sique
*** Retour ***
liste params = {a, b, f, e} avec a : demi grand axe, b : demi petit axe, f : aplatissement, e : première excentricité
liste helmert : coefficients de helmert
]]
if datum == 'WGS84' or datum == 'WGS 84' then
params = ellipsoide('WGS84', unit)
helmert = {0, 0, 0, 0, 0, 0, 0}
elseif datum == 'ETRS89' or datum == 'ETRS 89' then
params = ellipsoide('GRS80', unit)
helmert = {0, 0, 0, 0, 0, 0, 0}
elseif datum == 'OSGB36' or datum =='OSGB 36' then
params = ellipsoide('Airy1830', unit)
helmert = {-446.448, 125.157, -542.060, 20.4894, -0.1502, -0.2470, -0.8421}
elseif datum == 'Bessel' then
params = ellipsoide('Bessel', unit)
helmert = {-674.374, -15.056, -405.346, 0, 0, 0, 0}
end
return params, helmert
end
function p.toCartesian(lat, lon, height, datum)
--[[
convertit lat/lon en coordonnées cartésiennes
*** Paramètres ***
- lat, lon : latitude et longitude en radians
- height : hauteur (0 par défaut)
- datum : stsyème géodésique
*** Retour ***
x, y, z : coordonnées cartésiennes
]]
local datum = datum or 'WGS84'
local datumParams, _ = geodesicDatum(datum)
local eSq = datumParams[4]*datumParams[4] -- première excentricité au carré ≡ (a²−b²)/a²
local h = height or 0
local v = datumParams[1]/math.sqrt(1-eSq*math.sin(lat)*math.sin(lat))
local x = (v+h) * math.cos(lat) * math.cos(lon)
local y = (v+h) * math.cos(lat) * math.sin(lon)
local z = (v*(1-eSq)+h) * math.sin(lat)
return x, y, z
end
function p.helmertTransform(x, y, z, toDatum)
--[[
convertit les coordonnées cartésiennes dans un autre le datum avec la transformation de helmert
*** Paramètres ***
- x, y, z : coordonnées cartésiennes
- toDatum : système géodésique visé
*** Retour ***
- x1, y1, z1 : coordonnées cartésiennes
Pour plus d'informations, voir : https://fr.wikipedia.org/wiki/Transformation_de_Helmert
]]
local _, helmert = geodesicDatum(toDatum)
local tx = helmert[1] -- déplacement x en metres
local ty = helmert[2] -- déplacement y en metres
local tz = helmert[3] -- déplacement z en metres
local s = helmert[4]/1e6 + 1 -- echelle : ppm en (s+1)
local rx = math.rad(helmert[5]/3600) -- rotation en x : arcseconds en radians
local ry = math.rad(helmert[6]/3600) -- rotation en y : arcseconds en radians
local rz = math.rad(helmert[7]/3600) -- rotation en z : arcseconds en radians
local x1 = tx+x*s-y*rz+z*ry;
local y1 = ty+x*rz+y*s-z*rx;
local z1 = tz-x*ry+y*rx+z*s;
return x1, y1, z1
end
function p.cartesianToLatLon(x, y, z, datum)
--[[
convertit les coordonnées cartésiennes en lat/lon
*** Paramètres ***
- x, y, z : coordonnées cartésiennes
- datum : système géodésique
*** Retour ***
- lat, lon : latitude et longitude
- h : hauteur
]]
local datum = datum or 'WGS84'
local datumParams, _ = geodesicDatum(datum)
local a = datumParams[1]
local b = datumParams[2]
local f = datumParams[3]
local e2 = 2*f-f*f -- première excentricité au carré ≡ (a²−b²)/a²
local eps2 = e2/(1-e2) -- seconde excentricité au carré ≡ (a²−b²)/b²
local p = math.sqrt(x*x+y*y) -- distance de l'axe mineur
local R = math.sqrt(p*p+z*z) -- rayon polaire
-- latitude paramétrique
local tanBeta = (b*z)/(a*p)*(1+eps2*b/R)
local sinBeta = tanBeta/math.sqrt(1+tanBeta*tanBeta)
local cosBeta = sinBeta/tanBeta
-- latitude géodésique
local lat = math.atan2(z+eps2*b*sinBeta*sinBeta*sinBeta, p-e2*a*cosBeta*cosBeta*cosBeta)
-- longitude
local lon = math.atan2(y, x)
-- hauteur
local sinLat = math.sin(lat)
local cosLat = math.cos(lat)
local v = a/math.sqrt(1-e2*sinLat*sinLat)
local h = p*cosLat+z*sinLat-(a*a/v)
return lat, lon, h
end
function p.convertDatum(lat, lon, height, datum, toDatum)
--[[
Modifie les coordonnées lat/lon d'un système géodésique à un autre
*** Paramètres ***
- lat, lon : latitude et longitude en radians
- height : hauteur
- datum : système géodésique initial
- toDatum : système géodésique visé
*** Retour ***
- lat, lon : latitude et longitude dans le système géodésique visé en radians
- height : hauteur dans le système géodésique visé
]]
local x, y, z = p.toCartesian(lat, lon, height, datum)
local x, y, z = p.helmertTransform(x, y, z, toDatum)
local lat, lon, h = p.cartesianToLatLon(x, y, z, toDatum)
return lat, lon, h
end
function p.pseudoMercator(lat, lon, datum)
--[[
Projection Web Mercator (ou Pseudo Mercator)
*** Paramètres ***
- lat, lon : latitude et longitude en radians
- datum : système géodésique utilisé. Par défaut, WGS84
Pour plus d'informations, voir : https://en.wikipedia.org/wiki/Web_Mercator_projection
]]
local datum = datum or 'WGS84'
local datumParams, _ = geodesicDatum(datum)
local X = lon*datumParams[1]
local Y = math.log(math.tan(math.pi/4+lat/2))*datumParams[1] -- équivalent à a*0.5*math.log((1+math.sin(lat))/(1-math.sin(lat)))
return X, Y
end
function p.lambert(lat, lon, lat_0, lat_1, lat_2, lon_0, x_0, y_0, datum)
--[[
Projection conique conforme de Lambert
*** Paramètres ***
- lat, lon : latitude et longitude en radians
- lat_0, lon_0 : latitude et longitude du centre de la projection, en degrés
- lat_1 : first standard parallel
- lat_2 : second standard parallel
- x_0 : false easting
- y_0 : false northing
- datum : système géodésique utilisé. Par défaut, WGS84
Pour plus d'informations, voir : https://fr.wikipedia.org/wiki/Projection_conique_conforme_de_Lambert et https://proj.org/operations/projections/lcc.html?highlight=lcc
]]
local datum = datum or 'WGS84'
local datumParams, _ = geodesicDatum(datum)
local a = datumParams[1]
local e = datumParams[4]
local lat_0 = math.rad(lat_0)
local lat_1 = math.rad(lat_1)
local lat_2 = math.rad(lat_2)
local lon_0 = math.rad(lon_0)
local n_num = math.log(math.cos(lat_2)/math.cos(lat_1)) + 0.5*math.log((1-e^2*math.sin(lat_1)^2)/(1-e^2*math.sin(lat_2)^2))
local n_dem = math.log( math.tan(math.pi/4+lat_1/2)*(1-e*math.sin(lat_1))^(e/2)*(1+e*math.sin(lat_2))^(e/2) / (math.tan(math.pi/4+lat_2/2)*(1+e*math.sin(lat_1))^(e/2)*(1-e*math.sin(lat_2))^(e/2)) )
local n = n_num/n_dem
local init = a*math.cos(lat_1)/(n*math.sqrt(1-e^2*math.sin(lat_1)^2)) * ((math.tan(math.pi/4+lat_1/2))*((1-e*math.sin(lat_1))/(1+e*math.sin(lat_1)))^(e/2))^n
local rho = init * ((1/math.tan(math.pi/4+lat/2))*((1+e*math.sin(lat))/(1-e*math.sin(lat)))^(e/2))^n
local rho_0 = init * ((1/math.tan(math.pi/4+lat_0/2))*((1+e*math.sin(lat_0))/(1-e*math.sin(lat_0)))^(e/2))^n
local X = math.floor(x_0 + rho * math.sin(n * (lon - lon_0)))
local Y = math.floor(y_0 + rho_0 - rho * math.cos(n * (lon - lon_0)))
return X, Y
end
function p.utm(lat, lon, fuseau, k_0, datum, lat_0, E_0, N_0)
--[[
Projection transverse universelle de Mercator
*** Paramètres ***
- lat, lon : latitude et longitude en radians
- fuseau : fuseau UTM utilisé pour lon_0 (longitude de référence)
- k_0 : facteur d'échelle. Par défaut 0.9996
- datum : système géodésique utilisé. Par défaut, WGS84
- lat_0 : latitude de référence, en degrés. Par défaut, 0
- E_0 : false easting
- N_0 : false northing
Pour plus d'informations, voir : https://fr.wikipedia.org/wiki/Transverse_universelle_de_Mercator et https://proj.org/operations/projections/tmerc.html?highlight=tmerc
]]
local datum = datum or 'WGS84'
local datumParams, _ = geodesicDatum(datum, 'km')
local a = datumParams[1]
local e = datumParams[4]
local lon_0 = math.rad(-183+6*fuseau)
local E_0 = E_0 or 500
local N_0 = N_0 or 0
local v = 1/math.sqrt(1-e^2*math.sin(lat)^2)
local A = (lon-lon_0)*math.cos(lat)
local s1 = 1-e^2/4-3*e^4/64-5*e^6/256
local s2 = 3*e^2/8+3*e^4/32+45*e^6/1024
local s3 = 15*e^4/256+45*e^6/1024
local s4 = 35*e^6/3072
local s = s1*lat-s2*math.sin(2*lat)+s3*math.sin(4*lat)-s4*math.sin(6*lat)
local T = math.tan(lat)^2
local C = e^2/(1-e^2)
local k_0 = k_0 or 0.9996
local E = E_0+k_0*a*v*(A+(1-T+C)*A^3/6+(5-18*T+T^2)*A^5/120)
local E = 1000*E -- en mètres
if lat_0 then
local s0 = s1*math.rad(lat_0)-s2*math.sin(2*math.rad(lat_0))+s3*math.sin(4*math.rad(lat_0))-s4*math.sin(6*math.rad(lat_0))
N = N_0+k_0*a*(s-s0+v*math.tan(lat)*(A^2/2+(5-T+9*C+4*C^2)*A^4/24+(61-58*T+T^2)*A^6/720))
else
N = N_0+k_0*a*(s+v*math.tan(lat)*(A^2/2+(5-T+9*C+4*C^2)*A^4/24+(61-58*T+T^2)*A^6/720))
end
local N = 1000*N -- en mètres
return E, N
end
return p