
'FOSSGIS 2022
SQL für QGIS und SpatiaLite.
Verschiedene Übungen 
'
'INFOS IM NETZ'
https://sqlitebrowser.org/dl/ --sqlitebrowser
http://www.gaia-gis.it/gaia-sins/windows-bin-amd64/mod_spatialite-5.0.1-win-amd64.7z --Download Spatialite Library für sqlite3
https://github.com/sqlitebrowser/sqlitebrowser/wiki/SpatiaLite-on-Windows  -- Anbindung von Spatialite unter Windows

http://www.gaia-gis.it/gaia-sins/index.html  					-- Spatialite-Webseite
http://www.gaia-gis.it/gaia-sins/spatialite-sql-5.0.1.html			-- Spatialite-SQL-Funktionen

http://www.sqlite.org/lang.html							-- SQLite-Sql-Syntax
http://northredoubt.com/n/2012/10/19/spatialite-and-triggers-to-update-data/	-- Trigger in SpatiaLite

http://www.gaia-gis.it/gaia-sins/spatialite-gui-sources/spatialite_gui-2.1.0-beta1.zip	-- Download SpatiaLite-Gui

https://www.surfaces.co.il/spatialite-speedup-your-query-with-spatial-indexing/	--Zu räumlichen Idexen
https://www.gaia-gis.it/fossil/libspatialite/wiki?name=SpatialIndex		-- Zu räumlichen Indexen

https://sqlitebrowser.org/dl/ --sqlitebrowser
http://www.gaia-gis.it/gaia-sins/windows-bin-amd64/mod_spatialite-5.0.1-win-amd64.7z --Download Spatialite Library für sqlite2
-----



'LOS GEHTS!'
--------------------------------------------
'SÄMTLICHE SHAPES AUS gemeinden IN EINE SPATIALITE-DATENBANK SCHREIBEN'
'DEN TEXT IN EINE EIGENE DATEI KOPIEREN UND ALS shp2sqlite in das Verzeichnis bin im Nutzerverzeichnis speichern'
--Mit diesem Skript wird aus den Shapefiles im NRW-Ordner eine SpatiaLite-Datenbank erzeugt
#!/bin/bash
# shp2sqlite-Shapes zu SQLITE
for shps in *.shp; do
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS OFF  --config OGR_SQLITE_CACHE 8192 -f "SQLite" -dsco SPATIALITE=yes -lco SRID=25832 -nlt PROMOTE_TO_MULTI -skipfailures -append -gt unlimited -a_srs EPSG:25832 ziel.sqlite $shps
done
#csv-Dateien
for csvs in *.csv; do
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS OFF  --config OGR_SQLITE_CACHE 8192 -f "SQLite" -update -append -gt unlimited-skipfailures ziel.sqlite $csvs
done
# Excel-Dateien
for xlsxs in *.xlsx; do
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS OFF  --config OGR_SQLITE_CACHE 8192 -f "SQLite" -update -append -gt unlimited -skipfailures ziel.sqlite $xlsxs
done
---------------------------------------------------------------------------------------
-- Windows-Version des gleichen Skripts
# 'shp2sqlite.bat'
for %%x in (*.shp) do (
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS OFF  --config OGR_SQLITE_CACHE 8192 -f "SQLite" -dsco SPATIALITE=yes -lco SRID=25832 -nlt PROMOTE_TO_MULTI -skipfailures -append -gt unlimited -a_srs EPSG:25832 ziel.sqlite %%xx
)
#csv-Dateien
for %%x in (*.csv) do (
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS OFF  --config OGR_SQLITE_CACHE 8192 -f "SQLite" -update -append -gt unlimited-skipfailures ziel.sqlite %%x
)
# Excel-Dateien
for %%x in (*.xlsx) do (
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS OFF  --config OGR_SQLITE_CACHE 8192 -f "SQLite" -update -append -gt unlimited -skipfailures ziel.sqlite %%x
)
------------------------------------------------------------

'Neue Layer in eine SpatiaLite-Datenbank importieren:'
--Über den Speichern als-Dialog lassen sich Layer in eine SpatiaLite-Datenbank schreiben.
-- Weitere Möglichkeit DB-Verwaltung: Erlaubt auch Umprojizieren beim Import
Datenbank > DB-Verwaltung > DB-Verwaltung
Dort: Tabelle > Layer importieren
Primärschlüssel und Codierung utf.8 wählen
Räumlichen Index anlegen
Umprojizieren wenn möglich.

'SPATIALITE_DATENBANK iIN QGISANMELDEN'
Layer > Layer hinzufügen > SpatiaLite-Layer  hinzufügen
NEU klicken und Datenbank ziel.sqlite auswählen
Geht auch im Datenbankmanager: 
Rechtsklick auf SpatiaLite im  Datenbankbaum > Neue Verbindung


------
'EINGEBAMÖGLICHKEITEN FÜR SQL'
-- Virtuelle Layer
Layer > Layer hinzufügen > Virtuelle Layer
Zugriff auf alle geladenen Layer, virtueller layer wird dynamisch im Projekt erzeugt
Langsam! SpatiaLite-SQL und alle Funktionen des Ausdruckseditors sind verfügbar.

---Verarbeitungs-Toolbox
 Vektoren Allgemein > SQL-Anweisung ausführen
 Datenbank > spatialite-SQL ausführen

GDAL> Vektor Verschiedenes > SQL-Anweisung ausführen
basiert auf ogr2ogr -sehr schnell

