DE:SpatiaLite/Beispiele/Messtischblatt-TK25-Nummern-Bundesländern-zuordnen

From OpenStreetMap Wiki
Jump to navigation Jump to search
SpatiaLite/Examples/Mapping SurveySheets to FederalState names - Other languages

Einführung

Dieses Skript erzeugt eine SpatiaLite Datenbank mit zwei Nutzertabellen:

  • states (Bundesländer) mit den Spalten name (Name des Bundeslandes), source (Indikator für die Geometrieherkunft), g (Geometrie, WGS84, XY Geokoordinaten)
  • mtbmbr (Messtischblätter mit minimal umgebenden Rechtecken, d.h. Viereck des Kartenausschnitts) mit den Spalten name (Name des Messtischblatts, derzeit ungenutzt), nr (Rasternummer des MTB), g (Geometriespalte, Bessel, XY Geokoordinaten)

Es wurde mit dem Wunsch im Hinterkopf erstellt, die Messtischblatt- (fortan MTB) bzw. TK25-Kartennummern Bundesländern zuordnen zu können. Zudem sollen die Geometrien der Blattschnitte mit denen der Bundesländer spatial, d.h. räumlich verschnitten, bzw. verrechnet werden können. Je nach Abfrage sind unten genannte Probleme dann leicht zu lösen.

Die Kartennummern der MTB wurden ursprünglich aus einem zweidimensionalen Raster abgeleitet. Von einer Kartennummer lässt sich direkt das minimal umgebende Rechteck des Karten-/Blattschnitts ableiten. Ein Import von Geometrien externer Quellen nach mtbmbr ist daher nicht nötig.

Anders sieht es für die Geometrien der Bundeslandflächen, bzw. -grenzen, aus. Diese werden zum einen durch einen Import der Polygone approximiert, die auch von der GeoFabrik zum Segmentieren des osm planets verwendet werden, hier finden Überlappungen an den Landesgrenzen statt, und zum anderen mittels eines Abrufs der Landesgrenz-Relationen von Overpass API, die nicht überlappend aber umfangreicher sind.

Momentan ist die durch das Skript erzeugte DB nützlich, um folgende Probleme / Fragestellungen zu lösen:

  • alle MTB-Kartennummern ausgeben, die ein bestimmtes Bundesland abdecken
  • die Bundesländer ausgeben, die vom einem bestimmten MTB abgedeckt werden (falls eine Landesgrenze im Blattschnitt verläuft, die nicht zugleich Staatengrenze ist, sind das 2 oder mehr, sonst 1)
  • alle MTB-Kartennummern ausgeben, deren Inhalte eine Landesgrenze komplett abdecken (grenzmtb Sicht in der DB)
  • genau die MTB-Nummern ermitteln, welche die Grenzregion zweier benachbarter Bundesländer abdeckt
  • den abgebildeten Flächenanteil eines im Blattschnitt dargestellten Bundeslandes, relativ zur insgesamt abgebildeten Fläche des Blattschnittes; dies ist insbes. bei Blattschnitten interessant innerhalb derer eine Landesgrenze verläuft und so mehrere Bundesländer abgebildet sind
  • den Flächenunterschied/differenz zwischen dem für das Ausschneiden der Daten eines Bundeslandes genutzten Polygons und dessen tatsächlich in OSM erfasster Fläche errechnen
  • ...

states wird zweimal 'gefüllt'

  • der erste Lauf lädt die groben Geometrien der Schnittpolygone von Geofabrik's Servern und liest sie in die Tabelle ein (source=poly)
  • der zweite Lauf fragt mittels Overpass API nach den exakteren, umfangreicheren GrenzverlaufsRelationen (source=osm) und benutzt osmjs für das Einlesen der Daten in die Tabelle, evtl. ist auch eine Variante mit readosm möglich

mtbmbr wird mit den MBRs/MURs/BBOXes der Blattschnitte gefüllt, die für Deutschland in den heutigen Grenzen relevant sind

  • die bounding boxes werden dabei schrittweise durch das Skript erzeugt
  • insgesamt hat die Tabelle 89x59 Einträge, viele davon repräsentieren Blattschnitte, die historisch weder vermessen, noch als Kartenblatt publiziert wurden

Das Programm

  • wurde mit Kubuntu 14.04 LTS Live/Install bootbar von USB erstellt und getestet, Hardware-Eckdaten: corei3 CPU, 3G RAM
  • holt benötigte Abhängigkeiten mittels "sudo apt-get", darunter auch spatialite selbst, falls es nicht installiert ist
  • benötigt Schreibzugriff in /tmp' und dem Verzeichniss, in dem das ausführbare Bash-Skript liegt
  • braucht eine jüngere Version von osmjs
    • das deb Paket, das auf den Archivservern des Ubunutu-Projekts bereitgestellt wird, ist zu alt und stürzte im Test mit der Option -l sparsetable regelmäßig ab, da der RAM erschöpft war
    • auf (Live-)Systemen mit knappem Speicher empfehlen sich ausnahmslos Versionen von osmjs, welche die Option -l vector verstehen
    • ist osmjs nicht installiert, baut das Skript selbst eine Version aus den GIT Quellen von Osmium und extrahiert das binary osmjs, um diese Teilschritte bei wiederholten Aufrufen nicht wiederholen zu müssen

