• Convert HydroBASIN polygons from SHP to KML
  • Fix geometry errors
  • Simplify geometry (optionally)

In [46]:
%matplotlib inline

import logging
import sys
import os
import glob

import math
import ogr
import shapely.geometry, shapely.wkt
import pylab
import matplotlib.pyplot as plt
from utils.shapely_plot import draw
import shapely as sl
import fiona
import numpy as np

pylab.rcParams['figure.figsize'] = (17.0, 15.0)
logging.basicConfig(stream=sys.stderr, level=logging.INFO)

In [329]:
def convert(input_path):
    output_path = os.path.splitext(input_path)[0] + '.kml'

    filename = os.path.splitext(os.path.basename(input_path))[0]

    dst_ds = ogr.GetDriverByName('KML').CreateDataSource(output_path)
    dst_lyr = dst_ds.CreateLayer(filename)

    # create fields using OGR
    src_ds = ogr.Open(input_path)
    src_lyr = src_ds.GetLayerByIndex(0)
    f = src_lyr.GetFeature(0)
    [dst_lyr.CreateField(f.GetFieldDefnRef(i)) for i in range(f.GetFieldCount())]

    for feat in src_lyr:
        try:
            geom = shapely.wkt.loads(feat.GetGeometryRef().ExportToWkt())
        except Exception as e: 
            print('Error({0}), skipping geometry.'.format(e))
            continue

        #id = feat.GetField('HYBAS_ID')
        count = get_geom_coord_count(geom)

        # geom coord count
        # print('{0}: {1}'.format(id, count))

        if count > 10000:
            print ('simplifying ...')
            geom = geom.simplify(0.004)
        
        # uncomment in case of messed-up geometry
        #
        if not geom.is_valid:
            geom = geom.buffer(0.0)

        # geom = geom.simplify(0.004)

        # if id == 6030021870:
        #    print ('simplifying ...')
        #    geom = geom.simplify(0.004)

        f = ogr.Feature(dst_lyr.GetLayerDefn())

        # set field values
        for i in range(feat.GetFieldCount()):
            fd = feat.GetFieldDefnRef(i)
            f.SetField(fd.GetName(), feat.GetField(fd.GetName()))

        # set geometry    
        f.SetGeometry(ogr.CreateGeometryFromWkt(geom.to_wkt()))
        
        # create feature
        dst_lyr.CreateFeature(f)
        
        f.Destroy() 
        
    src_ds.Destroy()

    dst_ds.Destroy()

# files = glob.glob(r'D:\gis\HydroBASINS\without_lakes\hybas_sa_lev*_v1c.shp')
# files = glob.glob(r'D:\gis\HydroBASINS\without_lakes\simplified\hybas_au_lev08_v1c.shp')

# files = glob.glob(r'D:\gis\GRanD\GRanD_dams_v1_1.shp')
files = glob.glob(r'D:\gis\GRanD\GRanD_reservoirs_v1_1.shp')

# files = glob.glob(r'D:\gis\HydroBASINS\rivers\*_riv_15s.shp')
# files = [r'D:\gis\OpenStreetMaps\Asia\Australia\rivers_lines.shp',
#          r'D:\gis\OpenStreetMaps\Asia\Australia\rivers_multipolygons.shp']

# files = [r'D:\gis\OpenStreetMaps\simplified-land-polygons-complete-3857\simplified_land_polygons.shp']
    
# convert(r'D:\src\GitHub\data-utils\notebooks\shapefile.shp')

# convert(r'D:\gis\HydroBASINS\without_lakes\hybas_au_lev07_v1c.shp')

for f in files:
    convert(f)


