Blog der Heimetli Software AG

Umriss der Schweiz mit eingezeichneten Städten

D3.js kann auch Karten darstellen wenn man ein passendes GeoJSON-File hat.

© OpenStreetMap contributors

Daten beschaffen und konvertieren

Meine Suche brachte kein solches File zu Tage, aber auf planet.osm.ch fand ich ein Textfile mit den Koordinaten der Schweizer Grenze. Es handelt sich dabei um einen Auszug aus den Daten von openstreetmap.org.

Das Fileformat ist TSV (Tab Separated Values), also gut zu verarbeiten mit Linux-Tools.

Als erstes habe ich mir die Outline eines GeoJSON-Files aus der Wikipedia kopiert und aufs notwendigste reduziert:

{ "type": "FeatureCollection",
    "features": [
      { "type": "Feature",
         "geometry": {
           "type": "Polygon",
           "coordinates": [ [
#include switzerland-exact.poly
           ] ]
         },
         "properties": {}
      }
    ]
}

Das #include markiert den Ort wo die Daten aus dem .poly eingefügt werden sollen.

Ein awk-Script machte daraus ein GeoJSON:

awk '{
        if( $1 == "#include" )
        {
           while( (getline line < $2) > 0 )
           {
              if( split(line,coords,"\t") == 3 )
              {
                 print "[" coords[2] "," coords[3] "],";
              }
           }
           close($2) ;
        }
        else
        {
           print ;
        }
     }' geo.template

Bis darauf dass hinter dem letzten Array-Element kein Komma stehen sollte, ist das GeoJSON valid. Das habe ich mit dem Editor weggeputzt, weil es mir den Aufwand für ein besseres Script nicht wert war.

Das JSON funktioniert problemlos, aber es ist gewaltig! Es hat über 27000 Zeilen und ist mehr als 790kB gross.

Daten reduzieren

JSON ist zwar nicht gerade ein effizientes Format, aber daran liegt es in diesem Fall nicht. Die Daten von openstreetmap sind einfach zu detailliert.

Also ab in Internet und nachforschen was man da tun könnte. Es gibt Online-Konverter die versprechen GeoJSON-Files zu vereinfachen, aber einer davon hat mir mal ein File vermurkst. Selbstverständlich gibt es jede Menge Applikationen die ebenfalls solche Funktionen anbieten, aber ich kenne keine davon.

Bald einmal stiess ich auf den Ramer-Douglas-Peucker Algoritmus der anscheinend sehr beliebt ist. Die Wikipedia-Beschreibung ist ausführlich genug um ihn zu implementieren.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

typedef struct
{
   double x ;
   double y ;
   bool   valid ;
} POINT ;

POINT points[30000] ;

/*************************************************/
 double distance( POINT &p1, POINT &p2, POINT &p )
/*************************************************/

{
   return abs((p2.y-p1.y) * p.x - (p2.x - p1.x) * p.y + p2.x * p1.y - p1.x * p2.y) /
	  sqrt( (p2.y - p1.y) * (p2.y - p1.y) + (p2.x - p1.x) * (p2.x - p1.x) ) ;
}

/****************************************************/
 void simplify( int first, int last, double epsilon )
/****************************************************/

{
   int    index = 0 ;
   double max   = 0 ;

   for( int i = first + 1; i < last; i++ )
   {
      double d = distance( points[first], points[last], points[i] ) ;

      if( d > max )
      {
         max   = d ;
         index = i ;
      }
   }

   if( max > epsilon )
   {
      simplify( first, index, epsilon ) ;
      simplify( index, last,  epsilon ) ;
   }
   else
   {

      points[first].valid = true ;
      points[last].valid  = true ;
   }
}


/**********************************/
 int main( int argc, char *argv[] )
/**********************************/

{
   if( argc != 2 )
   {
      printf( "usage: rdp epsilon\n" ) ;
      return 1 ;
   }

   double epsilon = atof( argv[1] ) ;
 
   int    index = 0 ;

   while( scanf("%lf%lf",&points[index].x,&points[index].y) == 2 )
      points[index++].valid = false ;

   points[0].valid       = true ;
   points[index-1].valid = true ;

   simplify( 0, index - 1, epsilon ) ;

   printf( "{ \"type\": \"FeatureCollection\",\n" ) ;
   printf( "  \"features\": [\n"                  ) ;
   printf( "   { \"type\": \"Feature\",\n"        ) ;
   printf( "       \"geometry\": {\n"             ) ;
   printf( "         \"type\": \"Polygon\",\n"    ) ;
   printf( "        \"coordinates\": [ [\n"       ) ;

   for( int i = 0; i < index - 1; i++ )
   {
      if( points[i].valid )
	 printf( "[%f,%f],\n", points[i].x, points[i].y ) ;
   }

   printf( "[%f,%f]\n", points[0].x, points[0].y  ) ;

   printf( "          ] ]\n"                      ) ;
   printf( "      },\n"                           ) ;
   printf( "      \"properties\": {}\n"           ) ;
   printf( "   }\n"                               ) ;
   printf( " ]\n"                                 ) ;
   printf( "}\n"                                  ) ;

   return 0 ;
}

Das Programm erwartet Epsilon als Argument und das .poly als Standard-Input. Mit diesem Programm konnte ich die Filegrösse auf 14kB reduzieren. Wahrscheinlich enthält es immer noch zu viele Details, aber irgendwann wollte ich auch mal die Karte sehen ;-)

Aufgepasst: der erste und der letzte Punkt haben genau die gleichen Koordinaten. Das machte dem Algorithmus das Leben schwer, bis ich die letzte Stelle der Endkoordinaten um eins erhöhte.

Das Geo-JSON mit dem Umriss der Schweiz können Sie hier herunterladen. Aber nicht vergessen dass das Copyright an den Daten bei OpenStreetMap liegt.