---- Datenbankmanager
Datenbank > DB-Verwaltung > DB-Verwaltung
Dort SQL-Fenster öffnen
Über die kleine SQL-Schaltfläche ist eine grafische Hilfe verfügbar


-- SQLite-Studio und SQLitebrowser
Sehr gut gestaltete GUI für SQLite

https://sqlitebrowser.org/dl/

https://sqlitestudio.pl/index.rvt?act=download


Einbinden der SpatiaLite-Moduls nowendig!



'Lektion 1: LANDKREISE AUS GEMEINDEN GENERIEREN'

'Inhalt anzeigen'
-- In der DB-Verwaltung das SQL-Fenster aufrufen
-- Datenbank > SQL-Fenster (F2)
-- In der DB-Verwaltung und in QspatiaLite ist jeweils nur eine Abfrage ausführbar
-- Kein Skript, mit durch ; getrennten Abfragen.


-- Abfrage eingeben
select              -- Auswahl aus einer Tabelle
*                   -- sämtliche Spalten werden angezeigt
from gemeinden      -- die Tabelle von der gewählt wird

--Bei der Abfrage einzelne Spalten angeben und über as mit einem Alias versehen
select 
rs as regionalschluessel, --Spalte rs wird als regionalschluessel angezeigt
gen as name               --Spalte gen wird als name angezeigt
from gemeinden
-- Einwohnertabelle anzeigen
select
*
from ewz

'EINWOHNER AN DIE GEMEINDEN ANBINDEN'
-- Abfrage
select
h.geom as geom,
h.rs as rs,
h.gen as gen,
e.ewz as ewz
from gemeinden as h join ewz as e on (e.rs = h.rs)  
-- h.geom enthält die Polygon-Geometrie der Gemeinden // In Datenbanken gibt es Geometriespalten.
-- gemeinden bekommt den Alias h und ewz den Alias e.
-- h.rs bedeutet die Spalte rs der tabelle h, was der Alias von gemeinden ist.
-- Mit join wird die ewz-Tabelle über ein gemeinsames Schlüsselfeld angebunden
-- on (e.rs = h.rs)  verbindet beide Tabellen über die Schlüsselspalten h.rs und e.rs
-- In beiden Tabellen gibt es den Regionalschlüssel


'NEUER LAYER MIT EWZ ERZEUGEN'
-- Mit create table wirs aus der Abfrage eine neue Rabelle
create table gemeinden_ewz as --erzeugt eine neue Tabelle gemeinden_ewz aus der vorherigen Abfrage
select
casttomultipolygon(g.geom) as geom,
g.rs as rs,
g.gen as gen,
e.ewz as ewz
from gemeinden as g left join ewz as e on (e.rs = g.rs);

--In Spatialite muss die Geometrie in die Metadaten geschrieben werden
-- Dies lässt sich mit folgendem Aufruf umsetzen

Select RecoverGeometryColumn('gemeinden_ewz','geom',25833,'MULTIPOLYGON');  --legt die Geometriespalte an

--Im folgenden wird noch ein räumlicher Index angelegt, der ein schnelleres Verarbeiten des neuen Layers erlaubt.

Select CreateSpatialIndex('gemeinden_ewz','geom');  --erzeugt einen räumlichen Index


--------

'Die geometrielose Landkreistabelle mit den Landkreisnamen soll an die gemeinden gejoint werden.'
--INHALT VON GEMEINDEN UND LANDKREISEN ANZEIGEN'
select
*
from lkr
--dann
select
*
from gemeinden
------------------------------------

'KREISID UND REGIONALSCHLÜSSEL mit UNION VERGLEICHEN'
-- Beide Teile der UNION-Abfrage müssen die gleichen  Spalten haben, darum über Alias-Bildung umbennen!
-- Der Vergleich wird über where auf die kreisfreie Stadt Marburg beschränkt
-- Die Where-Klausel dient der Auswahl bestimmter Zeilen und entspricht dem, was im QGIS-Ausdruckseditor passiert.
select
'Kreise' as tabelle,		   -- Bezeichnung der Ursprungstabelle ueber text 'Kreise'
krid as id,                        -- krid wird als Spalte id ausgegeben
name as name                        -- name wird als Spalte Name ausgeben
from lkr                           -- gewählt von der Tabelle Landkreise  
where name like '%Marb%'            -- Nur Zeilen mit der Zeichenkette Marburg in der Spalte name

union -- Beide Abfragen werden untereinandergeschrieben

select
'Gemeinden' as tabelle,		   -- Bezeichnung der Ursprungstabelle ueber text 'Gemeinden'
rs as id,                          -- rs wird als Spalte id ausgegeben
gen as name                        -- gen wird als Spalte Name ausgebe
from gemeinden                     -- gewählt von der Tabelle gemeinden 
where gen like '%Marburg%'            -- Nur Zeilen mit der Zeichenkette Marburg in der Spalte gen

order by id desc  --Absteigende Sortierung nach id

'ERGEBNIS:'
--krid aus lkr besteht aus dem 2. bis 5. Zeichen des Regionalschlüssel rs
--------------------------------------------------------------


'Mit der Funktion substr(spaltenname,StelleDesErstenZeichen,AnzahlZeichen) KREISID AUS GeMEINDE ANZEIGEN'
select
gen,
substr(rs,2,4) as krid -- aus der Spalte rs werden jeweils ab dem zweiten Zeichen vier Zeichen in eine neue Spalte krid geschrieben
from gemeinden
order by krid;
--substr gibt es auch als Funktion im Ausdruckseditor