WARNING:shapely.geos:Ring Self-intersection at or near point -121.14047832043306 41.899166666666702
WARNING:shapely.geos:Ring Self-intersection at or near point -120.37725641358978 37.795833333333881
WARNING:shapely.geos:Ring Self-intersection at or near point -121.9911453024788 37.191388888889392
WARNING:shapely.geos:Ring Self-intersection at or near point -109.62555555555888 40.95676045735663
WARNING:shapely.geos:Ring Self-intersection at or near point -118.33805555555764 40.628888888889392
WARNING:shapely.geos:Ring Self-intersection at or near point -116.5180555555586 32.692777777778062
WARNING:shapely.geos:Ring Self-intersection at or near point -98.637777777778751 57.488611111111545
WARNING:shapely.geos:Ring Self-intersection at or near point -95.143611111111809 56.419444444444792
WARNING:shapely.geos:Ring Self-intersection at or near point -102.93277777777537 57.360833333333751
WARNING:shapely.geos:Ring Self-intersection at or near point -93.139128892686728 50.573333333333835
WARNING:shapely.geos:Ring Self-intersection at or near point -95.102777777772943 50.218200389144911
WARNING:shapely.geos:Ring Self-intersection at or near point -93.18333333333328 48.698611111111852
WARNING:shapely.geos:Ring Self-intersection at or near point -95.967436930334117 48.322375008519813
WARNING:shapely.geos:Ring Self-intersection at or near point -92.137499999999861 47.2977777777784
WARNING:shapely.geos:Ring Self-intersection at or near point -94.051388888888908 46.466388888889448
WARNING:shapely.geos:Ring Self-intersection at or near point -95.948825819223003 46.116541675186305
WARNING:shapely.geos:Ring Self-intersection at or near point -99.989722222218901 39.784166666667161
WARNING:shapely.geos:Ring Self-intersection at or near point -92.734166666666482 38.09507776472325
WARNING:shapely.geos:Ring Self-intersection at or near point -96.789444444440861 37.858333333333675
WARNING:shapely.geos:Ring Self-intersection at or near point -93.705555555555449 37.589244431389879
WARNING:shapely.geos:Ring Self-intersection at or near point -92.736944444444262 36.470355542500897
WARNING:shapely.geos:Ring Self-intersection at or near point -97.340833333329797 35.337222222222358
WARNING:shapely.geos:Ring Self-intersection at or near point -97.233888888885346 35.227222222222352
WARNING:shapely.geos:Ring Self-intersection at or near point -94.817371410451656 34.932222222222855
WARNING:shapely.geos:Ring Self-intersection at or near point -93.180704743784858 34.613055555556159
WARNING:shapely.geos:Ring Self-intersection at or near point -93.740704743784903 34.173055555556125
WARNING:shapely.geos:Ring Self-intersection at or near point -94.183204743784941 33.337777777778278
WARNING:shapely.geos:Ring Self-intersection at or near point -96.008211263020584 32.984722222222793
WARNING:shapely.geos:Ring Self-intersection at or near point -102.79541666666761 19.329305555555752
WARNING:shapely.geos:Ring Self-intersection at or near point -96.552361111108624 18.274305555555642
WARNING:shapely.geos:Ring Self-intersection at or near point -96.34763888888638 18.033750000000001
WARNING:shapely.geos:Ring Self-intersection at or near point -74.800555555555135 53.550833333333173
WARNING:shapely.geos:Ring Self-intersection at or near point -88.135000000000005 49.378512294558128
WARNING:shapely.geos:Ring Self-intersection at or near point -79.264754460011886 48.442257342868828
WARNING:shapely.geos:Ring Self-intersection at or near point -78.586976682234067 47.586979565090985
WARNING:shapely.geos:Ring Self-intersection at or near point -76.463643348900561 47.369201787313187
WARNING:shapely.geos:Ring Self-intersection at or near point -84.165833334133467 44.66728213310315
WARNING:shapely.geos:Ring Self-intersection at or near point -77.253333333333003 44.093244628906319
WARNING:shapely.geos:Ring Self-intersection at or near point -78.156388888888642 40.374355740017137
WARNING:shapely.geos:Ring Self-intersection at or near point -76.824869249131751 37.367222222222082
WARNING:shapely.geos:Ring Self-intersection at or near point -78.430702582465216 36.609444444444243
WARNING:shapely.geos:Ring Self-intersection at or near point -86.550555555555505 36.07504547860875
WARNING:shapely.geos:Ring Self-intersection at or near point -83.465474378797538 35.950277777777856
WARNING:shapely.geos:Ring Self-intersection at or near point -81.449363267686252 35.774999999999999
WARNING:shapely.geos:Ring Self-intersection at or near point -80.924918823241768 35.581111111111156
WARNING:shapely.geos:Ring Self-intersection at or near point -83.87797437879756 35.370277777777808
WARNING:shapely.geos:Ring Self-intersection at or near point -87.049166666667091 34.700277777777707
WARNING:shapely.geos:Ring Self-intersection at or near point -82.900277777777831 34.396240340675462
WARNING:shapely.geos:Ring Self-intersection at or near point -85.665555555555869 34.242499999999893
WARNING:shapely.geos:Ring Self-intersection at or near point -83.925277777777907 34.337073674008792
WARNING:shapely.geos:Ring Self-intersection at or near point -81.413611111111038 34.081795896230993
WARNING:shapely.geos:Ring Self-intersection at or near point -85.191944444444715 33.012222222222015
WARNING:shapely.geos:Ring Self-intersection at or near point -82.549722222221703 29.039305555555433
WARNING:shapely.geos:Ring Self-intersection at or near point -89.044583333333392 14.045972222222087
WARNING:shapely.geos:Ring Self-intersection at or near point -78.479861111110651 9.1009722222225129
WARNING:shapely.geos:Ring Self-intersection at or near point -79.645694444444203 -0.87208333333300403
WARNING:shapely.geos:Ring Self-intersection at or near point -69.277331068249978 54.440833333333252
WARNING:shapely.geos:Ring Self-intersection at or near point -73.304553290472526 54.127499999999891
WARNING:shapely.geos:Ring Self-intersection at or near point -65.430466280490279 53.960987781679975
WARNING:shapely.geos:Ring Self-intersection at or near point -63.859355169379036 53.587932226124387
WARNING:shapely.geos:Ring Self-intersection at or near point -65.236577391601372 53.581543337235502
WARNING:shapely.geos:Ring Self-intersection at or near point -70.869553290472325 50.779999999999625
WARNING:shapely.geos:Ring Self-intersection at or near point -68.901388889689031 49.696178229650265
WARNING:shapely.geos:Ring Self-intersection at or near point -74.533671569823923 48.545388868243066
WARNING:shapely.geos:Ring Self-intersection at or near point -73.905060458712768 46.725666646020699
WARNING:shapely.geos:Ring Self-intersection at or near point -71.906666666666425 44.348055555555568
WARNING:shapely.geos:Ring Self-intersection at or near point -69.854444444444098 44.179166666667335
WARNING:shapely.geos:Ring Self-intersection at or near point -74.184444444444381 41.956666666666486
WARNING:shapely.geos:Ring Self-intersection at or near point -63.026249999999628 6.9473611111112437
WARNING:shapely.geos:Ring Self-intersection at or near point -74.986249999999671 6.3209722222223608
WARNING:shapely.geos:Ring Self-intersection at or near point -57.309248491750338 48.351111111111628
WARNING:shapely.geos:Ring Self-intersection at or near point -55.764526269527991 48.175000000000502
WARNING:shapely.geos:Ring Self-intersection at or near point -54.959027777777266 4.6562500000003029
WARNING:shapely.geos:Ring Self-intersection at or near point -59.761805555555149 -1.781527777777568
WARNING:shapely.geos:Ring Self-intersection at or near point -49.457916666666286 -4.6684722222222019
WARNING:shapely.geos:Ring Self-intersection at or near point -55.726249999999212 -14.996805555555556
WARNING:shapely.geos:Ring Self-intersection at or near point -48.25486111111119 -20.083749999999515
WARNING:shapely.geos:Ring Self-intersection at or near point -48.712361111111228 -20.153472222221744
WARNING:shapely.geos:Ring Self-intersection at or near point -45.801527777777665 -20.562638888888443
WARNING:shapely.geos:Ring Self-intersection at or near point -52.323750000000224 -21.868750000000155
WARNING:shapely.geos:Ring Self-intersection at or near point -50.467916666666746 -27.859583333332939
WARNING:shapely.geos:Ring Self-intersection at or near point -36.872361111110521 -5.7040277777776547
WARNING:shapely.geos:Ring Self-intersection at or near point -8.1120833333332261 38.19513888888936
WARNING:shapely.geos:Ring Self-intersection at or near point -8.1695833333331471 11.231249999999756
WARNING:shapely.geos:Ring Self-intersection at or near point -7.02458333333312 6.3481944444445304
WARNING:shapely.geos:Ring Self-intersection at or near point -3.195138888888768 5.6065277777778473
WARNING:shapely.geos:Ring Self-intersection at or near point 12.794722222222402 59.046250000000647
WARNING:shapely.geos:Ring Self-intersection at or near point 4.6262500000005247 9.804305555555894
WARNING:shapely.geos:Ring Self-intersection at or near point 13.875138888889312 8.8770833333336405
WARNING:shapely.geos:Ring Self-intersection at or near point 13.080694444444804 6.510694444444562
WARNING:shapely.geos:Ring Self-intersection at or near point -0.76569444444435097 8.5931944444446859
WARNING:shapely.geos:Ring Self-intersection at or near point 0.096805555555718001 6.1604166666667126
WARNING:shapely.geos:Ring Self-intersection at or near point 11.345138888889108 6.0893055555556392
WARNING:shapely.geos:Ring Self-intersection at or near point 28.605682373046989 61.652885437011776
WARNING:shapely.geos:Ring Self-intersection at or near point 27.28250000000029 54.542638888889257
WARNING:shapely.geos:Ring Self-intersection at or near point 27.04402777777814 -10.88263888888812
WARNING:shapely.geos:Ring Self-intersection at or near point 28.296249999999812 -16.54541666666606
WARNING:shapely.geos:Ring Self-intersection at or near point 25.605694444444826 -30.643749999999624
WARNING:shapely.geos:Ring Self-intersection at or near point 37.719166666667334 58.414861111111719
WARNING:shapely.geos:Ring Self-intersection at or near point 30.781388888889062 50.17513888888859
WARNING:shapely.geos:Ring Self-intersection at or near point 32.360138888889523 49.482083333334131
WARNING:shapely.geos:Ring Self-intersection at or near point 45.917361111111326 51.333194444445084
WARNING:shapely.geos:Ring Self-intersection at or near point 33.697638888889628 48.929583333334087
WARNING:shapely.geos:Ring Self-intersection at or near point 38.437361111111102 37.589583333333906
WARNING:shapely.geos:Ring Self-intersection at or near point 32.869861111111639 23.719027777777864
WARNING:shapely.geos:Ring Self-intersection at or near point 31.911805555555599 -2.7223611111109598
WARNING:shapely.geos:Ring Self-intersection at or near point 30.706250000000001 -15.632638888888094
WARNING:shapely.geos:Ring Self-intersection at or near point 31.080031703795669 -24.867639990434554
WARNING:shapely.geos:Ring Self-intersection at or near point 53.135833333333437 55.789861111111293
WARNING:shapely.geos:Ring Self-intersection at or near point 73.345972222222827 25.716527777778168
WARNING:shapely.geos:Ring Self-intersection at or near point 73.797916666667177 21.274583333333773
WARNING:shapely.geos:Ring Self-intersection at or near point 74.400138888888677 18.777083333333191
WARNING:shapely.geos:Ring Self-intersection at or near point 81.622638888888687 22.424305555555637
WARNING:shapely.geos:Ring Self-intersection at or near point 79.257916666666674 21.670972222222467
WARNING:shapely.geos:Ring Self-intersection at or near point 77.320138888889531 17.832361111111339
WARNING:shapely.geos:Ring Self-intersection at or near point 77.877361111110872 11.907916666667006
WARNING:shapely.geos:Ring Self-intersection at or near point 76.552083333333002 11.389027777777933
WARNING:shapely.geos:Ring Self-intersection at or near point 76.809861111110806 10.380138888888963
WARNING:shapely.geos:Ring Self-intersection at or near point 104.50138888808591 30.192222222222572
WARNING:shapely.geos:Ring Self-intersection at or near point 92.624140258749307 25.502499999999706
WARNING:shapely.geos:Ring Self-intersection at or near point 91.818055555556356 23.493333333333943
WARNING:shapely.geos:Ring Self-intersection at or near point 95.516944310930469 20.391666666666477
WARNING:shapely.geos:Ring Self-intersection at or near point 100.55194444444115 17.811666666667005
WARNING:shapely.geos:Ring Self-intersection at or near point 99.420833333333235 15.561388888888931
WARNING:shapely.geos:Ring Self-intersection at or near point 98.526040858806013 15.043888888888926
WARNING:shapely.geos:Ring Self-intersection at or near point 100.97816813151019 13.194444444444468
WARNING:shapely.geos:Ring Self-intersection at or near point 98.562777777778365 9.0716666666669994
WARNING:shapely.geos:Ring Self-intersection at or near point 102.56166666666353 4.9400175214524538
WARNING:shapely.geos:Ring Self-intersection at or near point 111.16152777777565 34.824305555556165
WARNING:shapely.geos:Ring Self-intersection at or near point 114.85958333333039 29.79430555555572
WARNING:shapely.geos:Ring Self-intersection at or near point 114.84541666666372 29.607638888889038
WARNING:shapely.geos:Ring Self-intersection at or near point 119.10763888889144 29.541250000000804
WARNING:shapely.geos:Ring Self-intersection at or near point 115.31513888889114 29.28708333333412
WARNING:shapely.geos:Ring Self-intersection at or near point 115.80894388077891 27.827279005573374
WARNING:shapely.geos:Ring Self-intersection at or near point 112.67381123742808 24.954557156050186
WARNING:shapely.geos:Ring Self-intersection at or near point 115.6075689239198 23.944379676611096
WARNING:shapely.geos:Ring Self-intersection at or near point 114.51347222222702 23.880972222222866
WARNING:shapely.geos:Ring Self-intersection at or near point 109.57763888889153 22.0554167833245
WARNING:shapely.geos:Ring Self-intersection at or near point 104.95027777777372 21.843611111111585
WARNING:shapely.geos:Ring Self-intersection at or near point 105.20375000000229 20.727361227768839
WARNING:shapely.geos:Ring Self-intersection at or near point 109.6566666666637 19.304444444444364
WARNING:shapely.geos:Ring Self-intersection at or near point 105.37055555555224 14.952777777777348
WARNING:shapely.geos:Ring Self-intersection at or near point 107.03711554633128 11.116944444444616
WARNING:shapely.geos:Ring Self-intersection at or near point 111.91244120279792 1.173540701230178
WARNING:shapely.geos:Ring Self-intersection at or near point 107.36629846360907 -6.5484800160109327
WARNING:shapely.geos:Ring Self-intersection at or near point 125.58013888889116 39.996805555555333
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
WARNING:shapely.geos:Self-intersection at or near point 131.49400972199999 34.357413998378128
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
WARNING:shapely.geos:Ring Self-intersection at or near point 131.43357769981174 33.138324618101549
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
WARNING:shapely.geos:Self-intersection at or near point 130.060966111 32.817799444400002
WARNING:shapely.geos:Ring Self-intersection at or near point 130.49771414117069 31.993888722847824
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
WARNING:shapely.geos:Self-intersection at or near point 137.74883159917832 36.155175357375548
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
ERROR:shapely.geos:IllegalArgumentException: Points of LinearRing do not form a closed linestring
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
Error(Could not create geometry because of errors while reading input.), skipping geometry.
WARNING:shapely.geos:Ring Self-intersection at or near point 147.48083333333352 -9.5063914908505218
WARNING:shapely.geos:Ring Self-intersection at or near point 149.14694444444922 -35.30916666666608
WARNING:shapely.geos:Ring Self-intersection at or near point 151.76338934878515 -24.772866925011328


