Cleaned up province segmentation

stimmen/cbs.py

@@ -0,0 +1,93 @@
+from threading import Thread
+import time
+
+from pygeoif.geometry import mapping
+from shapely.geometry import shape, box
+
+from stimmen.shapefile import shapefiles_to_geojson
+
+import os.path
+
+
+def data_file(*args):
+    return os.path.join(
+        os.path.dirname(os.path.abspath(__file__)),
+        '..',
+        *args
+    )
+
+__cache = {}
+__last_access = None
+
+
+def clear_cache_poll():
+    global __cache, __last_access
+    while True:
+        time.sleep(60*60)
+        if (time.time() - __last_access) > 60*60:
+            __cache = {}
+
+
+Thread(target=clear_cache_poll, daemon=True).start()
+
+
+def get_available_provinces():
+    file = data_file('data', 'geo', 'cbs', "2018_Imergis_provinciegrenzen_met_water")
+    provinces = shapefiles_to_geojson(file)['features']
+    return {province['properties']['provincien'] for province in provinces}
+
+
+def province_geojson(province, with_water=False):
+    global __cache, __last_access
+    file = data_file(
+        'data', 'geo', 'cbs',
+        "2018_Imergis_provinciegrenzen_met_water"
+        if with_water else
+        "2018-Imergis_provinciegrenzen_kustlijn"
+    )
+
+    if (province, with_water) not in __cache:
+        provinces = shapefiles_to_geojson(file)['features']
+        available_provinces = get_available_provinces()
+
+        assert province in available_provinces, "Province {} not found, options: {}".format(
+            province, ', '.join(available_provinces))
+
+        __cache[(province, with_water)] = next(
+            province_ for province_ in provinces if province_['properties']['provincien'] == province)
+
+    __last_access = time.time()
+    return __cache[(province, with_water)]
+
+
+def gwb_in_province(
+        province='Friesland', region_level='wijk', region_year='2018',
+        polygon_simplification=0.001, province_dilation=0.0005
+):
+    assert region_level in {'gem', 'wijk', 'buurt'}, (
+        "region_level {} not supported, must be gem, wijk or buurt".format(region_level))
+    assert region_year in {2017, 2018}, (
+        "region_year {} not supported, must 2017 or 2018".format(region_year))
+
+    province_with_water = shape(province_geojson(province, with_water=True)['geometry'])
+    province_bounding_box = box(*province_with_water.bounds)
+
+    province_land_only = shape(province_geojson(province, with_water=False)['geometry'])
+    province_land_only_dilated = province_land_only.buffer(-province_dilation)
+
+    geojson = shapefiles_to_geojson(
+        data_file('data', 'geo', 'cbs', '{}_{}'.format(region_level, region_year)))
+    shapes = [shape(geojson_['geometry']) for geojson_ in geojson['features']]
+
+    shapes, geojson = map(list, zip(*(
+        (intersection.simplify(tolerance=polygon_simplification), geojson_)
+        for shape_, geojson_ in zip(shapes, geojson['features'])
+        if province_bounding_box.contains(shape_)
+        for intersection in [shape_.intersection(province_land_only_dilated)] # alias
+        if intersection.area > 0
+    )))
+
+    for geojson_, shape_ in zip(geojson, shapes):
+        geojson_['geometry'] = mapping(shape_)
+
+    return {"type": "FeatureCollection", "features": geojson}

stimmen/geojson.py

@@ -0,0 +1,43 @@
+from pygeoif.geometry import mapping
+from shapely.geometry import shape
+
+
+def merge_features(geojson, condition, aggregate={}):
+    """Merge the geometries using shapely's union for all the geojson's features that
+    meet the condition, condition get's passed an item of feature. Then aggregate the properties
+    in the aggregate dict using a function that get as input the list of property values off alle
+    matched features. Operates inplace."""
+    indices = [index for index, feature in enumerate(geojson['features']) if condition(feature)]
+    if len(indices) == 0: # also if there is one index, we
+        return geojson
+    properties = {
+        prop: agg([
+            properties_[prop]
+            for index in indices
+            for properties_ in [geojson['features'][index]['properties']]
+            if prop in properties_
+        ])
+        for prop, agg in aggregate.items()
+    }
+    properties.update({
+        key: value
+        for index in indices
+        for key, value in geojson['features'][index]['properties'].items()
+        if key not in aggregate
+    })
+    if len(indices) == 1:
+        geojson['features'][indices[0]]['properties'] = properties
+        return geojson
+
+    union = shape(geojson['features'][indices[0]]['geometry'])
+    for index in indices[1:]:
+        union = union.union(shape(geojson['features'][index]['geometry']))
+
+    for index in indices[::-1]: # reverse, such that the 'todo' indices willnot change.
+        del geojson['features'][index]
+
+    geojson['features'].append({
+        'geometry': mapping(union),
+        'properties': properties
+    })
+    return geojson

stimmen/latitude_longitude.py

@@ -0,0 +1,9 @@
+
+
+def reverse_latitude_longitude(rd_multipolygon):
+    if len(rd_multipolygon) == 2 and tuple(map(type, rd_multipolygon)) == (float, float):
+        return rd_multipolygon[::-1]
+    return [
+        reverse_latitude_longitude(element)
+        for element in rd_multipolygon
+    ]

stimmen/shapefile.py

@@ -0,0 +1,43 @@
+from osgeo.osr import SpatialReference, CoordinateTransformation
+import shapefile
+
+epsg28992 = SpatialReference()
+epsg28992.ImportFromEPSG(28992)
+
+epsg28992.SetTOWGS84(565.237 ,50.0087 ,465.658 ,-0.406857 ,0.350733 ,-1.87035 ,4.0812)
+
+epsg4326 = SpatialReference()
+epsg4326.ImportFromEPSG(4326)
+
+rd2latlon = CoordinateTransformation(epsg28992, epsg4326)
+
+
+def rd_to_longlat(rd_multipolygon):
+    if len(rd_multipolygon) == 2 and tuple(map(type, rd_multipolygon)) == (float, float):
+        return list(rd2latlon.TransformPoint(*rd_multipolygon)[:2])
+    return [
+        rd_to_longlat(element)
+        for element in rd_multipolygon
+    ]
+
+
+def shapefiles_to_geojson(filename):
+    shape_file = shapefile.Reader(filename)
+    fields = shape_file.fields[1:]
+    field_names = [field[0] for field in fields]
+
+    buffer = []
+    for shape_record in shape_file.shapeRecords():
+        properties = dict(zip(field_names, shape_record.record))
+        geometry = shape_record.shape.__geo_interface__
+        geometry['coordinates'] = rd_to_longlat(geometry['coordinates'])
+
+        buffer.append({
+            'type': "Feature",
+            'geometry': geometry,
+            'properties': properties
+        })
+    return {
+        "type": "FeatureCollection",
+        "features": buffer
+    }