'GRUPPIEREN NACH LANDKREISBEZUG, damit es nur eine Zeile je kreisid gibt: group by steht immer hinter from und where'
select
substr(rs,2,4) as krid
from gemeinden
group by substr(rs,2,4) -- Alle Zeilen mit dem gleichen krid werden zusammengefasst es gibt nur eine Zeile je landkreis

'WEITERE ZEILEN AUS DER GEMEINDETABELLE MITNEHMEN UND AGGREGIEREN UM DIE ANDEREN SPALTEN SINNVOLL ZUSAMMENZUFASSEN'
select
substr(rs,2,4) as krid,	
st_union(geom) AS geom,        -- st_union (oder gunion) verschmilzt die Geometrien der gruppierten Zeilen also DISSOLVE!
astext(st_union(geom)) AS WKT,     -- noch einmal als WKT-text ausgeben, um den Geometrietyp zu prüfen 
count(rs) as anzahl,                   -- die anzahl der zusammengefassten rs-Zeilen wird gezählt (Anzahl der Gemeinden im LKR)
sum(ewz) as ewz_kreis,		       -- Die Einwohnerzahl der Gemeinden wird summiert
group_concat(gen,'; ') as gem_liste   -- Die Gemeindenamen aus der Spalte gen werden als Liste mit dem Trennzeichen Semikolon geschrieben                  
from gemeinden_ewz 		       -- es wird von gemeinden_ewz, der Tabelle mit den gejointen Einwohner gewählt!  
group by substr(rs,2,4) 	       -- Alle Zeilen mit dem gleichen krid werden zusammengefasst es gibt nur eine Zeile je landkreis

--DEN GEOMETRIETYP GENERELL AUF MULTIPOLYGON SETZEN

select
casttomulti(st_union(geom)) AS geom, --- casttomultipolygon sorgt für durchgehend Multipolygon-Geometrien
astext(casttomulti(st_union(geom))) AS WKT,
count(rs) as anzahl,
group_concat(gen,'; ') as gem_liste,
sum(ewz) as ewz_kreis,
substr(rs,2,4) as krid
from gemeinden_ewz
group by substr(rs,2,4)
--LANDKREISNAMEN ANJOINEN

select
casttomulti(st_union(geom)) as geom,
substr(g.rs,1,5) as krid,
l.name as name,
count(g.gen) as anzahl,
sum(g.ewz) as ewz,
group_concat(g.gen,'; ') as liste
from gemeinden_ewz as g join lkr as l on (substr(g.rs,1,5) = l.krid)   -- die Landkreistabelle wird über das neue krid-feld angehängt
group by substr(rs,2,4)

--Hiermit haben wir eine vollständige DISSOLVE-FUNKTION mit Agreggierung von Spalten.

create table landkreise as
select
casttomulti(st_union(geom)) as geom,
substr(g.rs,1,5) as krid,
l.name as name,
count(g.gen) as anzahl,
sum(g.ewz) as ewz,
group_concat(g.gen,'; ') as liste
from gemeinden_ewz as g join lkr as l on (substr(g.rs,1,5) = l.krid)  -- die Landkreistabelle wird angehängt
group by substr(rs,2,4); -- Alle Zeilen mit dem gleichen krid werden zusammengefasst es gibt nur eine Zeile je landkreis


Select RecoverGeometryColumn'landkreise','geom',25833,'MULTIPOLYGON');  --legt die Geometriespalte an
Select CreateSpatialIndex('landkreise','geom');  --erzeugt einen räumlichen Index



-- Die ganze Abfrage lässt sich auch ogr2ogr übergeben
ogr2ogr --config OGR_SQLITE_SYNCHRONOUS off --config OGR_SQLITE_CACHE 8192  -f SQLite -dsco SPATIALITE=yes -lco SRID=25833 -nln kreise -skipfailures -gt unlimited -a_srs EPSG:25832 datei.sqlite \
-sql"create table landkreise as select
casttomulti(st_union(g.geom)) AS geom, 
count(g.rs) as anzahl,
group_concat(g.gen,'; ') as gem_liste,
sum(g.ewz) as ewz_kreis,
l.name as name,
substr(rs,2,4) as krid
from gemeinden_ewz as g inner join lkr as l on (l.krid = substr(g.rs,1,5)) -- die Landkreistabelle wird angehängt
group by substr(rs,2,4)"


select 'AUFGABE WIE IM MODELLER-KURS'
/* Waldanteil und Waldfläache je EW ermitteln 
in mehreren Abfragen hintereinander  */

-- Einwohnertabelle an Gemeinden anbinden, über einen Join an den Regionalschlüssel
create table gemewz as
Select
g.geom AS geom,
g.RS AS RS,
g.GEN AS GEN,
e.EWZ AS EWZ
from
gemeinden as g  inner join EWZ as e ON (e.RS = g.RS) -- Einwohner an gemeinden anjoinen
;