In [ ]:
import fiona

c = fiona.collection(files[0], "r")
print(c.schema)

type_map = { 'float' : ogr.OFTReal, 'int' : ogr.OFTInteger }

fields = c.schema['properties']
for p in fields:
    print p + ' ' + fields[p]

In [3]:
# Select sub-basins that fall within a larger basin and then group them by using a regular grid

# Amazon
# main_basin_shp = 'D:\gis\HydroBASINS\without_lakes\hybas_sa_lev02_v1c.shp'
# main_basin_id = 6020006540
# sub_basin_shp = 'D:\gis\HydroBASINS\without_lakes\hybas_sa_lev06_v1c.shp'

# Murray & Darling
main_basin_shp = 'D:\gis\HydroBASINS\without_lakes\hybas_au_lev03_v1c.shp'
main_basin_id = 5030073410
sub_basin_shp = 'D:\gis\HydroBASINS\without_lakes\hybas_au_lev07_v1c.shp'

pylab.rcParams['figure.figsize'] = (17.0, 20.0)

fig = plt.figure()
axes = plt.axes()
axes.set_aspect('equal', 'datalim')

main_basin_geom = None

with fiona.collection(main_basin_shp, "r") as input:
    for f in input:
        geom = sl.geometry.shape(f['geometry'])

        id = f['properties']['HYBAS_ID']
        if id == main_basin_id:
            main_basin_geom = geom
        draw(geom)

