2

I am trying to read in openstreetmap data (osm.pbf) files into python. Initially, I've always used the ogr2ogr command line approach to convert the data I want from the osm.pbf to a shapefile and then read the shapefile using geopandas. However, I don't really need shapefiles so I thought it would be nice to load the data directly into python.

Unfortunately, I seem to miss ways/segments when directly reading it into python. Please see below the code I am using for both approaches:

# path to osm.pbf file
osm_path = 'C:\\planet_osm\\tanzania-latest.osm.pbf'

# import through python
driver=ogr.GetDriverByName('OSM')
data = driver.Open(osm_path)
layer = data.GetLayerByName('lines')

features=[x for x in layer]
railway = []
for feature in features:
    if feature.GetField('railway') is not None:
        osm_id = feature.GetField('osm_id')
        shapely_geo = shapely.wkt.loads(feature.geometry().ExportToWkt()) 
        service = feature.GetField('service')
        maxspeed = feature.GetField('maxspeed')
        highway=feature.GetField('railway')
        railway.append([osm_id,highway,service,maxspeed,shapely_geo])

rail_dir_import = gpd.GeoDataFrame(railway,columns=['osm_id','infra_type','service','maxspeed','geometry'],crs={'init': 'epsg:4326'})  

# import through ogr2ogr      
shape_in = 'C:\\planet_osm\\tanzania-railway.shp'

os.system("ogr2ogr -progress -overwrite -f \"ESRI Shapefile\" -sql \
               \"SELECT osm_id,railway FROM lines WHERE railway IS NOT NULL\" \
               -lco ENCODING=UTF-8 -skipfailures -nlt geometry "+shape_in+" "+osm_path)       

rail_com_line = gpd.read_file(shape_in)

The python approach only gives me 563 segments, whereas the ogr2ogr command line approach (and then loading it into python) gives me 1060 segments. See below the output in visual form.

Python approach: https://www.dropbox.com/s/kuxcaiiigdn1u58/python_import.png?dl=0

ogr2ogr command line approach: https://www.dropbox.com/s/a854i0345e62fc4/ogr2ogr_import.png?dl=0

Does anyone know why the python approach misses data?

Update

If I use the following, I do get the right amount of railway segments:

    layer = data.ExecuteSQL("SELECT osm_id,service,maxspeed,railway FROM lines WHERE railway IS NOT NULL")

Still, that does not solve the problem above. It makes me think that there is some sort of caching problem when loading all features from a layer?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.