Přeskočit obsah

10. PostGIS - import dat a jejich zpracování

Import a export dat

SQL

Import:

CREATE TABLE poi (fid serial primary key,
                  nazev varchar(100),
                  geom geometry(point, 5514));

INSERT INTO poi(nazev, geom) VALUES('Litoměřice',
   ST_GeomFromText('POINT(-756371 -991031)', 5514));

Export:

COPY (select fid,nazev,st_astext(geom) from poi) TO /tmp/poi.csv DELIMITER ',' CSV HEADER;

-- vs.

\copy (select fid,nazev,st_astext(geom) from poi) TO '/tmp/poi.csv' DELIMITER ',' CSV HEADER;

Nástroje PostGISu

(pouze Esri Shapefile)

Import:

shp2pgsql -s 5514 -I ChraneneUzemi.shp chranene_uzemi > chranene_uzemi.sql
psql uzpr -h gislab -U <user> -f chranene_uzemi.sql

# vs.

shp2pgsql -s 5514 -d -I ChraneneUzemi.shp chranene_uzemi | psql uzpr -h gislab -U <user> 

Export:

pgsql2shp -h gislab -u <user> -P <password> uzpr osm.zeleznice

Nástroje knihovny GDAL

(pod MS Windows dostupné z OSGeo4W Shell)

Import:

ogr2ogr -f PostgreSQL PG:"dbname=uzpr host=gislab user=<user> password=<password>" arccr500.sqlite kladyzakladnichmap
# vs.

ogr2ogr -f PGDump lesy.sql arccr500.sqlite lesy
psql uzpr -h gislab -U <user> -f lesy.sql

Export:

ogr2ogr -f GPKG zeleznice.gpkg PG:"dbname=uzpr host=gislab user=<user> password=<password>" osm.zeleznice

Nástroje QGIS

Ve správci databáze (DB Manager) použijte nástroj Import Layer/File a Export to File.

Cvičení

Pokračování úloh z GIS1

Opakování úloh z předmětu GIS1 2. cvičení a 3. cvičení. Dokončení úloh z osmé lekce.

Otázky jsou postaveny nad datasetem ArcČR v3.3 (S:\K155\Public\data\GIS\ArcCR500 3.3).

Praktická úloha

Porovnejme shodu třídy půdního pokryvu datasetu LUCAS na první úrovni s vrstvou lesa (zdroj: Data50).

Nejprve data naimportujeme (dejte si pozor na velká písmena v názvech tabulek). V případě datasetu LUCAS můžeme využít specializovaný QGIS plugin.

Import z příkazové řádky

    ogr2ogr -f PostgreSQL PG:"dbname=uzpr host=gislab user=xxx password=xxx" lucas2018_cz.gpkg lucas_points
    ogr2ogr -f PostgreSQL -lco precision=NO PG:"dbname=uzpr host=gislab user=xxx password=xxx" data50/Les.shp

Nejprve zjistíme počet bodů LUCAS s půdním pokryvem "Woodland":

SELECT count(*) FROM lucas2018_cz;
SELECT count(*) FROM lucas2018_cz WHERE lc1 LIKE 'C%';

Před operací překrytí s vrstvou lesů transformujeme LUCAS body do stejného souřadnicového systému:

SELECT * FROM geometry_columns WHERE f_table_name = 'les';
ALTER TABLE lucas2018_cz ADD COLUMN geom_5514 geometry(point, 5514);
UPDATE lucas2018_cz SET geom_5514 = st_transform(geom, 5514);
CREATE INDEX ON lucas2018_cz USING GIST (geom_5514);

Zjistíme počet bodů se shodou půdního pokryvu "Woodland" (LUCAS) a vrstvy lesů (Data50):

SELECT count(*) FROM lucas2018_cz AS p
JOIN les AS l
ON st_within(p.geom_5514, l.geom)
WHERE p.lc1 LIKE 'C%';

Podívejme se detailně na body, u kterých ke shodě nedošlo:

SELECT point_id,p.geom_5514 FROM lucas2018_cz AS p
LEFT JOIN les AS l
ON st_within(p.geom_5514, l.geom)
WHERE p.lc1 LIKE 'C%' AND l.geom IS NULL;

Pro takto identifikované body spočítejme vzdálenost k nejkratšímu lesu:

WITH lucas_nomatch AS (
 SELECT point_id,p.geom_5514 FROM lucas2018_cz AS p
 LEFT join les AS l
 ON st_within(p.geom_5514, l.geom)
 WHERE p.lc1 LIKE 'C%' AND l.geom IS NULL
)
SELECT point_id,min(p.geom_5514 <-> l.geom) AS dist_m
FROM lucas_nomatch AS p, les AS l
GROUP BY point_id order BY dist_m DESC;

Úkol

Dotaz lze optimalizovat použitím cross join lateral, viz zde.

WITH lucas_nomatch AS (
 SELECT point_id,p.geom_5514 FROM lucas2018_cz AS p
 LEFT JOIN les AS l
 ON st_within(p.geom_5514, l.geom)
 WHERE p.lc1 LIKE 'C%' AND l.geom IS NULL
)
SELECT p.point_id,dist_m FROM lucas_nomatch AS p
CROSS JOIN LATERAL (
 SELECT p.point_id, p.geom_5514 <-> l.geom AS dist_m
 FROM les AS l
 ORDER BY dist_m LIMIT 1
) AS l ORDER BY dist_m DESC;

Na základě této analýzy můžeme zahrnout i body, které leží do 50m od nejbližšího lesa:

SELECT count(*) FROM lucas2018_cz AS p
JOIN les AS l
ON st_dwithin(p.geom_5514, l.geom, 50)
WHERE p.lc1 LIKE 'C%';

Úkol

Proveďte podobnou analýzu pro vodní plochy:

SELECT count(*) FROM lucas2018_cz WHERE lc1 LIKE 'G%';
...