bounds = main_basin_geom.bounds
axes.set_xlim(bounds[0], bounds[2])
axes.set_ylim(bounds[1], bounds[3])
                                
with fiona.collection(sub_basin_shp, "r") as input:
    for f in input:
        geom = sl.geometry.shape(f['geometry'])

        if main_basin_geom.contains(geom):
            draw(geom, fill='#aaaaff', alpha=0.5)
        else:
            pass
            # draw(geom, fill='#ffffff', lw=0.1)

plt.show()


C:\Anaconda\lib\site-packages\numba\dataflow.py:52: RuntimeWarning: inconsistent stack offset for block(offset:326, outgoing: [], incoming: [21, 61, 110, 323])
  RuntimeWarning)

In [7]:
def get_basin(id, path):
    with fiona.collection(path, "r") as input:
        for f in input:
            id = f['properties']['HYBAS_ID']
            if id == main_basin_id:
                return f
    return None

def get_child_basins(parent_pfaf_id, child_basins_path):
    with fiona.collection(child_basins_path, "r") as input:
        for feature in input:
            pfaf_id = feature['properties']['PFAF_ID'] 
            if str(pfaf_id).startswith(str(parent_pfaf_id)): 
                yield feature

# Amazon
# files = glob.glob(r'D:\gis\HydroBASINS\without_lakes\hybas_sa_lev*_v1c.shp')[1:]
# main_basin_id = 6020006540
# start = 0

