On the fly Darstellung von Grenzlängen und Richtungswinkeln
Aus kvwmap
--Rainer Gareis 12:19, 26. Aug 2010 (CEST) Darstellung der Grenzlängen und Winkel in der Liegenschaftskarte Um in einem Tooltip-fenster die Grenzlänge(Gl) und den Richtungswinkel(Rw) zu sehen, muß man eine neue Tabelle mit oid und Geometrie erzeugen. In dieser Tabelle steht dann die Flurstücksgeometrie aufgeschlüsselt als Line und nicht mehr als Multipoygon.
Tabelle erzeugen:
CREATE TABLE grenzabschnitte2( the_geom geometry, CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2), CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'LINESTRING'::text OR the_geom IS NULL), CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 2398) ) WITH (OIDS=TRUE);
um die Tabelle füllen zu können benötigt man eine neue Funktion "linefrompoly"
CREATE OR REPLACE FUNCTION linefrompoly(geometry)
RETURNS geometry AS
$BODY$SELECT
geomfromtext(
replace(
replace(
replace(
replace(
replace(
asText($1),'MULTIPOLYGON','MULTILINESTRING'
),'POLYGON','MULTILINESTRING'
), '(((', '(('
), ')))', '))'
), ')),((', '),('
), srid($1)
)$BODY$
LANGUAGE 'sql' IMMUTABLE STRICT;
mit folgendem Befehl wird dann die Tabelle gefüllt.
insert into grenzabschnitte2 SELECT makeline(pointn(foo.linestring, foo.count1), pointn(foo.linestring, foo.count1 + 1)) AS the_geom
FROM ( SELECT generate_series(1, npoints(foofoo.linestring)) AS count1, foofoo.linestring
FROM ( SELECT linefrompoly(alkobj_e_fla.the_geom) AS linestring
FROM alkobj_e_fla where folie = '001') foofoo) foo
WHERE (foo.count1 + 1) <= npoints(foo.linestring)
Gist-Index erzeugen, damit die räumliche Abfrage auch schnell durchgeführt wird:
CREATE INDEX grenzabschnitte_gist ON grenzabschnitte2 USING gist (the_geom);
(Dauer ca. 15 Minuten für einen Kompletten Landkreis)
Jetzt kann man sich noch eine Sicht erzeugen, in der die Länge und der Richtingswinkel berechnet wird.
CREATE OR REPLACE VIEW grenzabschnitte AS
SELECT grenzabschnitte2.oid, grenzabschnitte2.the_geom, round(length(grenzabschnitte2.the_geom)::numeric, 2) AS laenge,
CASE
WHEN ((y(endpoint(grenzabschnitte2.the_geom)) - y(startpoint(grenzabschnitte2.the_geom))) = 0 AND (x(endpoint(grenzabschnitte2.the_geom)) > x(startpoint(grenzabschnitte2.the_geom))))THEN '100.0000 gon'::text
WHEN ((y(endpoint(grenzabschnitte2.the_geom)) - y(startpoint(grenzabschnitte2.the_geom))) = 0 AND (x(endpoint(grenzabschnitte2.the_geom)) < x(startpoint(grenzabschnitte2.the_geom))))THEN '300.0000 gon'::text
WHEN round((degrees(atan((x(endpoint(grenzabschnitte2.the_geom)) - x(startpoint(grenzabschnitte2.the_geom))) / (y(endpoint(grenzabschnitte2.the_geom)) - y(startpoint(grenzabschnitte2.the_geom))))) * 1.1111111::double precision)::numeric, 4)::double precision < 0::double precision THEN (round((degrees(atan((x(endpoint(grenzabschnitte2.the_geom)) - x(startpoint(grenzabschnitte2.the_geom))) / (y(endpoint(grenzabschnitte2.the_geom)) - y(startpoint(grenzabschnitte2.the_geom))))) * 1.1111111::double precision)::numeric, 4)::double precision + 400::double precision) || ' gon'::text
ELSE round((degrees(atan((x(endpoint(grenzabschnitte2.the_geom)) - x(startpoint(grenzabschnitte2.the_geom))) / (y(endpoint(grenzabschnitte2.the_geom)) - y(startpoint(grenzabschnitte2.the_geom))))) * 1.1111111::double precision)::numeric, 4)::double precision::text || ' gon'::text
END AS winkel
FROM grenzabschnitte2;
Nun benötigt man nur noch einen abfragbaren Layer ohne Klassen.
Mit freundlicher Unterstützung durch STEFAN RAHN ist dieses Kunststück gelungen.
--Markus Hentschel 08:29, 5. Okt 2010 (CEST) Da zwei benachbarte Flurstücke auch zwei gemeinsame Grenzabschnitte mitbringen, muss man die doppelten noch löschen. Dafür bietet sich folgender Befehl an:
delete from grenzabschnitte2 where oid in ( select a.oid from grenzabschnitte2 a left join grenzabschnitte2 b on length(a.the_geom)=length(b.the_geom) and a.the_geom=b.the_geom and x(pointn(a.the_geom,1)) > x(pointn(b.the_geom,1)) )
Problem: Es gibt zu viele Linien, die Abfrage dauert unendlich. Lösung: ein Script. In diesem Script kann dann auch das Füllen der Tabelle passieren. Das Script geht kilometerquadratweise vor und löscht die doppelten Linien immer nur in diesem Bereich, wodurch die Abfragezeit auf der Datenbank deutlich verringert wird.
Dazu muss eine weitere Tabelle angelegt werden, die ebenfalls nur ein Geometriefeld enthält. Angenommen, sie heißt grenzabschnitte1. Dann sieht das Script so aus:
#!/bin/ksh
# Füllt die Tabelle grenzabschnitte2 mit den Grenzlängen aus den Flurstücksgeometrien
cd /usr/local/pgsql/bin
LOGFILE=/pfad/zu/load_grenzabschnitte.log; export LOGFILE
DBNAME=kvwmapsp; export DBNAME
DBUSER=kvwmap; export DBUSER
# Definition des zu kachelnden Gebiets im Ziel-Koordinatenreferenzsystem S42/83
min_rechts=4517000 # Rechtswert der Ecke links unten
min_hoch=5982000 # Hochwert der Ecke links unten
max_rechts=70000 # maximale Ausdehnung des Bereichs in Metern (min_rechts + max_rechts = Rechtswert rechts oben)
max_hoch=59000 # maximale Ausdehnung des Bereichs in Metern (min_hoch + max_hoch = Hochwert rechts oben)
kachel=1000 # Größe der Kachel in Metern
echo "############## Beginn `date +%c` ################" >> $LOGFILE 2>&1
# Leeren und neues Füllen der temporären Tabelle
./psql -d $DBNAME -U $DBUSER -c "
-- Vorhandene Indices löschen
DROP INDEX g1ixoid;
DROP INDEX g1ixgeom;
TRUNCATE grenzabschnitte1;
-- temporäre Tabelle füllen
INSERT INTO grenzabschnitte1
SELECT makeline(pointn(foo.linestring, foo.count1), pointn(foo.linestring, foo.count1 + 1)) AS the_geom
FROM ( SELECT generate_series(1, npoints(foofoo.linestring)) AS count1, foofoo.linestring
FROM ( SELECT linefrompoly(alkobj_e_fla.the_geom) AS linestring
FROM alkobj_e_fla where folie = '001') foofoo) foo
WHERE (foo.count1 + 1) <= npoints(foo.linestring);
" # Ende SQL
if test $? -eq 0
then
echo "`date +%c` grenzabschnitte1 aktualisiert" >> $LOGFILE 2>&1
# Tabelle analysieren für einen schnellen Zugriff
./psql -d $DBNAME -U $DBUSER -c "
VACUUM ANALYZE grenzabschnitte1;
" # Ende SQL
./psql -d $DBNAME -U $DBUSER -c "
-- Indices wieder aufbauen für einen noch schnelleren Zugriff
CREATE INDEX g1ixoid ON grenzabschnitte1 (oid);
CREATE INDEX g1ixgeom ON grenzabschnitte1 USING gist (the_geom);
" # Ende SQL
# Lösche die doppelt vorhandenen wieder raus
# portionsweiseweise, weil sonst die Prüfung zu lange dauert
# äußere Schleife - erst die Zeilen von unten nach oben ...
i=0
while [ $i -lt $max_hoch + $kachel ]
do
z2=$(($min_hoch + $i))
# innere Schleife - ... dann die Kacheln von links nach rechts
j=0
while [ $j -lt $max_rechts + $kachel ]
do
z1=$(($min_rechts + $j))
z3=`expr $z1 + $kachel`
z4=`expr $z2 + $kachel`
echo $z1 " " $z2 ", " $z3 " " $z4
sql="delete from grenzabschnitte1 where oid in (select a.oid from grenzabschnitte1 a left join grenzabschnitte1 b on length(a.the_geom)=length(b.the_geom) and a.the_geom=b.the_geom and x(pointn(a.the_geom,1)) > x(pointn(b.the_geom,1)) where a.the_geom && GeomFromText('POLYGON((${z1} ${z2},${z3} ${z2},${z3} ${z4},${z1} ${z4},${z1} ${z2}))', 2398) and b.the_geom && GeomFromText('POLYGON((${z1} ${z2},${z3} ${z2},${z3} ${z4},${z1} ${z4},${z1} ${z2}))', 2398));"
#echo $sql
./psql -d $DBNAME -U $DBUSER -c "$sql"
j=`expr $j + $kachel`
done # Ende innere Schleife
i=`expr $i + $kachel`
done # Ende äußere Schleife
echo "`date +%c` grenzabschnitte1 auf doppelte Grenzen geprüft" >> $LOGFILE 2>&1
# Übertragung der Grenzlängen in die endgültige Tabelle
./psql -d $DBNAME -U $DBUSER -c "
-- Starten der Transaktion
BEGIN;
-- Index entfernen
DROP INDEX g2ixgeom;
-- Leeren der Tabelle
TRUNCATE grenzabschnitte2;
-- Füllen der Tabelle
INSERT INTO grenzabschnitte2 SELECT * FROM grenzabschnitte1;
-- Beenden der Transaktion
END;
" # Ende SQL
if test $? -eq 0
then
echo " grenzabschnitte2 erfolgreich aktualisiert" >> $LOGFILE 2>&1
else
echo ">> grenzabschnitte2 konnte nicht aktualisiert werden" >> $LOGFILE 2>&1
fi
# Tabelle analysieren
./psql -d $DBNAME -U $DBUSER -c "
VACUUM ANALYZE grenzabschnitte2;
" # Ende SQL
if test $? -eq 0
then
echo " VACUUM erfolgreich durchgeführt" >> $LOGFILE 2>&1
else
echo ">> VACUUM konnte nicht durchgeführt werden" >> $LOGFILE 2>&1
fi
# Tabelle indizieren
./psql -d $DBNAME -U $DBUSER -c "
CREATE INDEX g2ixgeom
ON grenzabschnitte2
USING gist
(the_geom);
" # Ende SQL
if test $? -eq 0
then
echo " INDEX erfolgreich angelegt" >> $LOGFILE 2>&1
else
echo ">> INDEX konnte nicht angelegt werden" >> $LOGFILE 2>&1
fi
else
echo "`date +%c`>> grenzabschnitte1 konnte nicht aktualisiert werden" >> $LOGFILE 2>&1
echo "`date +%c`>> grenzabschnitte2 behält den alten Stand" >> $LOGFILE 2>&1
fi
echo " " >> $LOGFILE 2>&1
echo "############## Ende `date +%c` ################" >> $LOGFILE 2>&1
echo " " >> $LOGFILE 2>&1
exit 0
Das Script benötigt auf meiner Maschine etwas unter einer halben Stunde bei 2,5 Millionen Ausgangslinien.
