bezirksteil-plz.py

#!/usr/local/bin/python3

from urllib.parse import quote_plus
from urllib.request import URLopener
import os, subprocess
import json
from shapely.geometry import shape
import csv  


def download_data():
    bbox = '48.07303233901773,11.348190307617188,48.25028349849019,11.73614501953125'
    query = 'rel(' +bbox + ')[boundary=administrative][admin_level=10]; out geom;' + \
            'rel(' +bbox + ')[boundary=postal_code][postal_code~"8[01].+"]; out geom;'


    file = URLopener()
    file.retrieve('http://overpass-api.de/api/interpreter?data=' + quote_plus(query), 'bezirksteile-plz.osm')
    pass

def convert_to_geojson():
    # install via `npm install -g osmtogeojson`
    
    
    with open('bezirksteile-plz.geojson', 'w') as out:
        
        os.environ['PATH'] += ':/usr/local/bin/'
        return_code = subprocess.call(['osmtogeojson', 'bezirksteile-plz.osm'], stdout=out)
    return return_code 


bezirksteile = []
def get_bezirksteile():
    global bezirksteile
    
    with open('bezirksteile-plz.geojson', encoding='utf-8') as data_file:    
        data = json.load(data_file)
        for item in data['features']:
            if 'boundary' in item['properties']:
                if item['properties']['boundary']  == 'administrative':
                    bezirksteile.append({
                        'name': item['properties']['name'].replace('Bezirksteil ', ''),
                        'plz' : [],
                        'shape': shape(item['geometry'])
                    });
                elif item['properties']['boundary']  == 'postal_code':
                    plz_shape = shape(item['geometry'])
                    for bezirksteil in bezirksteile:
                        if plz_shape.intersects(bezirksteil['shape']):
                            
                            #calculate how many percent of this bezirksteil is covered by this postal code area
                            percent = (plz_shape.intersection(bezirksteil['shape']).area)/bezirksteil['shape'].area
                            
                            bezirksteil['plz'].append([item['properties']['postal_code'], percent])
                    
def write_to_csv():                 
    with open('plz-bezirksteile.csv', 'wt', encoding='utf-8') as csvfile:
        w = csv.writer(csvfile) 
        w.writerow(['plz', 'abdeckung', 'bezirksteil']) 
        for bezirksteil in bezirksteile:
            for [plz, percent] in bezirksteil['plz']:
                #print(plz, "{0:.6f}".format(percent), bezirksteil['name'])
                w.writerow([plz, "{0:.6f}".format(percent), bezirksteil['name']])
                
    pass

def main():

    download_data()
    convert_to_geojson()
    get_bezirksteile()
    write_to_csv()
    
    print('done')
    pass

if __name__ == '__main__':
    main()
    pass

plz-bezirksteile.csv

TODO: replace plz_shape.intersects() with within() or something simular – compare http://toblerity.org/shapely/manual.html#object.intersects -> should remove 0.000 entries