OpenStreetMap / GDAL

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

Tato stránka shrnuje poznámky pro práci s daty OpenStreetMap pomocí knihovny GDAL.

Stažení dat

Data pro ČR stáhneme ve formátu OSM/PBF:

wget http://download.geofabrik.de/europe/czech-republic-latest.osm.pbf

Pro kontrolu vypíšeme seznam vrstev:

ogrinfo czech-republic-latest.osm.pbf
Had to open data source read-only.
INFO: Open of `czech-republic-latest.osm.pbf'
      using driver `OSM' successful.
1: points (Point)
2: lines (Line String)
3: multilinestrings (Multi Line String)
4: multipolygons (Multi Polygon)
5: other_relations (Geometry Collection)

Příklad výběru dálnic (viz [1]) pomocí konzolových nástrojů knihovny GDAL:

ogrinfo czech-republic-latest.osm.pbf lines -so
ogrinfo czech-republic-latest.osm.pbf lines -where "highway = 'motorway'"

Převod do Esri Shapefile:

ogr2ogr -f 'Esri Shapefile' dalnice.shp czech-republic-latest.osm.pbf lines -where "highway = 'motorway'"

Ukázka skriptu v programovacím jazyku Python

Níže uvedený skript vytvoří na základě dat OSM vektorovou vrstvu ve formátu Esri Shapefile obsahující silniční síť (silnice.shp) v souřadnicovém systému S-JTSK (EPSG 5514).

V atributové tabulce je uveden typ silnice (dálnice, rychlostní silnice, silnice 1.třídy, ...).

  1 #!/usr/bin/env python
  2 # -*- coding: utf-8 -*-
  3 
  4 import os
  5 import sys
  6 
  7 from osgeo import ogr, osr
  8 
  9 osm_file = "czech-republic-latest.osm.pbf"
 10 osm_layer = "lines"
 11 shp_file = "silnice.shp"
 12 
 13 # otevrit vstupni OSM soubor
 14 idrv = ogr.GetDriverByName("OSM")
 15 if idrv is None:
 16     sys.exit("Nelze cist OSM data")
 17 ids = idrv.Open(osm_file, False)
 18 if ids is None:
 19     sys.exit("Nelze otevrit OSM soubor '%'" % osm_file)
 20 ilayer = ids.GetLayerByName(osm_layer)
 21 if ilayer is None:
 22     sys.exit("Nevalidni OSM soubor")
 23 os.environ['OGR_INTERLEAVED_READING'] = 'YES' # nutne pro cteni vetsiho objemu dat
 24 
 25 # vytvorit prazdny shapefile
 26 odrv = ogr.GetDriverByName("ESRI Shapefile")
 27 if odrv is None:
 28     sys.exit("Nelze zapisovat do formatu Esri Shapefile")
 29 
 30 # prepsat existujici shapefile
 31 if os.path.exists(shp_file):
 32     odrv.DeleteDataSource(shp_file)
 33 ods = odrv.CreateDataSource(shp_file)
 34 if ods is None:
 35     sys.exit("Nelze vytvorit Shapefile soubor '%s'" % shp_file)
 36 
 37 # nastavit souradnicove system pro vstup a vystup (S-JTSK)
 38 srs_source = ilayer.GetSpatialRef()
 39 srs_target = osr.SpatialReference()
 40 srs_target.ImportFromEPSG(5514)
 41 
 42 # nastavit transformaci
 43 coord_trans = osr.CoordinateTransformation(srs_source, srs_target)
 44 olayer = ods.CreateLayer(shp_file.split('.')[0],
 45                          srs_target, ogr.wkbLineString, ['ENCODING=UTF-8'])
 46 
 47 # vytvorit atributy pro vystupni vrstvu
 48 field_name = ogr.FieldDefn("typ", ogr.OFTString)
 49 field_name.SetWidth(254)
 50 olayer.CreateField(field_name)
 51 
 52 featDefn = olayer.GetLayerDefn()
 53 
 54 # slovnik pro konverzi typu komunikace
 55 stype = {'motorway'  : 'dálnice',
 56          'trunk'     : 'rychlostní silnice',
 57          'primary'   : '1. třída',
 58          'secondary' : '2. třída',
 59          'tertiary'  : '3. třída'}
 60 
 61 read_osm = True # opetovne cteni je ryze specialita OSM driveru !!!
 62 while read_osm:
 63     read_osm = False
 64     
 65     for idx in range(ids.GetLayerCount()):
 66         layer = ids.GetLayer(idx)
 67         layer_name = layer.GetName()
 68         print("Reading '%s'..." % layer_name)
 69         while True:
 70             # nacteni dalsiho prvku ze vstupu
 71             ifeature = layer.GetNextFeature()
 72             if ifeature is None:
 73                 break
 74             
 75             read_osm = True
 76             if layer_name != 'lines': # preskocit ostatni vrstvy
 77                 continue
 78 
 79             # manualni atributovy filtr
 80             ltype = ifeature.GetField('highway')
 81             if not ltype or ltype not in ('motorway', 'trunk', 'primary', 'secondary', 'tertiary',
 82                                           'motorway_link', 'trunk_link', 'primary_link',
 83                                           'secondary_link', 'tertiary_link'):
 84                 continue
 85             
 86             # transformace geometrie do S-JTSK
 87             geom = ifeature.GetGeometryRef()
 88             geom.Transform(coord_trans)
 89             
 90             # zapis noveho prvku do vystupni vrstvy
 91             ofeature = ogr.Feature(featDefn)
 92             ofeature.SetGeometry(geom)
 93             ofeature.SetField("typ", stype.get(ifeature.GetField('highway').split('_',1)[0], 'neklasifikováno'))
 94             olayer.CreateFeature(ofeature)
 95         
 96 print("%d fearures written to output" % olayer.GetFeatureCount())
 97 
 98 # zavreni vstupni a vystupni vrstvy
 99 ids.Destroy()
100 ods.Destroy()

Kontrola výsledku:

ogrinfo silnice.shp silnice -so
INFO: Open of `silnice.shp'
      using driver `ESRI Shapefile' successful.

Layer name: silnice
Geometry: Line String
Feature Count: 98745
Extent: (-905139.251535, -1225746.274104) - (-430013.772685, -934787.131005)
Layer SRS WKT:
PROJCS["S_JTSK_Krovak_East_North",
    GEOGCS["GCS_S-JTSK",
        DATUM["System_Jednotne_Trigonometricke_Site_Katastralni",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.83333333333333],
    PARAMETER["azimuth",30.28813972222222],
    PARAMETER["pseudo_standard_parallel_1",78.5],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
typ: String (254.0)
Vizualizace silniční sítě v QGISu