files = glob.glob(r'D:\gis\HydroBASINS\without_lakes\hybas_au_lev*_v1c.shp')[1:]
main_basin_id = 5030073410
start = 2

main_basin = get_basin(main_basin_id, files[start])

pfaf_id = str(main_basin['properties']['PFAF_ID'])
main_basin_geom = sl.geometry.shape(main_basin['geometry'])

print('Main basin: ' + files[0])
print('Area: ' + str(main_basin_geom.area))

total_count = 0

def get_children(main_basin_id, child_basins_path):
    print('Child basins: ' + child_basins_path)
    child_basins = list(get_child_basins(main_basin_id, child_basins_path))
    child_basins_ids = [str(f['properties']['PFAF_ID']) for f in child_basins]
    child_basins_geoms = [sl.geometry.shape(f['geometry']) for f in child_basins]
    count = len(child_basins)
    global total_count
    total_count += count
    print(count)
    
    return (child_basins, child_basins_ids, child_basins_geoms)

basins3, ids3, geoms3 = get_children(pfaf_id, files[start + 1])
basins4, ids4, geoms4 = get_children(pfaf_id, files[start + 2])
basins5, ids5, geoms5 = get_children(pfaf_id, files[start + 3])
basins6, ids6, geoms6 = get_children(pfaf_id, files[start + 4])
basins7, ids7, geoms7 = get_children(pfaf_id, files[start + 5])
basins8, ids8, geoms8 = get_children(pfaf_id, files[start + 6])