Beispiele für Abfragen sind am Ende des Skriptes zu finden. Wer selbst experimentieren und eigene Abfragen entwickeln bzw. testen möchte, kann entweder die CLI-Shell oder spatialite_gui mit der DB verbinden. Letzteres enthält u.a. einen dialogbasierten Abfragegenerator und hilft evtl. Einsteigern.

Ausgaben

Das Skript erzeugt folgende Dateien

  • die SQLite DB (~30mb) mit
    • den Tabellen samt Geometriespalten und räumlichem Index
    • gespeicherten Sichten, die das Lösen o.a. Problemstellungen erleichtern
  • eine Textdatei, welche für jede MTB-Kartennummer das zugehörige Bundesland listet (~70kb unkomprimiert), ein Eintrag pro Zeile und Zeichen | als Feldtrenner; Grundlage für die Wikiliste
    • erzeugt mittels der Sicht mtb2state

Siehe auch

Listing

#!/bin/bash

# urls
NADURL="http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/BETA2007.gsb"
POLYURLS="http://download.geofabrik.de/europe/germany/"
OVERPASS="http://overpass-api.de/api/interpreter"

# permanent files
MTBMAP="${0%.sh}.mtbmap"
NAD="${0%.sh}.beta2007"
DB="${0%.sh}.sqlite"

# temporary files
QF="$(mktemp --dry-run --suffix=.sql)"
OF="${QF//./-}.osm" # else osmjs deducts pbf from filename


which spatialite &> /dev/null \
  || sudo apt-get install spatialite-{bin,gui} proj-bin \
  || exit 1


