SpatiaLite/Examples/Mapping SurveySheets to FederalState names

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

Intro

This script builds a SpatiaLite database from scratch containing two tables:

  • states with columns name, source, g
  • mtbmbr with columnr name (just a placeholder), nr, g

It was build with the goal in mind to determine which federal state territory is covered by which survey map sheet(s) of the TK 25 or resp. Messtischblätter map series. The map sheets in the series are numbered uniquely by a grid scheme which is leveraged by the script to calculate, rather than import, the bounding boxes of these map sheets. It's currently useful to

  • determine all map sheets covering a certain federal state completely, and, vice versa:
  • determine the federal state(s) a map sheet has information about
  • retrieve all map sheet numbers covering the border region of a federal state (grenzmtb view)
  • get map sheets that cover the border region along two (or more) certain adjacent federal states
  • calculate the difference between the polygon areas and the boundary as entered in OSM
  • ...

states is filled twice, using

  • rough geometry contained in Geofabrik's polygon files (source=poly)
  • precise boundary geometry pulled from Overpass API (source=osm)

mtbmbr is filled with minimum bounding rectangles of Messtischblatt / TK25 map rectangles

  • bbox'es are synthesized / computed by the script
  • the table totals 89x59 entries, some cover ocean or territory a paper map never was published for

The code

  • was written and tested on Kubuntu 14.04 LTS Live/Install ISO booted off an USB key, corei3 CPU, 3G RAM
  • pulls needed deps using "sudo apt-get", among them spatialite
  • needs write access within /tmp and the directory the script is put in
  • needs a recent osmjs, if it is uninstalled the script makes an attempt to build it from source
    • the deb distributed with trusty was not stable and ran out of memory too quick, even with -l sparsetable.
    • since the live system resides in memory, -l vector is used to constrain osmjs to a lower memory footprint.
    • the result of a potential build is cached, consecutive runs of the script will not rebuild osmjs

Some exemplary queries are found near the end of the script. Feel free to reuse this code or the db it produces for your own purposes.

Produced files

  • the sqlite database to work with later on, containing
    • the tables with geometry column and spatial index
    • saved views that make use of the spatial index
  • an mtb to state map which is a dump of the view mtb2state

See also

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=SpatiaLite/Examples/Mapping_SurveySheets_to_FederalState_names&oldid=1049329"