basins_per_level = [basins3, basins4, basins5, basins6, basins7, basins8]

print('Total basins: ' + str(total_count))


Main basin: D:\gis\HydroBASINS\without_lakes\hybas_au_lev01_v1c.shp
Area: 100.489618539
Child basins: D:\gis\HydroBASINS\without_lakes\hybas_au_lev04_v1c.shp
10
Child basins: D:\gis\HydroBASINS\without_lakes\hybas_au_lev05_v1c.shp
27
Child basins: D:\gis\HydroBASINS\without_lakes\hybas_au_lev06_v1c.shp
137
Child basins: D:\gis\HydroBASINS\without_lakes\hybas_au_lev07_v1c.shp
434
Child basins: D:\gis\HydroBASINS\without_lakes\hybas_au_lev08_v1c.shp
1439
Child basins: D:\gis\HydroBASINS\without_lakes\hybas_au_lev09_v1c.shp
3809
Total basins: 5856

Put all catchments into a graph


In [8]:
import networkx as nx
n = nx.DiGraph()

ids = [pfaf_id]+ids3+ids4+ids5+ids6+ids7
areas = [g.area for g in [main_basin_geom]+geoms3+geoms4+geoms5+geoms6+geoms7]
basins = [b for b in [main_basin]+basins3+basins4+basins5+basins6+basins7]

n.add_nodes_from([(i,{'area':a, 'feature':f}) for (i,a,f) in zip(ids, areas, basins)])
# n.add_nodes_from([(i,{'area':a}) for (i,a,f) in zip(ids, areas, basins)])

n.add_edges_from([(pfaf_id, i) for i in ids3])
n.add_edges_from([(i[:-1], i) for i in ids4])
n.add_edges_from([(i[:-1], i) for i in ids5])
n.add_edges_from([(i[:-1], i) for i in ids6])
n.add_edges_from([(i[:-1], i) for i in ids7])

# max_level = 6
# n = n.subgraph([kv[0] for kv in list(n.degree().iteritems()) if len(kv[0]) < max_level + 2])

In [9]:
nodes=dict(n.nodes(data=True))

shape=sl.geometry.shape(nodes['5643001']['feature']['geometry'])
shape


Out[9]:

In [10]:
def get_wh(geom):
    bounds = geom.bounds

    w = bounds[2]-bounds[0]
    h = bounds[3]-bounds[1]
    
    return (w, h) 
    
(w, h) = get_wh(shape)

print('w', w)
print('h', h)


('w', 2.4458333333333258)
('h', 1.3172339545355882)

In [11]:
def geometry_larger_than(geom, w, h):
    bounds = geom.bounds
    gw = bounds[2]-bounds[0]
    gh = bounds[3]-bounds[1]
    return (gw > w) or (gh > h)

