User:Maxbe/Kartenversuch/paesseausrichten.php

<?php

// srtmele() liefert zu einem Punkt (x,y in EPSG 900913) die passende Hoehe. Dazu sucht man in der Datenbank $dbs die naechstgelegenen Hoehenwerte
// und holt macht eine bilineare Interpolation aus den 4 angrenzenden Werten. Falls es am Kartenrand keine 4 Werte gibt, nimmt man den einen naechsten
// Dass man erst 6 Werte holt und dann die 4 angrenzenden daraus aussucht, liegt daran, dass wir zwar immer von "90m-Raster" reden, aber in echt
// ein "3-Sekunden-Raster" haben. Die Hoehenwerte liegen nicht im Quadrat und bei unguenstiger Lage des gsuchten Wertes koennten die 4 naechsten
// Hoehenwerte auch T-foermig um diesen Punkt verteilt sein.

function srtmele($x,$y,$dbs){
 $srtmquery="SELECT ele,distance(ST_GeometryFromText('POINT($x $y)',900913),way) as entfernung, st_x(way) as x, st_y(way) as y FROM elepoint WHERE way && expand(ST_GeometryFromText('POINT($x $y)',900913),300) order by entfernung limit 6;";
 $resultsrtm=pg_query($dbs,$srtmquery);
 $h1=$h2=$h3=$h4=$h5=$hz=0;
 while($hoehe=pg_fetch_array($resultsrtm)){
  if(!($h1&&$h2&&$h3&&$h4)){$h5=$hoehe['ele'];}
  if((!$h1)&&($hoehe['x']<=$x)&&($hoehe['y']<=$y)){$h1=$hoehe['ele'];$d1=$hoehe['entfernung'];$x1=$hoehe['x'];$y1=$hoehe['y'];$hz++;}
  if((!$h2)&&($hoehe['x']> $x)&&($hoehe['y']<=$y)){$h2=$hoehe['ele'];$d2=$hoehe['entfernung'];$x2=$hoehe['x'];$y2=$hoehe['y'];$hz++;}
  if((!$h3)&&($hoehe['x']<=$x)&&($hoehe['y']> $y)){$h3=$hoehe['ele'];$d3=$hoehe['entfernung'];$x3=$hoehe['x'];$y3=$hoehe['y'];$hz++;}
  if((!$h4)&&($hoehe['x']> $x)&&($hoehe['y']> $y)){$h4=$hoehe['ele'];$d4=$hoehe['entfernung'];$x4=$hoehe['x'];$y4=$hoehe['y'];$hz++;}
 }
 if($hz==4){
  $h=$h1/(($x2-$x1)*($y3-$y1))*($x2-$x)*($y3-$y)+
     $h2/(($x2-$x1)*($y3-$y1))*($x-$x1)*($y3-$y)+
     $h3/(($x2-$x1)*($y3-$y1))*($x2-$x)*($y-$y1)+
     $h4/(($x2-$x1)*($y3-$y1))*($x-$x1)*($y-$y1);
 }else{$h=$h5;}
 return $h;
}


// --------------- Hier gehts los -------------------------------------------------

// dbs: Datenbank mit einer Tabelle elepoint, die alle paar Meter die SRTM-Hoehen enthaelt
// db:  OSM-Datenbank mit osm_poi, die die interessanten POIs enthaelt

$dbs=pg_pconnect('host=localhost port=5432 dbname=srtm user=osm password=osm') or die("Keine DB, sonst alles in Ordnung");
$db=pg_pconnect('host=localhost port=5432 dbname=osm user=osm password=osm') or die("Keine DB, sonst alles in Ordnung");

// Wir holen alle Saettel, Paesse und Cols und machen eine Schleife

$saddlequery="select *,st_x(way) as x,st_y(way) as y from osm_poi where (nature='saddle' or nature='mountain_pass' or nature='col');";
$resultsaddle=pg_query($db,$saddlequery);
while ($s=pg_fetch_array($resultsaddle)){
 $osmid=$s['osm_id'];

// holen die Hoehe des Sattels aus der DB
// laufen in 15-Grad-Schritten in 80 ("Meter" in EPSG:900913) Abstand um den Punkt
// und bestimmen die Differenz aus dieser Hoehe und der Sattelhoehe

 $h=srtmele($s['x'],$s['y'],$dbs);
 for($w=0;$w<24;$w++){
  $co=cos($w*15*2*3.1415927/360);
  $si=sin($w*15*2*3.1415927/360);
  $x1=$s['x']+$co*80;
  $y1=$s['y']+$si*80;
  $ho[$w]=srtmele($x1,$y1,$dbs)-$h;
 }
 $mh=100000;$wi=45;

// Wir laufen nochmal einen Halbkreis um den Punkt und bestimmen
// Das Minimum aus den Hoehendifferenzen (vorne-Sattel)+(hinten-Sattel)-(rechts-Sattel)-(links-Sattel)
// Die Nachbarwerte (vorne+-1)... werden mit 1/3 Gewichtung einbezogen, die Nachbarn (vorne+-2) mit 1/7
//
// Gewonnen hat die Richtung mit dem kleinsten Wert.

 for($w=0;$w<12;$w++){
    $dh=($ho[$w]       +$ho[($w+12)%24]
         +0.33*($ho[($w+23)%24]+$ho[($w+1)%24] +$ho[($w+11)%24] +$ho[($w+13)%24])+
         +0.14*($ho[($w+22)%24]+$ho[($w+2)%24] +$ho[($w+10)%24] +$ho[($w+14)%24]))-
        ($ho[($w+6)%24]+$ho[($w+18)%24]
         +0.33*($ho[($w+5)%24] +$ho[($w+7)%24] +$ho[($w+17)%24] +$ho[($w+19)%24])+
         +0.14*($ho[($w+4)%24] +$ho[($w+8)%24] +$ho[($w+16)%24] +$ho[($w+20)%24]));
  if($dh<$mh){$mh=$dh;$wi=$w*15;}
 }
 
// Die Richtung stecken wir in das Feld "population", weil das ist da und wird bei Saetteln sonst nicht gebraucht ;)
 
 $saddleupdquery="update osm_poi set population=$wi where osm_id=$osmid and (nature='saddle' or nature='mountain_pass' or nature='col') and (population is null or population!='$wi');";
 echo "$saddleupdquery\n";
}
 
?>