--Wald und Gemeindeflächen Verschneiden, um Waldflchen den gemeinden zuordnen zu können
create table intersectum as
Select
    g.RS AS RS,
    g.GEN AS GEN,
    ST_Intersection(g.geom,d.geom) AS geom -- Intersect Geometrie-verschneidung, um die gemeinsamen Wald-Gemeinde-Geometrien zu ermitteln.
    from
    gemeinden as g  left join dlm250 as d  ON intersects(g.geom,d.geom) -- Raeumliche Abfrage: nur Gemeinden verschneiden, die sich berühren
    where  d.nutzung = 'AX_Wald'   -- Nur objekte aus dem DLM nutzen, die tatschlich Wald sind
         and d.ROWID IN
            (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=g.geom)  -- Räumlichen Index benutzen (optional)
;

-- Wald-Gemeinde-Geometrien nach regionalschlüssel verschmelzen, um die Gesamtwaldfläche jeder Gemeinde zu ermitteln
-- Danach gibt es jeweils eine Multipart-Geometrie des Waldes je Gemeinde
create table wdissolve as
select
RS,
(st_union(geom)) as geom -- st_union ist dissolve in kombination mit group by, hier nach RS gruppieren und die Geometrie verschmelzen
from
intersectum
group by RS;

-- Ermittelte Wald-Gemeinde-Flächen anjoinen und für die Berechnung des Waldanteils nutzen

create table waaanteil as
Select
g.geom AS geom, --die Gemeindegeometrie
g.RS AS RS, 		--Regionalschlüssel
g.GEN AS GEN, 		-- Gemeindename
g.EWZ AS EWZ,		-- Einwohnerzahl
round((area(i.geom) / 1000000),2) AS WFL, -- Die Waldfläche wird berechnet aus der wdissolve-Geometrie
round((area(i.geom)  / (area(g.geom))* 100)) AS anteil, -- Der Waldanteil wird aus der Wald- und Gesamtgeometrie berechnet
round(area(i.geom)  / g.EWZ)  AS WF_EW		-- Waldfläche je Einwohner wird berechnet.
from wdissolve as i left join gemewz as g on (i.RS = g.RS);

Select RecoverGeometryColumn("waaanteil",'geom',25833,'MULTIPOLYGON');
Select CreateSpatialIndex("waaanteil",'geom');

SELECT UpdateLayerStatistics('waaanteil');
-- Hilfstabellen löschen
select droptable(NULL,'gemewz');
select droptable(NULL,'intersectum');
select droptable(NULL,'wdissolve');

'ALTERNATIVE MIT VERSCHACHTELTET ABFRAGE'
/* Waldanteil und Waldfläache je EW ermitteln 
verschachtelt  */
-- Hier sind alle Abfragen in verschachtelten Funktionen zusammengefasst
---Verschachtelete Funktionen mit räumlichen Index
create table waaanteil as
Select
g.geom AS geom,
g.RS AS RS,
g.GEN AS GEN,
e.EWZ AS EWZ,
round((area(st_union(ST_Intersection(g.geom,d.geom))) / 1000000),2) AS WFL,
round((area(st_union(ST_Intersection(g.geom,d.geom)))) / (area(g.geom))* 100) AS anteil,
round((area(st_union(ST_Intersection(g.geom,d.geom)))) / e.EWZ)  AS WF_EW
from
(gemeinden as g left join dlm250 as d ON intersects(g.geom,d.geom)) left join EWZ as e ON e.RS = g.RS
where d.ROWID IN
            (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=g.geom)
   and
    d.nutzung = 'AX_Wald'
Group by g.RS;

Select RecoverGeometryColumn('waaanteil','geom',25833,'MULTIPOLYGON');
Select CreateSpatialIndex('waaanteil','geom');
select updatelayerstatistics('waaanteil');

-- ohne räumlichen Index
create table waaanteil as
Select
g.geom AS geom,
g.RS AS RS,
g.GEN AS GEN,
e.EWZ AS EWZ,
round((area(st_union(ST_Intersection(g.geom,d.geom))) / 1000000),2) AS WFL,
round((area(st_union(ST_Intersection(g.geom,d.geom)))) / (area(g.geom))* 100) AS anteil,
round((area(st_union(ST_Intersection(g.geom,d.geom)))) / e.EWZ)  AS WF_EW
from
(gemeinden as g left join dlm250 as d ON intersects(g.geom,d.geom)) left join EWZ as e ON e.RS = g.RS
where
    d.nutzung = 'AX_Wald'
Group by g.RS;

---------------
----------------
------------
--
'VERSCHNEIDUNGEN'
-- INTERSECT ZWEIER GEOMETRIEN
create table verschneidung as   --Erzeugend der neuen Tabelle über create table
select                          -- Die Tabelle wird aus der Abfrage erstellt
  casttomultipolygon(st_intersection(a.geom,b.geom)) AS geom, --Die Funktion ST_intersect entspricht der Intersect / Verschneiden-Funktion in GIS-Pprogrammen
  a.rs,								   -- casttomultipolygon ermöglicht eine Multipolygon-Geometrie im Ergebnis
  a.gen,
  a.bez,                                                             -- weitere Spalten können mit abgefragt werden (Aufzählung mit Komma)       
  b.nutzung as nutzung 
  
from gemeinden as a inner join dlm250 as b on st_intersects(a.geom,b.geom) = 1 -- Die Verschneidung wird nur für Geometrien durchgeführt, die sich wirklich schneiden
   where nutzung = 'AX_Wald' and  a.ROWID IN
              (SELECT ROWID FROM SpatialIndex WHERE f_table_name='gemeinden' AND search_frame=b.geom); -- Ein räumlicher Index wird verwendet, dieser Teil ist optional. Bei großen Datenmengen, aber für schnelles Verschneiden unerlässslich