def geometry_smaller_than(geom, w, h):
    bounds = geom.bounds
    gw = bounds[2]-bounds[0]
    gh = bounds[3]-bounds[1]
    return (gw < w) and (gh < h)

def get_leaves_smaller_than(node, w, h):
    for child_node in n.successors(node):
        geom = sl.geometry.shape(nodes[child_node]['feature']['geometry'])
        
        last = len(n.successors(child_node)) == 0
    
        if(geometry_smaller_than(geom, w, h) or last):
            yield child_node
        else:
            for child_child_node in get_leaves_smaller_than(child_node, w, h):
                yield child_child_node

In [65]:
max_w = 2.0
max_h = 2.0

In [66]:
leaves = list(get_leaves_smaller_than('564', max_w, max_h))
print(len(leaves))
s = n.subgraph(leaves)


469

In [53]:
basins_wh = []
basins_wh_min = []
basins_wh_max = []
basins_count = []


for basins in basins_per_level:
    whs = [get_wh(sl.geometry.shape(basin['geometry'])) for basin in basins]
    max_wh = np.max([max(whi) for whi in whs])
    min_wh = np.min([min(whi) for whi in whs])
    count = len(basins)
    wh = 0.5*(max_wh+min_wh)
    basins_wh_min.append(min_wh)
    basins_wh_max.append(max_wh)
    basins_wh.append(wh)
    basins_count.append(count)
    print(min_wh, wh, max_wh, count)


(1.4041666666666401, 5.3000028822156935, 9.1958390977647468, 10)
(0.021097988552526203, 3.0334685431586337, 6.0458390977647412, 27)
(0.021097988552526203, 2.4917989942762677, 4.9625000000000092, 137)
(0.0041666666666628771, 1.7315336439344549, 3.4589006212022468, 434)
(0.0041666666666628771, 1.1895833333333314, 2.375, 1439)
(0.0041666666666628771, 0.75236697726779944, 1.500567287868936, 3809)

In [59]:
# show number of catchments per level + maximum Width / Height
min_wh_values = []
max_wh_values = []
wh_values = []
average_wh_values = []
catchment_count = []
for i in np.arange(0.2, 9.0, 0.2):
    selection = list(get_leaves_smaller_than('564', i, i))
    
    wh = [get_wh(sl.geometry.shape(nodes[basin]['feature']['geometry'])) for basin in selection]
    max_wh = np.max([max(whi) for whi in wh])
    min_wh = np.min([min(whi) for whi in wh])
    average_wh_values.append(0.5*(min_wh+max_wh))
    count = len(basins)

    count = len(selection)
    print(min_wh, i, max_wh, count)
    catchment_count.append(count)
    max_wh_values.append(max_wh)
    min_wh_values.append(min_wh)
    
    wh_values.append(i)


