Téléchargement et conversion en GPX de la BDCavité (liste des cavités souterraines du BRGM)

cavites.php

<?php

// Mode d'emploi : créer un répertoire vierge et lancer "php cavites.php"

const JSON_URL = 'https://www.georisques.gouv.fr/webappReport/ws/telechargement/cavites?anneemin=2003';

const PROJECTIONS = [
	1 => 'LambertI',
	2 => 'LambertII',
	3 => 'LambertIII',
	4 => 'LambertIV',
	5 => 'LambertIIExtend',
	6 => 'WGS84',
	7 => 'Lambert93',
];

if (!file_exists('csv')) {
	$json = file_get_contents(JSON_URL);

	$json = json_decode($json);

	@mkdir('csv');

	foreach ($json->departemental as $key => $data) {
		$target = sprintf('csv/%s.csv', $key);

		if (!file_exists($target)) {
			copy($data->lien, $target);
			echo ".";
		}
	}

	echo "\n";
	@mkdir("gpx");
}

foreach (glob('csv/*.csv') as $file) {
	echo basename($file) . "...";
	csv_to_gpx($file, sprintf('gpx/%s.gpx', str_replace('.csv', '', basename($file))));
	echo "\n";
}


echo "\n";

function csv_to_gpx(string $src, string $target) {
	$target = fopen($target, 'w');
	$src = fopen($src, 'r');

	fwrite($target, '<?xml version="1.0" encoding="UTF-8"?>
		<gpx
			xmlns="http://www.topografix.com/GPX/1/1"
			xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
			xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd"
			version="1.1">' . "\n");

	$i = 0;
	$header = null;

	while (!feof($src)) {
		$line = fgetcsv($src, null, ';');

		if (!$line) {
			continue;
		}

		// Ignore headers
		if ($i++ == 0) {
			$line = array_filter($line, fn($a) => trim($a) === '' ? false : true);
			$header = $line;
			continue;
		}

		$line = array_slice($line, 0, count($header), true);

		if (count($line) != count($header)) {
			printf("Line %d: field count mismatch\n%s\n", $i, print_r($line, true) . print_r($header, true));
			continue;
		}

		$line = array_combine($header, $line);

		if (!isset($line['confidentialite'], $line['nomCavite'], $line['zOuvrage'], $line['xOuvrage'], $line['yOuvrage'], $line['lambertOuvrage'])) {
			printf("Line %d: Missing mandatory field\n%s\n", $i, print_r($line, true));
			continue;
		}

		if ($line['confidentialite'] != 'public') {
			continue;
		}

		unset($line['confidentialite'], $line['raisonConfidentialite'], $line['organismeChxConfidentialite']);

		$desc = [];

		foreach ($line as $key => $value) {
			$desc[] = sprintf('%s: %s', $key, trim($value));
		}

		$desc = implode("\n", $desc);

		$name = $line['nomCavite'];

		$latlon = convert_coordinates($line);

		if ($latlon === null) {
			printf("Line %d: Invalid ref: %s for %s\n", $i, $line['lambertOuvrage'], $line['numCavite']);
			continue;
		}

		list($lat, $lon) = $latlon;

		fwrite($target, sprintf('<wpt lat="%f" lon="%f"><sym>Tunnel</sym><name>%s</name><desc>%s</desc><ele>%d</ele></wpt>' . "\n",
			$lat, $lon, htmlspecialchars($name, ENT_XML1), htmlspecialchars($desc, ENT_XML1), $line['zOuvrage']));
	}

	fwrite($target, "</gpx>\n");
	fclose($target);
}

function convert_coordinates(array $line) {
	$ref = PROJECTIONS[$line['lambertOuvrage']] ?? null;

	if (!$ref) {
		return null;
	}

	if ($ref == 'LambertIIExtend') {
		$c = new Convert($line['xOuvrage'], $line['yOuvrage']);
		return $c->convert();
	}
	elseif ($ref == 'Lambert93') {
		return lambert93ToWgs84($line['xOuvrage'], $line['yOuvrage']);
	}
	elseif ($ref == 'WGS84') {
		return [$line['yOuvrage'], $line['xOuvrage']];
	}
	else {
		return lambert2gps($line['xOuvrage'], $line['yOuvrage'], $ref);
	}
}

/**
 * https://codes-sources.commentcamarche.net/source/52227-convertisseur-lambert2-etendu-en-coordonnee-geographique-longitude-latitude
 * @author Florent Cardot
 * Y = Latitude
 * X = Longitude
 */
/*
$lon = 5802906.829;
$lat = 6453674.479;

$L = new Convert($lon, $lat);
$L->convert();
*/
class Convert
{
	private $X;
	private $Y;
	private $Coord;
	private $Cm;
	private $n;
	private $XSm;
	private $YSm;
	private $a;
	private $f1;

	/**Contructeur (initialise les variables qui doivent l'etre)**/
	function __construct($X, $Y)
	{
		
		$this->Cm =	11745793.393435;
		$this->n =	0.728968627421412;
		$this->XSm = 600000;
		$this->YSm = 8199695.76800186;
		$this->a = 6378249.2000;
		$this->f1 = 6356515.0000;
		$this->X =	$X - $this->XSm;
		$this->Y =	$Y - $this->YSm;
	}//end function
	
	public function convert()
	{
		return [self::convertY(),self::convertX()];
		$this->Coord[0] = self::ConvertX();  //X
		$this->Coord[1] = self::ConvertY();  //Y
		return $this->Coord;
	}//end function
	
	public function ConvertX()
	{
		$longitude = atan(-($this->X)/($this->Y));
		$longitude = $longitude / $this->n;
		$longitude = $longitude * 180 / pi();
		$constante = 2 + (20 / 60) + (14.025 / 3600);
		$longitude = $longitude + $constante;
		return($longitude);
	}//end function
	
	public function ConvertY()
	{
		$latitude = sqrt(pow($this->X, 2) + pow($this->Y, 2));
		$f  = ($this->a - $this->f1) / $this->a;
		$e² = 2 * $f - pow($f, 2);
		$e  = sqrt($e²);
		$Latiso = log($this->Cm / $latitude) / $this->n;
		$latitude = tanh($Latiso);
		for ($i = 0; $i < 6; $i++)
		{
			$latitude = tanh($Latiso + $e * self::atanh($e * $latitude));
		}
		$latitude = asin($latitude);
		$latitude = $latitude / pi();
		$latitude = $latitude * 180;
		return($latitude);
		//43,87602354
	}//end function
	
	public function atanh($x)
	{
		$resultat = log((1 + $x) / (1 - $x)) / 2;
		return $resultat;
	}//end function

}//end classs


// https://gist.github.com/will83/5920606
function lambert93ToWgs84($x, $y)
{
	$b8  = 1 / 298.257222101;
	$b10 = sqrt(2 * $b8 - $b8 * $b8);
	$b16 = 0.7256077650532670;
	$x = number_format(floatval($x), 10, '.', '') - 700000;
	$y = number_format(floatval($y), 10, '.', '') - 12655612.0499;
	$gamma = atan(-$x / $y);
	$latiso = log(11754255.426096 / sqrt(($x * $x) + ($y * $y))) / $b16;
	$sinphiit = tanh($latiso + $b10 * atanh($b10 * sin(1)));

	for ($i = 0; $i != 6 ; $i++) {
		$sinphiit = tanh($latiso + $b10 * atanh($b10 * $sinphiit));
	}
	return [
		asin($sinphiit) / pi() * 180, // latitude
		($gamma / $b16 + 3 / 180 * pi()) / pi() * 180, // longitude
	];
}

// https://gist.github.com/lovasoa/096d8be82520ea6b0abe
function lambert2gps($x, $y, $lambert) {
      $lamberts = array(
       "LambertI" => 0,
       "LambertII" => 1,
       "LambertIII" => 2,
       "LambertIV" => 3,
       "LambertIIExtend" => 4,
       "Lambert93" => 5
      );
      $index = $lamberts[$lambert];
      $ntabs =  array(0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322, 0.7289686274, 0.7256077650);
      $ctabs =  array(11603796.98, 11745793.39, 11947992.52, 12136281.99, 11745793.39, 11754255.426);
      $Xstabs = array(600000.0, 600000.0, 600000.0, 234.358, 600000.0, 700000.0);
      $Ystabs = array(5657616.674, 6199695.768, 6791905.085, 7239161.542, 8199695.768, 12655612.050);

      $n  = $ntabs [$index];
      $c  = $ctabs [$index];            // En mètres
      $Xs = $Xstabs[$index];          // En mètres
      $Ys = $Ystabs[$index];          // En mètres
      $l0 = 0.0;                    //correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich
      $e = 0.08248325676;           //e du NTF (on le change après pour passer en WGS)
      $eps = 0.00001;     // précision


      /***********************************************************
      *  coordonnées dans la projection de Lambert 2 à convertir *
      ************************************************************/
      $X = $x;
      $Y = substr($y, 1);

      /*
       * Conversion Lambert 2 -> NTF géographique : ALG0004
       */
        $R = Sqrt((($X - $Xs) * ($X - $Xs)) + (($Y - $Ys) * ($Y - $Ys)));
        $g = Atan(($X - $Xs) / ($Ys - $Y));

        $l = $l0 + ($g / $n);
        $L = -(1 / $n) * Log(Abs($R / $c));


        $phi0 = 2 * Atan(Exp($L)) - (pi() / 2.0);
        $phiprec = $phi0;
        $phii = 2 * Atan((Pow(((1 + $e * Sin($phiprec)) / (1 - $e * Sin($phiprec))), $e / 2.0) * Exp($L))) - (pi() / 2.0);

        while (!(Abs($phii - $phiprec) < $eps)) {
                $phiprec = $phii;
                $phii = 2 * Atan((Pow(((1 + $e * Sin($phiprec)) / (1 - $e * Sin($phiprec))), $e / 2.0) * Exp($L))) - (pi() / 2.0);
        }

        $phi = $phii;

    /*
     * Conversion NTF géogra$phique -> NTF cartésien : ALG0009
     */
    $a = 6378249.2;
    $h = 100;         // En mètres

    $N = $a / (Pow((1 - ($e * $e) * (Sin($phi) * Sin($phi))), 0.5));
    $X_cart = ($N + $h) * Cos($phi) * Cos($l);
    $Y_cart = ($N + $h) * Cos($phi) * Sin($l);
    $Z_cart = (($N * (1 - ($e * $e))) + $h) * Sin($phi);

    /*
     * Conversion NTF cartésien -> WGS84 cartésien : ALG0013
     */

    // Il s'agit d'une simple translation
    $XWGS84 = $X_cart - 168;
    $YWGS84 = $Y_cart - 60;
    $ZWGS84 = $Z_cart + 320;


    /*
     * Conversion WGS84 cartésien -> WGS84 géogra$phique : ALG0012
     */
    
    $l840 = 0.04079234433;    // 0.04079234433 pour passer dans un référentiel par rapport au méridien
    // de Greenwich, sinon mettre 0
    
    $e = 0.08181919106;              // On change $e pour le mettre dans le système WGS84 au lieu de NTF
    $a = 6378137.0;

    $P = Sqrt(($XWGS84 * $XWGS84) + ($YWGS84 * $YWGS84));

    $l84 = $l840 + Atan($YWGS84 / $XWGS84);

    $phi840 = Atan($ZWGS84 / ($P * (1 - (($a * $e * $e))
                                / Sqrt(($XWGS84 * $XWGS84) + ($YWGS84 * $YWGS84) + ($ZWGS84 * $ZWGS84)))));

    $phi84prec = $phi840;

    $phi84i = Atan(($ZWGS84 / $P) / (1 - (($a * $e * $e * Cos($phi84prec))
            / ($P * Sqrt(1 - $e * $e * (Sin($phi84prec) * Sin($phi84prec)))))));

    while (!(Abs($phi84i - $phi84prec) < $eps))
    {
        $phi84prec = $phi84i;
        $phi84i = Atan(($ZWGS84 / $P) / (1 - (($a * $e * $e * Cos($phi84prec))
                / ($P * Sqrt(1 - (($e * $e) * (Sin($phi84prec) * Sin($phi84prec))))))));

    }

    $phi84 = $phi84i;

    return array($phi84 * 180.0 / pi(), $l84 * 180.0 / pi());
}