---------------
Select RecoverGeometryColumn('verschneidung','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('verschneidung','geom');
SELECT UpdateLayerStatistics('verschneidung');


/* Abfragen aus der Live-Demo */

--##################################################

drop table  if exists waldschnitt;			-- Vorhande Tabelle löschen

create table waldschnitt as				-- Neue Tabelle erzeugen
	with  verschnitt as							-- Die Unterabfrage bekommt den Namen verschnitt		
      ( --Start Unterabfrage
      select
	  row_number() over() as id,			-- Eindeutige ID für das Ergebnis
      g.RS as RS,							-- Der Reginalschlüssel wird aus der Gemeindetabelle geholt (g ist Alias der Gemeinde)
	  g.gen as gemeinde,					-- gen ist der Gemeindename, umbenannt zu gemeinde
      round((st_area(g.geom)/ 1000000),3) as qkm,  -- Berechnung der Gemeindefläche aus der Gemeindegeometrie g.geom
       round(		-- Runden						-- mit der Funktion ST_area Fläche berechnen, dann Division durch 1000000 für qkm
			sum(				-- Summierung der Fläche sämtlicher Waldschnipsel (Gruppierung nach RS)
				st_area(  		-- Flächenberechnung aus Geometrie
					st_intersection(g.geom,w.geom)  --as geom	-- Dies ist die Intersection-Funktion (Verschneidung von Wald und Gemeinden wie im GIs)
					    ) / 1000000			-- St_area schließen und durch 1000000 dividieren für qkm
					)		                 -- Summierung schließen
			 ,3)                             -- auf dritte Stelle runden und schließen
		  as qkm_wald,			-- Ergebnis ist die Waldfläche, die aus der verschnittenen Gemeinde/Waldgeometrie berechnet wird
		  g.geom as geom		-- Die einfache Gemeindegeometrie
	
        from  gemeinden  as g -- dlm250 as w			             -- Eingabe Gemeindelayer
				left join 					                 -- Die Anbindung eines anderen Layers über eine Beziehung, alle Objekte von Layer 1
				dlm250 as w on st_intersects(g.geom,w.geom)  --Angebunden wird DLM 250, wenn sich beide Geometrien überschneiden.
        where w.nutzung = 'AX_Wald'					        -- Wie im Ausdruckseditor, die Waldflächen werden ausgewählt
                         AND w.ROWID IN                    -- für schnellere Verarbeitung: Räumlichen Index ansprechen
                        (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=g.geom)
		group by g.rs	      					             --Gruppieren nach gemeindeschlüssel rs
		) -- Ende Unterabfrage
select					-- Hauptabfrage, damit es bei der Flächenberechnung übersichtlich bleibt
--*
id, RS, gemeinde, qkm, qkm_wald,
round((ifnull(qkm_wald,0) / qkm) * 100) as anteil,  --Hierwird der prozentuale Anteil berechnet
geom as geom
from verschnitt		-- Eingabe ist das Ergebnis der Unterabfrage verschnitt, die with eingeleitet wurde
;

Select RecoverGeometryColumn('waldschnitt','geom',25832,'MULTIPOLYGON');  -- geom anmelden, damit die Tabelle als Layer im QGIS angezeigt werden kann.
Select CreateSpatialIndex('waldschnitt','geom');
SELECT UpdateLayerStatistics('waldschnitt');
--Select disableSpatialIndex('waldschnitt','geom');

/* Verschneiden und verschmelzen*/
select
g.RS,
g.gen,
-- Dies ist eine Kombi von verschmelzen und verschneiden 
--(Verschneidung mit st_intersection, verschmelzen mit st_union, Geometriesplitter ausschliessen mit st_collectionextract
casttomulti(ST_union(st_collectionextract(st_intersection(g.geom,w.geom),3))) as geom		
from  gemeinden  as g -- dlm250 as w			             -- Eingabe Gemeindelayer
            left join 					                 -- Die Anbindung eines anderen Layers über eine Beziehung, alle Objekte von Layer 1
            dlm250 as w on st_intersects(g.geom,w.geom)  --Angebunden wird DLM 250, wenn sich beide Geometrien überschneiden.
where w.nutzung = 'AX_Wald'					        -- Wie im Ausdruckseditor, die Waldflächen werden ausgewählt
            AND w.ROWID IN                    -- für schnellere Verarbeitung: Räumlichen Index ansprechen
                (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=g.geom)
group by g.rs
;--Gruppieren nach gemeindeschlüssel rs                                                                        
                                                


/* Mit Ermittlung der Walfläche je Einwohner ohne Kommentare */
drop table  if exists wald_ant_mit_ew
create table wald_ant_mit_ew aS
 with  verschnitt as	
      (
      select
      g.RS as RS,
      round((st_area(g.geom)/ 1000000),3) as qkm,
       round(
			sum(st_area(
					st_intersection(g.geom,w.geom) 
                        ) / 1000000	
                ),3)
    as qkm_wald
    
    from  gemeinden  as g left join dlm250 as w on st_intersects(g.geom,w.geom) 
    where w.nutzung = 'AX_Wald'
    AND w.ROWID IN
            (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=g.geom)
    group by g.rs	      					  
    )
select					
a.fid, a.RS, a.gen as gemeinde, b.qkm, b.qkm_wald,
round((ifnull(b.qkm_wald,0) /b.qkm) * 100) as anteil, 
round((b.qkm_wald * 1000000) / c.ewz) as wfqm_ew,
c.ewz as ewz,
a.geom as geom
from gemeinden as a left join verschnitt as b on (a.RS = b.RS)
                    left join ewz as c on (a.rs = c.rs)
;


Select CreateSpatialIndex('wald_ant_mit_ew','geom');
Select RecoverGeometryColumn('wald_ant_mit_ew','geom',25832,'MULTIPOLYGON');
SELECT UpdateLayerStatistics('wald_ant_mit_ew');
--Select disableSpatialIndex('wald_ant_mit_ew','geom');





/* Strukturiert mit zwei Unterabfragen */
with 
   inters as  -- Verschneidung und Flächenberechnung
        (
        select
        a.RS as RS,
        round((st_area(a.geom)/ 1000000),3) as qkm,
        round((st_area(st_intersection(a.geom,b.geom)) / 1000000),3) as qkm_wald

        from gemeinden as a left join dlm250 as b on st_intersects(a.geom,b.geom)
        where b.nutzung = 'AX_Wald'
        	 AND b.ROWID IN
            (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=a.geom)			
        )
        ,
    wafl as   --- Summieren je Gemeinde
        (
        select
        RS,
        qkm,
        sum(qkm_wald) as qkm_wald
        from inters
        group by RS
        )
select
row_number() over() as id,
a.RS as RS,
a.GEN as gemeinde,
b.qkm as qkm,
b.qkm_wald as qkm_wald,
round((b.qkm_wald / b.qkm) * 100) as anteil,
a.geom as geom
from gemeinden as a left join wafl as b on (a.RS = b.RS)
;       
--##################################################            
    

/* geometryverschneidung SQL */

--DIFFERENZ / ERASS ZWEIER GEOMETRIEN IST KOMPLIZIERTER, WEIL ZUM AUSSTANZEN NUR EINZELGEOMETRIEN VERWENDET WERDEN KÖNNEN

-- Schritt 1: aus beiden Layern eine vollstndig verschmolzene Geometrie erzeugen
-----------------
create table diff_vorne			-- Die Tabelle aus der gestanzt wird
as
select
st_union(a.geom) AS geom		-- verschmelzen der ersten Geometrie mit Pufferung zum Lücken schließen
      from gemeinden as a;      -- DIES IST DER URSPRUNGSLAYER A, DER IM LETZTEN SCHRITT DES SKRIPTES NOCH EINMAL EINGETRAGEN WERDEN MUSS

Select RecoverGeometryColumn('diff_vorne','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('diff_vorne','geom');
SELECT UpdateLayerStatistics('diff_vorne');

----------------
create table diff_stanz 		-- Die Tabelle die stanzt
as
select
st_union(a.geom) AS geom		-- verschmelzen der ersten Geometrie mit Pufferung zum Lücken schließen
      from dlm250 as a
where a.nutzung = 'AX_Wald';  -- DIES IST DER URSPRUNGSLAYER B, DER DIE LÖCHER IN DEN LAYER A STANZT

Select RecoverGeometryColumn('diff_stanz','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('diff_stanz','geom');
SELECT UpdateLayerStatistics('diff_stanz');

--Schritt 2: Mit einem Layer ein Loch in den anderen Layer stanzen
-----------------
create table diffdiss as		--diffdiss ist der Ergebnis des Verschneidungsprozesses mmit st_difference		
select
      casttomultipolygon(st_difference(diff_vorne.geom,diff_stanz.geom)) AS geom 			--Der Layer diffstanz schneidet Löcher in die Geometrie von diffvorne
      from diff_vorne inner join diff_stanz on st_intersects(diff_vorne.geom,diff_stanz.geom) = 1
       where  diff_stanz.ROWID IN
              (SELECT ROWID FROM SpatialIndex WHERE f_table_name='diff_stanz' AND search_frame=diff_vorne.geom);

Select RecoverGeometryColumn('diffdiss','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('diffdiss','geom');
SELECT UpdateLayerStatistics('diffdiss');

---Alternative SCHRITT 1 UND SCHRITT 2 DURCH VERSCHACHTELUNG ZUSAMMENFASSEN



create table diffdiss2 as		--diffdiss ist der Ergebnis des Verschneidungsprozesses mmit st_difference		
select
         casttomultipolygon(st_difference(st_union(a.geom),st_union(a.geom))) AS geom --Der Layer diffstanz schneidet Löcher in die Geometrie von diffvorne
      --casttomultipolygon(st_difference(st_union(buffer(a.geom,0.01)),st_union(buffer(b.geom,0.01)))) AS geom 			--Der Layer diffstanz schneidet Löcher in die Geometrie von diffvorne
      from gemeinden as a inner join dlm250 as b on st_intersects(a.geom,b.geom) = 1
       where b.nutzung = 'AX_Wald' and  b.ROWID IN
              (SELECT ROWID FROM SpatialIndex WHERE f_table_name='dlm250' AND search_frame=a.geom);

Select RecoverGeometryColumn('diffdiss2','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('diffdiss2','geom');
SELECT UpdateLayerStatistics('diffdiss2');


-------------
--ATTRIBUTE ÜBER EIN INTERSET AUF DEN DIFFERENZLAYER ÜBERTRAGEN
---------------
create table differenzfertig as
select
  casttomultipolygon(st_intersection(a.geom,diffdiss.geom)) AS geom,   -- Über ein Intersect werden Einzelgeometrien und  Attribute wieder hergestellt
  a.* 
  from gemeinden as a inner join diffdiss on st_intersects(a.geom,diffdiss.geom) = 1  -- HIER MUSS FÜR DIFF_EINGABE DER URSPRUNGSLAYER A EINGETRAGEN WERDEN
   where  a.ROWID IN
              (SELECT ROWID FROM SpatialIndex WHERE f_table_name='gemeinden' AND search_frame=diffdiss.geom);

Select RecoverGeometryColumn('differenzfertig','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('differenzfertig','geom');
SELECT UpdateLayerStatistics('differenzfertig');

-------
--------
---EINE VEREINIGUNG ERREICHT MAN NICHT ÜBER GUNION - DIES FÜHRT ZU DOPPELTEN GEOMETRIEN
-- SONDERN IN DEM ÜBER EINE NORMALE ZEILEN-UNION DIE ERGEBNISSE DES INTERSECT UND DES DIFFERENZ-LAYERS VEREINIGT WERDEN
create table vereinigung as
select
geom,
rs,
gen,
bez,
nutzung
from verschneidung
union
select
geom,
rs,
gen,
bez,
'andere' as nutzung
from differenzfertig;

Select RecoverGeometryColumn('vereinigung','geom',25832,'MULTIPOLYGON');
Select CreateSpatialIndex('vereinigung','geom');
SELECT UpdateLayerStatistics('vereinigung');




/*##############################################
/* Effiziente Geometrieverschneidung im PostGis  */       
       
       /* Differenz mit ST_Differenz dauerte bei beim Testprojekt 3h   */
             
             drop table if exists ausgesstanzt;
                create table ausgesstanzt as
                with
                    stanze as
                     (
                     select
                     ST_Multi(st_union(geom))::geometry('MULTIPOLYGON',25832) as geom
                     from tbl_rechts
                     )

                select
                a.gid as gid,
                --ST_multi(ST_CollectionExtract((ST_difference(st_makevalid(a.geom),st_makevalid(b.geom))),3))::geometry('MULTIPOLYGON',25832) as geom,
                ST_MULTI((ST_Dump(st_difference(a.geom,b.geom))).geom)::geometry('MULTIPOLYGON',25832) as geom

                from  tbl_links as a, stanze as b
                ;
                        ALTER TABLE ausgesstanzt	                -- 
                            ADD COLUMN gid SERIAL PRIMARY KEY; 		-- Primärschlüssel hinzufügen
                    
                    CREATE INDEX idx_ausgesstanzt_gix ON ausgesstanzt USING GIST (geom);  --Räumlichen Index anlegen
                    
                    
--#################################################################################################################################################
/* Differenz- und Union-Alternative über das Zerlegen der Geometrien, dauerte beim testprojekt 10 Minuten */
                   
                    drop table if exists pol_aus_segment;
                    create table pol_aus_segment as      -- Stufe 1: Eine Geometrievereingung beider Layer ohne Attribute 
                        with 
                            segmente as                  -- Erzeugen der Segmente aus den Polygonen mit ST_Boundary (Kanten extrahieren) und ST_Dump (Zerleggen in Einzelsegmente)
                            (
                                select
                                (ST_Dump(ST_Boundary(geom))).geom::geometry('Linestring',25832) as geom  -- Kanten der Polygone des Layers tbl_links extrahieren und in Segmente zerlegen
                                from tbl_links
                                    ----
                                    UNION      --vereinigen  des Segmente des linken und rechen Layers             
                                    -----
                                select
                                (ST_Dump(ST_Boundary(geom))).geom::geometry('Linestring',25832) as geom  -- Kanten der Polygone des Layers tbl_rechst extrahieren und in Segmente zerlegen
                                from tbl_rechts
                            )
                            ,
                            verei as
                            (
                                select 
                                (ST_Dump(ST_UNION(geom))).geom::geometry('Linestring',25832) as geom   -- Zusammenführen und Knoten erzeugen mit ST_Union / In Segmente zerlegen mit ST_Dump
                                from segmente
                            )
                            select                                                                                            --Segmente zu Polygonen mit ST_polygonize                                                               
                            (ST_Dump(ST_CollectionExtract(ST_POLYGONIZE(geom),3))).geom::geometry('POLYGON',25832) as geom   -- Mit ST_collection Nicht-Pols ausschließen
                            from verei
                        ;
                        
                    ALTER TABLE pol_aus_segment	                -- 
                        ADD COLUMN gid SERIAL PRIMARY KEY; 		-- Primärschlüssel hinzufügen
                    
                    CREATE INDEX idx_pol_aus_segment_gix ON pol_aus_segment USING GIST (geom);  --Räumlichen Index anlegen
                    

                    /* Eine neue Geometriespalte anlegen und eine Punktgeometrie auf die Polygonoberfläche für die spätere Attributübertragung anlegen */
                    alter table pol_aus_segment
                        add column poi geometry('POINT',25832)
                    ;
                    update pol_aus_segment
                        set 
                        poi = ST_PointonSurface(geom)
                    ;
                    
                    CREATE INDEX idx_seg_poi_gix ON pol_aus_segment  USING GIST (poi);  -- Räumlicher Index für die Punkte
                    

/* Unions-Layer mit übertragen der Attribute. Aus diesem kann anschließend der Differenz-Layer gebildet werden. */
                    drop table if exists pol_union;
                    create table pol_union as
                        with 
                            cte1 as   -- Anbinden von Attributen des linken Layers über räumliche Anfrage: Hier der Primärschlüssel gid
                                (
                                select
                                a.gid as gid,                           --Primärschlüssel aus der Geometrivereinigung
                                coalesce(b.gid,0) as id_links,                      -- Primärschlüssel aus dem linken Layer, wenn keine Überschneidung dann 0
                                0 as id_rechts,                         -- 0 für den Primärschlüssel von rechts, der hier nicht abgefragt wird
                                a.geom::geometry('POLYGON',25832) as geom

                                from pol_aus_segment as a left join tbl_links as b on st_intersects(a.poi,b.geom) = true  -- Räumliche Beziehung über die Punktgeometrie der Geometrivereinigung
                                )
                                ,
                            cte2 as    --  Anbinden von Attributen des rechten Layers über räumliche Anfrage: Hier der Primärschlüssel gid
                                (
                                select
                                a.gid as gid,             --Primärschlüssel aus der Geometrivereinigung
                                0 as id_links,             -- 0 für den Primärschlüssel von links, der hier nicht abgefragt wird
                                coalesce(b.gid) as id_rechts,           -- Primärschlüssel aus dem rechten Layer, wenn keine Überschneidung dann 0
                                a.geom::geometry('POLYGON',25832) as geom                             
                                                            
                                from pol_aus_segment as a left join tbl_rechts as b on st_intersects(a.poi,b.geom) = true  -- Räumliche Beziehung über die Punktgeometrie der Geometrivereinigung
                                )
                                ,
                            ver as    -- Vereinigung der Attributtierungn von links und rechts
                                (
                                select * from cte1
                                    ----
                                    UNION      --vereinigen  des Segmente des linken und rechen Layers             
                                    -----
                                select * from cte2
                                )
                        --Gruppierung nach gid und geom der Gemetrievereinigung, um doppelte Objekte nach dem Anbinden der Attribute von beiden Seiten zu vermeinden
                        select
                            gid,
                            MAX(id_links) as id_links,    --Primärschlüssel erzeugen-- Es muss eine Aggregierungsfunktion verwendet werden anders als bei SpatiaLite, darum Max
                            MAX(id_rechts) as id_rechts, -- Es muss eine Aggregierungsfunktion verwendet werden anders als bei SpatiaLite, darum Max
                            ST_MULTI(ST_union(geom))::geometry('MULTIPOLYGON',25832) as geom 
                        from ver
                        where id_links > 0 or id_rechts > 0  -- Wenn beide 0 sind, dann handelt es sich um Polygone, die durch zufällige geschlossene Linien entstanden sind
                        group by gid,geom
                        ;
       /*   Layer mit der Geometrievereinigung löschen */
                  drop table pol_aus_segment;
					
    /*   Primärschlüssel etv */
                        ALTER TABLE    pol_union
                            ADD CONSTRAINT pol_union_prim PRIMARY KEY(gid); --Primärschlüssel erzeugen
                        CREATE INDEX idx_pol_union_gix ON pol_union USING GIST (geom);	
		                    
/* Differenz-Layer erzeugen */
                        drop table if exists pol_diff;
                        Create table pol_diff
                        as
                            select
                            gid,
                            id_links,
                            geom::geometry('MULTIPOLYGON',25832) as geom 
                            from pol_union
                            where id_links > 0 and id_rechts = 0
                        ;
                        
                        ALTER TABLE    pol_diff
                            ADD CONSTRAINT pol_diff_prim PRIMARY KEY(gid); --Primärschlüssel erzeugen
                        
                        CREATE INDEX idx_pol_diff_gix ON pol_diff USING GIST (geom);	
                        
 /* Intersect-Layer erzeugen */                       
                        drop table if exists pol_inter;
                        Create table pol_inter
                        as
                            select
                            gid,
                            id_links,
                            id_rechts,
                            geom::geometry('MULTIPOLYGON',25832) as geom 
                            from pol_union
                            where id_links > 0 and id_rechts > 0
                        ;
                        
                        ALTER TABLE    pol_inter
                            ADD CONSTRAINT pol_inter_prim PRIMARY KEY(gid); --Primärschlüssel erzeugen
                            
                        CREATE INDEX idx_pol_inter_gix ON pol_inter USING GIST (geom);	

/* SymDifferenz-Layer erzeugen */
                        drop table if exists pol_symdiff;
                        Create table pol_symdiff
                        as
                            select
                            gid,
                            id_links,
                            id_rechts,
                            geom::geometry('MULTIPOLYGON',25832) as geom 
                            from pol_union
                            where (id_links = 0 and id_rechts > 0) or (id_links > 0 and id_rechts = 0 )
                        ;
                        
                        ALTER TABLE    pol_symdiff
                            ADD CONSTRAINT pol_symdiff_prim PRIMARY KEY(gid); --Primärschlüssel erzeugen
                            
                        CREATE INDEX idx_pol_symdiff_gix ON pol_symdiff USING GIST (geom);
                        
                 


/* 
-----------------------------------------
GKG-Kassel - Dr.-Ing. Claas Leiner
QGIS-Support und mehr

Geodatenservice, Kartenwerkstatt &
GIS-Schule Kassel

Wilhelmshöher Allee 304 E
34131 Kassel
Tel. 0561/56013445
claas.leiner@gkg-kassel.de
http://www.gkg-kassel.de
----------------------------------------
Unterstützen Sie QGIS
QGIS-DE e.V. | http://qgis.de
*/