(0.0041666666666628771, 0.20000000000000001, 2.375, 1439)
(0.0041666666666628771, 0.40000000000000002, 2.375, 1439)
(0.0041666666666628771, 0.60000000000000009, 2.375, 1439)
(0.0041666666666628771, 0.80000000000000004, 2.375, 1402)
(0.0041666666666628771, 1.0, 2.375, 1180)
(0.0041666666666628771, 1.2, 2.375, 1013)
(0.0041666666666628771, 1.4000000000000001, 2.375, 796)
(0.0041666666666628771, 1.6000000000000001, 2.375, 672)
(0.0041666666666628771, 1.8, 2.375, 546)
(0.0041666666666628771, 2.0, 2.375, 469)
(0.0041666666666628771, 2.2000000000000002, 2.375, 411)
(0.0041666666666628771, 2.4000000000000004, 2.375, 346)
(0.0041666666666628771, 2.6000000000000005, 2.5256442599826414, 305)
(0.0041666666666628771, 2.8000000000000003, 2.7203645494249145, 272)
(0.0041666666666628771, 3.0000000000000004, 2.9875000000000114, 222)
(0.0041666666666628771, 3.2000000000000002, 3.1883319430881158, 173)
(0.0041666666666628771, 3.4000000000000004, 3.2166666666666757, 165)
(0.0041666666666628771, 3.6000000000000005, 3.5790659586588447, 141)
(0.021097988552526203, 3.8000000000000003, 3.6041666666666856, 135)
(0.021097988552526203, 4.0, 3.8964006212022575, 110)
(0.021097988552526203, 4.2000000000000002, 3.8964006212022575, 110)
(0.021097988552526203, 4.4000000000000004, 4.3342163085937493, 78)
(0.021097988552526203, 4.6000000000000005, 4.3342163085937493, 78)
(0.021097988552526203, 4.8000000000000007, 4.6958333333333258, 69)
(0.021097988552526203, 5.0000000000000009, 4.9625000000000092, 60)
(0.021097988552526203, 5.2000000000000002, 4.9625000000000092, 60)
(0.021097988552526203, 5.4000000000000004, 5.2928892347547958, 52)
(0.021097988552526203, 5.6000000000000005, 5.2928892347547958, 52)
(0.021097988552526203, 5.8000000000000007, 5.7420199924045221, 44)
(0.021097988552526203, 6.0000000000000009, 5.8130672878689218, 36)
(0.021097988552526203, 6.2000000000000002, 6.0458390977647412, 27)
(0.021097988552526203, 6.4000000000000004, 6.0458390977647412, 27)
(0.021097988552526203, 6.6000000000000005, 6.0458390977647412, 27)
(0.18425394694011032, 6.8000000000000007, 6.6583333333333243, 19)
(0.18425394694011032, 7.0000000000000009, 6.6583333333333243, 19)
(0.18425394694011032, 7.2000000000000002, 6.6583333333333243, 19)
(0.18425394694011032, 7.4000000000000004, 6.6583333333333243, 19)
(0.18425394694011032, 7.6000000000000005, 6.6583333333333243, 19)
(0.18425394694011032, 7.8000000000000007, 6.6583333333333243, 19)
(0.18425394694011032, 8.0, 6.6583333333333243, 19)
(0.18425394694011032, 8.1999999999999993, 6.6583333333333243, 19)
(0.18425394694011032, 8.4000000000000004, 6.6583333333333243, 19)
(0.18425394694011032, 8.5999999999999996, 6.6583333333333243, 19)
(0.18425394694011032, 8.7999999999999989, 6.6583333333333243, 19)

In [60]:
plt.plot(catchment_count, wh_values, 'ob-', linewidth=2.0)
plt.plot(catchment_count, average_wh_values, 'b-')
plt.plot(catchment_count, min_wh_values, 'b--')
plt.plot(catchment_count, max_wh_values, 'b--')


# plt.errorbar(wh_values, catchment_count,xerr=[min_error_wh_values, max_error_wh_values], fmt='.')


# plt.errorbar(wh_values, catchment_count, xerr=error_wh_values)

plt.plot(basins_count, basins_wh, 'or-', linewidth=2.0)
plt.plot(basins_count, basins_wh_min, 'r--')
plt.plot(basins_count, basins_wh_max, 'r--')

# plt.errorbar(basins_wh, basins_count, xerr=[basins_wh_error_min,basins_wh_error_max])


Out[60]:
[<matplotlib.lines.Line2D at 0x1b635c18>]

In [198]:
large_nodes = list((node for (node, data) in n.nodes(data=True) 
                    if geometry_larger_than(sl.geometry.shape(data['feature']['geometry']), max_w, max_h)))
s = n.subgraph(large_nodes)
print(len(large_nodes))
leaves=[node for node,data in s.out_degree().items() if data==0]
print(len(list(leaves)))
s = n.subgraph(leaves)


136
Out[198]:
'5640017'

In [14]:
nx.write_dot(n, r"..\output\catchment_pfafids.dot")
print(len(n.nodes()))


2048

In [69]:
pos = nx.spring_layout(nx.Graph(n), iterations=500)

In [70]:
import math

node_areas=[data['area'] for node,data in n.nodes(data=True)]

node_sizes=[math.log(a)*100 for a in node_areas]

nx.draw_networkx_nodes(n, pos, node_color=node_areas, node_size=node_sizes, cmap=plt.cm.Blues, alpha=0.85)
_ = nx.draw_networkx_edges(n, pos, alpha=0.2)
# _ = nx.draw_networkx_labels(n, pos, font_size=15)

# nx.draw_networkx_nodes(s, pos, node_size=20, cmap=plt.cm.Reds, alpha=0.85)


C:\Anaconda\lib\site-packages\matplotlib\collections.py:764: RuntimeWarning: invalid value encountered in sqrt
  scale = np.sqrt(self._sizes) * dpi / 72.0 * self._factor

In [335]:
leaf_basins = list(s.nodes(data=True))

fig = plt.figure()
axes = plt.axes()
axes.set_aspect('equal', 'datalim')

for f in leaf_basins:
    geom = sl.geometry.shape(f[1]['feature']['geometry'])
    draw(geom, alpha=0.5)

plt.show()