# ----------------------------------------------------------------
# do not recreate the db, if it exists
# ----------------------------------------------------------------
[[ -f "$DB" ]] || { # START DB creation

cat <<EOF > "$QF"
CREATE TABLE IF NOT EXISTS states(
  id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT,
  name TEXT,
  source TEXT
);
SELECT AddGeometryColumn(
  'states', 'g', 4326, 'MULTIPOLYGON', 'XY'
);
SELECT CreateSpatialIndex(
  'states', 'g'
);
DELETE FROM states;
EOF


wget -qO - "$POLYURLS" |\
  sed -e 's/href=./\n/g' |\
  sed -e 's/".*//' |\
  grep ".poly$" |\
  while read p
  do
    echo "INSERT INTO states VALUES(NULL, \"${p%.poly}\", \"poly\", " >> "$QF"
    echo -n " GeomFromText('MULTIPOLYGON(((" >> "$QF"
    wget -qO - "$POLYURLS/$p" |\
      while read l k
      do
        [[ -z "$k" && "$l" -gt 1 ]] && echo -en ")),(("
        [[ -n "$k" ]] && LANG=C printf "%f %f," "$l" "$k"
      done |\
      sed -e 's/$/))/' -e 's/,)/)/g' >> "$QF"
    echo -e ")', 4326)\n);" >> "$QF"
  done


OSMJS="$(which osmjs)"
[[ -x "$OSMJS" ]] || {
  OSMJS="${0%.sh}-osmjs"
  sudo apt-get install g++ git \
    lib{boost,expat1,gdal1,geos++,icu,osmpbf,protobuf,shp,sparsehash,v8,z}-dev \
    || exit 1
}

[[ -x "$OSMJS" ]] || {
  CO="$(mktemp --directory --dry-run --suffix=.git)"
  git clone -b master https://github.com/joto/osmium.git "$CO"
  make -C "$CO/osmjs" \
    && cp -a "$CO/osmjs/osmjs" "$OSMJS" \
    && rm -rf "$CO"
}

[[ -x "$OSMJS" ]] && {
  wget -O "$OF" "$OVERPASS" --post-data '
    relation
      ["type"="boundary"]
      ["boundary"="administrative"]
      ["admin_level"="4"]
      ["ISO3166-2"~"^DE-"]
      (47.5,6.0,55.0,15.0);
    (
      ._;
      way(r);
      node(w);
    );
    out;'

  "$OSMJS" -m -l vector -j <(echo '
    Osmium.Callbacks.area = function() {
      print( "INSERT INTO states VALUES(NULL, \"" +
             this.tags.name +
             "\", \"osm\", GeomFromText(\"" +
             this.geom.toWKT() +
             "\", 4326));"
      );
    }'
  ) "$OF" >> "$QF"
}


cat <<EOF >> "$QF"
CREATE TABLE IF NOT EXISTS mtbmbr(
  id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT,
  nr TEXT,
  name TEXT
);
SELECT AddGeometryColumn(
  'mtbmbr', 'g', 4326, 'POLYGON', 'XY'
);
SELECT CreateSpatialIndex(
  'mtbmbr', 'g'
);
DELETE FROM mtbmbr;
EOF


[[ -f "$NAD" ]] || wget -O "$NAD" "$NADURL"


latd=56 latm=0 latdb=55 latmb=54
for r in {01..89}
do
  if [[ "$latm" -eq 0 ]]
  then
    latd="$((latd-1))"
    latm=54
  else
    latm="$((latm-6))"
  fi

  if [[ "$latm" -eq 0 ]]
  then
    latdb="$((latd-1))"
    latmb=54
  else
    latmb="$((latm-6))"
  fi

  lond=5 lonm=40 londr=5 lonmr=50
  for c in {01..59}
  do
    if [[ "$lonm" -eq 50 ]]
    then
      lond="$((lond+1))"
      lonm=0
    else
      lonm="$((lonm+10))"
    fi

    if [[ "$lonm" -eq 50 ]]
    then
      londr="$((lond+1))"
      lonmr=0
    else
      lonmr="$((lonm+10))"
    fi

    ul="${lond}d${lonm}'E ${latd}d${latm}'N"
    ur="${londr}d${lonmr}'E ${latd}d${latm}'N"
    lr="${londr}d${lonmr}'E ${latdb}d${latmb}'N"
    ll="${lond}d${lonm}'E ${latdb}d${latmb}'N"
    mbr=$(
      # nadgrids only support translation up to column 59
      echo -e "$ul \n$ur \n$lr \n$ll \n$ul" |\
      cs2cs -f "%.12f" +proj=lonlat +ellps=bessel +nadgrids="$NAD" \
        +to +init=epsg:4326 |\
      while read lon lat h
      do echo -n "$lon $lat, "
      done
    )

    cat <<EOF >> "$QF"
    INSERT INTO mtbmbr VALUES(NULL, "$r$c", "",
      GeomFromText("POLYGON(( ${mbr%,*} ))", 4326)
    );
EOF
  done
done


for view in grenzmtb=Overlaps mtb2state=Intersects
do
  cat <<EOF >> "$QF"
  CREATE VIEW IF NOT EXISTS ${view%=*} AS
  SELECT s.name Bundesland, m.nr MTB_Nr, m.g MTB_WGS84_Geom,
    Area(Intersection(s.g, m.g))*100/Area(m.g) AS Prozentanteil
  FROM states s, mtbmbr m
  WHERE s.name NOT LIKE "undef%"
    AND s.source="osm"
    AND ${view#*=}(s.g, m.g)
    AND s.ROWID IN (
      SELECT ROWID
      FROM SpatialIndex
      WHERE f_table_name = "states"
        AND search_frame = m.g
    )
  ;
EOF
done


spatialite "$DB" <<EOF
.echo ON
.stats ON
.read "$QF" UTF-8
.quit
EOF


# cleanup
#echo "Keep temporary files ${QF} and ${OF}? [y/N]"
#read ans
#[[ "${ans,}" != "y" ]] && rm -rf "$QF" "$OF"
rm -rf "$QF" "$OF"

} # END DB creation
# ----------------------------------------------------------------


# do a test query if nothing more specific was passed as arguments
GEOMTOUSE="${1:-GeomFromText('POINT(12 51)', 4326)}"
MATCHMODE="${2:-Contains}"

echo "Running some examplary queries.."
echo "This may take some time, safely abort this process with Ctrl+C."
spatialite "$DB" <<EOF |& cat
.headers ON
.stats ON
.echo ON
SELECT AsText($GEOMTOUSE) AS coord, name AS lies_within, source AS says_source
FROM states
WHERE $MATCHMODE(g, $GEOMTOUSE)
  AND id IN (
    SELECT ROWID FROM SpatialIndex
    WHERE f_table_name = 'states'
      AND search_frame = $GEOMTOUSE
  )
;
SELECT *
FROM mtb2state
WHERE Bundesland="Berlin"
  AND Prozentanteil < 30
;
SELECT *
FROM (
  SELECT * FROM grenzmtb WHERE Bundesland="Saarland"
) a
LEFT JOIN (
  SELECT * FROM grenzmtb WHERE Bundesland LIKE "Rheinl%"
) b
ON (
  a.MTB_Nr=b.MTB_Nr
)
;
.quit
EOF


echo "Dumping mtb2state mapping to $MTBMAP."
echo "This may take some time, safely abort this process with Ctrl+C."
spatialite "$DB" <<EOF > "$MTBMAP"
SELECT MTB_Nr, Bundesland FROM mtb2state;
.quit
EOF
echo "Bye."
Retrieved from "https://wiki.openstreetmap.org/w/index.php?title=DE:SpatiaLite/Beispiele/Messtischblatt-TK25-Nummern-Bundesländern-zuordnen&oldid=1049336"