from matplotlib import pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch

import numpy
import geopandas as gpd
import shapely
from shapely.wkt import loads
import cartagen as c4

lines = [
    loads('LineString (-187744.55296041004476137 5374879.20447004493325949, -187741.77143371716374531 5374866.09463881049305201, -187738.47241368613322265 5374847.8529986385256052, -187735.94963366238516755 5374821.84895839355885983, -187735.1733936551027 5374804.57761823106557131, -187735.56151365875848569 5374779.54387799557298422, -187737.8902336806349922 5374753.73389775305986404)'),
    loads('LineString (-187785.26243634850834496 5374925.43387492373585701, -187800.57161427021492273 5374917.1324192900210619, -187813.76769439433701336 5374908.98189921304583549, -187827.15783452027244493 5374899.47295912355184555, -187837.8311346206755843 5374889.76995903253555298, -187849.47473473017453216 5374876.76793891005218029, -187857.0430748013604898 5374863.7659187875688076, -187865.38765487985801883 5374848.04705864004790783, -187870.04509492366923951 5374834.65691851451992989)'),
    loads('LineString (-187601.03480794903589413 5375111.73147667571902275, -187603.79477241938002408 5375124.48553123883903027, -187607.2878524522529915 5375137.77864136453717947, -187613.10965250700246543 5375151.84799149632453918, -187620.48393257637508214 5375164.1708016125485301, -187629.3136626594059635 5375176.00846172403544188, -187639.8899327589024324 5375185.90552181657403708, -187657.16127292133751325 5375199.48972194455564022)'),
    loads('LineString (-187133.91082133317831904 5375205.57026867009699345, -187113.69623780969413929 5375215.49967209622263908, -187082.64663751763873734 5375226.17297219671308994, -187050.04455721098929644 5375233.35319226421415806, -187009.87413683318300173 5375237.81657230667769909, -186969.50965645350515842 5375238.39875231217592955, -186925.84615604282589629 5375233.54725226573646069)'),
    loads('LineString (-187372.99274358194088563 5374917.32647929154336452, -187354.52470007486408576 5374904.0333691667765379, -187334.53651988686760888 5374890.64322904124855995, -187313.96615969340200536 5374878.61150892730802298, -187289.90271946706343442 5374866.19166881125420332, -187257.30063916041399352 5374851.63716867379844189)'),
]

network = gpd.GeoDataFrame(geometry=lines)
polygon = loads('Polygon ((-187729.7181513816176448 5375043.76735159195959568, -187713.50336011801846325 5375059.29215173795819283, -187683.83374206116423011 5375082.06185861956328154, -187644.50424835790181533 5375100.00162767712026834, -187601.03480794903589413 5375111.73147667571902275, -187565.84526095140608959 5375108.97151220589876175, -187523.06581166014075279 5375100.6916187945753336, -187478.90638013367424719 5375104.14157438278198242, -187430.60700190160423517 5375113.11145891156047583, -187365.05784572951961309 5375140.71110361535102129, -187333.31825431986362673 5375160.72084602620452642, -187301.57866291023674421 5375184.87053514178842306, -187274.66900932378484868 5375205.57026867009699345, -187241.54943567892769352 5375226.27000219747424126, -187200.83995974049321376 5375236.61986896209418774, -187169.1003683308372274 5375234.54989560879766941, -187149.78061703802086413 5375222.13005549181252718, -187133.91082133317831904 5375205.57026867009699345, -187128.39089239237364382 5375186.94050849415361881, -187131.84084798037656583 5375168.31074831914156675, -187147.02065256761852652 5375146.23103255592286587, -187167.7203860956360586 5375124.84130791015923023, -187189.11011074125417508 5375107.59152997005730867, -187216.01976432770607062 5375095.86168097145855427, -187247.06936461976147257 5375090.34175203088670969, -187285.70886720539419912 5375089.65176091250032187, -187320.89841420305310749 5375084.82182308938354254, -187342.28813884867122397 5375073.09197409078478813, -187362.98787237671785988 5375051.7022494450211525, -187378.16767696393071674 5375033.07248927000910044, -187382.82511700777104124 5375017.84956023655831814, -187381.79013033135561273 5374999.86666673421859741, -187376.61519694936578162 5374979.03755987156182528, -187367.30031686174334027 5374959.24343968648463488, -187365.23034350891248323 5374936.99122614320367575, -187372.99274358194088563 5374917.32647929154336452, -187394.20997044816613197 5374903.35415916051715612, -187417.49717066719313152 5374901.80167914554476738, -187442.33685090084327385 5374918.36146596819162369, -187454.23919767944607884 5374944.23613287881016731, -187468.72901114908745512 5374962.86589305382221937, -187492.53370470632216893 5374983.04813324380666018, -187517.54588271933607757 5374999.95291562471538782, -187541.00558071775594726 5375008.23280903603881598, -187578.95509218581719324 5375010.30278238840401173, -187603.79477241943823174 5375003.40287121292203665, -187629.32444377068895847 5374989.60304886102676392, -187654.16412400430999696 5374967.52333309780806303, -187675.55384864995721728 5374941.30367062892764807, -187701.77351111877942458 5374920.60393710061907768, -187725.2332091172283981 5374886.79437233787029982, -187744.55296041004476137 5374879.20447004493325949, -187760.42275611485820264 5374886.79437233787029982, -187776.98254293730133213 5374901.97417692560702562, -187785.26243634850834496 5374925.43387492373585701, -187785.95242746610892937 5374946.82359956949949265, -187776.98254293730133213 5374979.25318209640681744, -187757.66279164445586503 5375010.30278238840401173, -187743.51797373362933286 5375027.89755588676780462, -187729.7181513816176448 5375043.76735159195959568))')

skeleton = c4.skeletonize_network(polygon, network, sigma=10, blend_smoothing=True)

fig = plt.figure(1, (8, 8))

sub1 = fig.add_subplot(211)
sub1.set_aspect('equal')
sub1.set_title('a) The provided polygon with incoming network', pad=10, family='sans-serif')
sub1.axes.get_xaxis().set_visible(False)
sub1.axes.get_yaxis().set_visible(False)

poly1 = Path.make_compound_path(Path(numpy.asarray(polygon.exterior.coords)[:, :2]),*[Path(numpy.asarray(ring.coords)[:, :2]) for ring in polygon.interiors])
sub1.add_patch(PathPatch(poly1, facecolor="lightgray", edgecolor='none'))

for line in lines:
    path1 = Path(numpy.asarray(line.coords)[:, :2])
    sub1.add_patch(PathPatch(path1, facecolor="none", edgecolor='gray', linewidth=1))

sub1.autoscale_view()

#########################################################

sub2 = fig.add_subplot(212)
sub1.set_aspect('equal')
sub2.set_title('b) The blended skeleton with a blended gaussian smoothing (sigma=10)', pad=10, family='sans-serif')
sub2.axes.get_xaxis().set_visible(False)
sub2.axes.get_yaxis().set_visible(False)

poly1 = Path.make_compound_path(Path(numpy.asarray(polygon.exterior.coords)[:, :2]),*[Path(numpy.asarray(ring.coords)[:, :2]) for ring in polygon.interiors])
sub2.add_patch(PathPatch(poly1, facecolor="lightgray", edgecolor='none'))

for line in skeleton:
    path1 = Path(numpy.asarray(line['geometry'].coords)[:, :2])
    sub2.add_patch(PathPatch(path1, facecolor="none", edgecolor='red', linewidth=1))

sub2.autoscale_view()

plt.tight_layout()
plt.show()