RUIAN / GDAL / Python

Z FreeGIS portál
Přejít na: navigace, hledání

Tato stránka shromažďuje příklady uživatelských skriptů v jazyku Python pro práci s daty ve výměnném formátu RÚIAN (VFR).

Export vybrané tématické VFR vrstvy do formátu Esri Shapefile

Vizualizace vrstvy 'ulice' z VFR ve formátu Esri Shapefile s podkladovou vrstvou ortofota (ČÚZK)
  1 #!/usr/bin/env python
  2 
  3 import os
  4 import sys
  5 import datetime
  6 
  7 try:
  8     from osgeo import ogr
  9 except:
 10     sys.exit('Knihovna GDAL neni dostupna')
 11 
 12 # otevrit vstupni VFR soubor
 13 def open_vfr(input_file, layer_name):
 14     driver = ogr.GetDriverByName('GML')
 15     if driver is None:
 16         sys.exit("Nelze cist VFR data, chybi OGR GML driver")
 17 
 18     datasource = driver.Open(input_file, False)
 19     if datasource is None:
 20         sys.exit("Nelze otevrit VFR soubor")
 21 
 22     layer = datasource.GetLayerByName(layer_name)
 23     if layer is None:
 24         sys.exit("Vrstva '%s' nenalezena" % layer_name)
 25 
 26     return datasource, layer
 27 
 28 # vytvorit prazdny Shapefile soubor
 29 def create_shp(layer_vfr, geom_idx = 0):
 30     driver = ogr.GetDriverByName('Esri Shapefile')
 31     if driver is None:
 32         sys.exit("Nelze zapisovat data do formatu Shapefile")
 33 
 34     shp_name = layer_vfr.GetName() + '.shp'
 35     if os.path.exists(shp_name):
 36         driver.DeleteDataSource(shp_name)
 37 
 38     datasource = driver.CreateDataSource(shp_name)
 39     if datasource is None:
 40         sys.exit("Nelze vytvorit vystupni soubor '%s'" % shp_name)
 41     
 42     defn_vfr = layer_vfr.GetLayerDefn()
 43     defn_geom = defn_vfr.GetGeomFieldDefn(geom_idx)
 44     if defn_geom is None:
 45         sys.exit("Atribut geometrie s index %d neexistuje" % geom_idx)
 46         
 47     layer_shp = datasource.CreateLayer(layer_vfr.GetName(),
 48                                        layer_vfr.GetSpatialRef(),
 49                                        defn_geom.GetType(),
 50                                        ["ENCODING=UTF-8"])
 51 
 52     return datasource, layer_shp
 53 
 54 # prekopirovat prvky z vybrane VFR vrstvy do ciloveho Shapefile
 55 # souboru (zachovat pouze vybrany atribut s geometrii)
 56 def copy_data(layer_vfr, layer_shp, geom_idx = 0):
 57 
 58     # vytvorit atributy pro vystupni vrstvu
 59     layer_defn = layer_vfr.GetLayerDefn()
 60     for i in range(0, layer_defn.GetFieldCount()):
 61         layer_shp.CreateField(layer_defn.GetFieldDefn(i))
 62     
 63     # pridat atribut geometrie (GDAL >= 2.0)
 64     ### layer_shp.CreateGeomField(layer_defn.GetGeomFieldDefn(geom_idx))
 65     
 66     # cist VFR prvky
 67     layer_vfr.ResetReading()
 68     while True:
 69         feature = layer_vfr.GetNextFeature()
 70         if feature is None:
 71             break
 72 
 73         # vytvorit kopii prvku ze vstupni vrstvy
 74         ofeature = ogr.Feature(layer_shp.GetLayerDefn())
 75         ofeature.SetFrom(feature, True)
 76         
 77         # ziskat geometrii z pozadovaneho atributu
 78         geom = feature.GetGeomFieldRef(geom_idx)
 79         if geom is None:
 80             continue # preskocit prvky bez geometrie
 81         
 82         # geometrii nastavit pro danou kopii prvku
 83         ofeature.SetGeometry(geom)
 84         
 85         # novy prvek pridat to vystupni vrstvy
 86         layer_shp.CreateFeature(ofeature)
 87 
 88 def datum():
 89     today = datetime.date.today()
 90     if today.month == 12:
 91         day = today.replace(day=31)
 92 
 93     day = (today.replace(month=today.month, day=1) - datetime.timedelta(days=1))
 94 
 95     return day.strftime("%Y%m%d")
 96        
 97 def main(vfr_file, layer_name, geom_idx = 0):
 98     ids, layer_vfr = open_vfr(vfr_file, layer_name)
 99     ods, layer_shp = create_shp(layer_vfr, geom_idx)
100     copy_data(layer_vfr, layer_shp, geom_idx)
101     
102     ids.Destroy()
103     ods.Destroy()
104 
105 if __name__ == "__main__":
106     main("/vsicurl/http://vdp.cuzk.cz/vymenny_format/soucasna/{}_OB_554821_UKSH.xml.gz".format(datum()),
107          "Ulice")

Příklad výstupu do geodatabáze PostGIS

 1 # vytvorit novou vrstvu (tj. feature table) v databazi PostGIS
 2 def create_pg(layer_vfr, geom_idx = 0):
 3     driver = ogr.GetDriverByName('PostgreSQL')
 4     if driver is None:
 5         sys.exit("Nelze zapisovat data do formatu Shapefile")
 6 
 7     conninfo = "PG:dbname=pgis_student user=postgis password=***** host=geo102.fsv.cvut.cz"
 8     datasource = driver.Open(conninfo)
 9     if datasource is None:
10         sys.exit("Nelze se pripojit k DB")
11     
12     defn_vfr = layer_vfr.GetLayerDefn()
13     defn_geom = defn_vfr.GetGeomFieldDefn(geom_idx)
14     if defn_geom is None:
15         sys.exit("Atribut geometrie s index %d neexistuje" % geom_idx)
16         
17     table_name = layer_vfr.GetName().lower() # velka pismena v PG delaji neplechu
18     layer_pg = datasource.CreateLayer(table_name,
19                                       layer_vfr.GetSpatialRef(),
20                                       defn_geom.GetType())
21     if layer_pg is None:
22         sys.exit("Nelze vytvorit tabulku '%s'" % table_name)
23     
24     return datasource, layer_pg

Užitečné odkazy