<jupyter_start><jupyter_text>pyflwdir 測試主要是看能不能使用上頭前溪的資料，看會不會用，能不能達到專案的需求雖然還有問題，但似乎我們專案要用的東西，具備了目前可以產生頭前溪的1.線圖2.子流域3.算一點上游的集水區4.算兩點之間的集水區5.取得通過幾點下游的線圖<jupyter_code>#畫圖要用的函示
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm, colors
import cartopy.crs as ccrs
import descartes
import numpy as np
np.random.seed(seed=101)
matplotlib.rcParams['savefig.bbox'] = 'tight'
matplotlib.rcParams['savefig.dpi'] = 256
plt.style.use('seaborn-whitegrid')

def quickplot(gdfs=[], maps=[], hillshade=True, title='', filename='flw', save=False):
    fig = plt.figure(figsize=(8,15))
    ax = fig.add_subplot(projection=ccrs.PlateCarree())
    # plot hillshade background
    if hillshade:
        ls = matplotlib.colors.LightSource(azdeg=115, altdeg=45)
        hillshade = ls.hillshade(np.ma.masked_equal(elevtn, -9999), vert_exag=1e3)
        ax.imshow(hillshade, origin='upper', extent=flw.extent, cmap='Greys', alpha=0.3, zorder=0)
    # plot geopandas GeoDataFrame
    for gdf, kwargs in gdfs:
        gdf.plot(ax=ax, **kwargs)
    for data, nodata, kwargs in maps:
        ax.imshow(np.ma.masked_equal(data, nodata), origin='upper', extent=flw.extent, **kwargs)
    ax.set_aspect('equal')
    ax.set_title(title, fontsize='large')
    ax.text(0.01, 0.01, 'created with pyflwdir', transform=ax.transAxes, fontsize='large')
    if save:
        plt.savefig(f'output/{filename}.png')
    return ax
# 載入資料
import rasterio
import pyflwdir
import geopandas as gpd

my=1

if my:
    basin_id='1300'
    flwdir_filename ='data/catchment/C%s_20m_ldd.tif' %(basin_id)
    elevtn_filename = 'data/catchment/C%s_20m_elv0.tif' %(basin_id)

else:
    flwdir_filename = "data/catchment/rhine_d8.tif"
    elevtn_filename = 'data/catchment/rhine_elv0.tif'

with rasterio.open(flwdir_filename, 'r') as src1:
    flwdir = src1.read(1)
    #flwdir[flwdir==-32768]=-9999 # debug
    #flwdir[flwdir==-1]=-999 # debug
    transform = src1.transform
    crs = src1.crs
    latlon = crs.to_epsg() == 4326
    print("%s info:%s" %(flwdir_filename,src1))
with rasterio.open(elevtn_filename, 'r') as src2:
    elevtn = src2.read(1)
    if my:
        elevtn[elevtn==-99999]=-9999
    prof = src2.profile
    print("%s info:%s" %(elevtn_filename,src2))

# 產生 pyflwdir 主要的物件
if my:
    ftype='ldd'
else:
    ftype='d8'
flw = pyflwdir.from_array(flwdir,ftype=ftype, transform=transform, latlon=latlon, cache=True) #d8
print(flw)
# calculate upstream area
uparea = flw.upstream_area(unit='km2')
# plot with lognorm scale
if my:
    upa = (uparea, -9999, dict(cmap='Greens', alpha=0.8,norm=colors.Normalize(vmin=1,vmax=1e5))) #-9999
else:
    upa = (uparea, -9999, dict(cmap='Greens', alpha=1.0, norm=colors.LogNorm(vmin=1, vmax=1e5))) #-9999
    ax = quickplot(maps=[upa], hillshade=True, title='Upstream area', filename='flw_uparea') #True
#河道
feats = flw.streams(strord=None,min_sto=9)
gdf = gpd.GeoDataFrame.from_features(feats, crs=crs)

if my:
    pass
else:
    # plot
    # create nice colormap of Blues with less white
    cmap_streams = colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7)))
    streams = (gdf, dict(column='strord', cmap=cmap_streams))
    ax = quickplot(gdfs=[streams], title='Streams', filename='flw_streams')
# 函式
# we can use rasterio.features to vectorize basins
from rasterio import features

def vectorize(data, nodata, transform, crs=crs):
    feats_gen = features.shapes(
        data, mask=data!=nodata, transform=transform, connectivity=8,
    )
    feats = [
        {"geometry": geom, "properties": {"value": val}}
        for geom, val in list(feats_gen)
    ]

    # parse to geopandas for plotting / writing to file
    gdf = gpd.GeoDataFrame.from_features(feats, crs=crs)
    return gdf
# 算子流域
# define output locations
if my:
    x, y = np.array([260993, 261000]), np.array([2735861, 2736000])
else:
    x, y = np.array([4.67916667, 7.60416667]), np.array([51.72083333, 50.3625])
    
# delineate subbasins
subbasins = flw.basins(xy=(x,y), streams=flw.stream_order()>=4)
# vectorize
gdf_bas = vectorize(subbasins.astype(np.int32), 0, flw.transform)
# plot
streams = (gdf[gdf['strord']>=6], dict(color='grey'))
bas = (gdf_bas, dict(edgecolor='black', facecolor='none', linewidth=0.8))
subbas = (subbasins, 0, dict(cmap='Set3', alpha=0.5))
if not my:
    ax = quickplot([streams, bas], [subbas], title='Basins from point outlets', filename='flw_basins')

# calculate subbasins with a minimum stream order 7
subbas = flw.subbasins_streamorder(min_sto=7, mask=None)
gdf_subbas = vectorize(subbas.astype(np.int32), 0, flw.transform)

#plot
if my:
    gdf_subbas.plot()
    gdf_subbas.to_file('output/river_c1300_subbas.geojson', driver='GeoJSON')
    
else:
    cmap = colors.ListedColormap(cm.Set3(np.random.rand(68))) # create a nice random colormap
    bas1 = (gdf_subbas, dict(edgecolor='black', facecolor='none', linewidth=0.8))
    subbas1 = (subbas, 0, dict(cmap=cmap, alpha=0.5))
    title='Subbasins based on a minimum stream order'
    ax = quickplot([streams, bas1], [subbas1], title=title, filename='flw_subbasins')
# 有 Error , 目前不會解決
# get the first level nine pfafstetter basins
pfafbas1 = flw.subbasins_pfafstetter(depth=1, uparea=uparea)
gdf_pfaf1 = vectorize(pfafbas1.astype(np.int32), 0, flw.transform)

# plot
if my:
    gdf_pfaf1.plot()
    
else:
    pfafbas = (pfafbas1, 0, dict(cmap='Set3', alpha=0.5))
    pfaf1 = (gdf_pfaf1, dict(edgecolor='black', facecolor='none', linewidth=0.8))
    title='Subbasins based on pfafstetter coding (level=1)'
    #ax = quickplot([streams, pfaf1], [pfafbas], title=title, filename='flw_pfafbas1')
# 算通過點的下游路線
# flow paths return the list of linear indices
if my:
    #xy=[[260993,2735861,'油羅上坪匯流'],[253520,2743364,'隆恩堰'],[247785,2746443,'湳雅取水口']]
    xy=([260993, 253520, 247785], [2735861, 2743364, 2746443])
    #xy=([253520], [2743364])
    
else:
    xy=([8.92, 5.55, 8.50], [50.28, 49.80, 47.3])
flowpaths, dists = flw.path(xy=xy, max_length=400e3, unit='m')
# which we than use to vectorize to geofeatures
feats = flw.geofeatures(flowpaths)
gdf_paths = gpd.GeoDataFrame.from_features(feats, crs=crs).reset_index()
gdf_pnts = gpd.GeoDataFrame(geometry=gpd.points_from_xy(*xy)).reset_index()
# and plot
if my:
    gdf_paths.to_file('output/river_c1300_path.geojson', driver='GeoJSON')
    gdf_paths.plot()
    #gdf_pnts.plot()
else:
    pnt = (gdf_pnts, dict(column='index', cmap='tab10', s=60, marker='<', zorder=4))
    fp = (gdf_paths, dict(column='index', cmap='tab10', linewidth=2))
    title='Flow path from source points (<) with max. distance of 400 km)'
    ax = quickplot([streams, fp, pnt], title=title, filename='flw_path2')<jupyter_output><empty_output><jupyter_text>有用的<jupyter_code># 河道不同細緻度
for i in range(4,12):
    feats = flw.streams(min_sto=i)
    print("min_sto=%i count %i" %(i,len(feats)))
    gdf = gpd.GeoDataFrame.from_features(feats, crs=crs)
    gdf.to_file('output/river_c1300_%i.geojson' %(i), driver='GeoJSON')

#單點到 stream 資訊: 距離，最近點，哪一個線段
from shapely.geometry import *
from shapely.ops import nearest_points

def point_with_streams(point_src): #[253520,2743364]

    point = Point(point_src[0],point_src[1])
    dist_min=5000
    idx_min=None
    for index, row in gdf.iterrows():
        line = row['geometry']
        dist = line.distance(point)
        if dist<dist_min:
            dist_min=dist
            idx_min=index
    #print("index=%i, minimal distance=%f" %(idx_min,dist_min))
    line_ori = gdf.loc[idx_min]['geometry']

    pt_in = nearest_points(line_ori, point)[0]
    #print(pt_in.coords[0][0])
    xy=pt_in.coords[0]

    return [idx_min,dist_min,xy[0],xy[1]] #index, distance, point_x, point_y

result = point_with_streams([253520,2743364])
print(result)
print("%s,%s" %(xy[0],xy[1]))
# 通過點的上游流域
points=[[260993,2735861,'油羅上坪匯流'],[253520,2743364,'隆恩堰'],[247785,2746443,'湳雅取水口']]
bas=[]
for p in points:
    x, y = np.array([p[0], p[0]+250]), np.array([p[1], p[1]+250])
    subbasins = flw.basins(xy=(x,y), streams=flw.stream_order()>=4)
    gdf_bas = vectorize(subbasins.astype(np.int32), 0, flw.transform)
    
    gdf_bas.to_file('output/river_c1300_basin_%s.geojson' %(p[2]), driver='GeoJSON')
    bas.append(gdf_bas)

#求  隆恩堰 - 油羅上坪匯流 的集水區差
a=bas[1].iloc[1]['geometry'] # iloc[0],iloc[1] 不是一樣的，要用 iloc[1] 才對
b=bas[0].iloc[1]['geometry']

c=a.difference(b)

gs = gpd.GeoSeries([c])
if 0: # 兩種方式輸出，效果相同
    gs.to_file('output/basin_diff_2point.geojson', driver='GeoJSON')
if 1:
    geojson = gs.to_json()
    with open("output/basin_diff_2point.geojson", "w") as out:
        out.write(geojson)
# 觀察 stream 的連接性
coords={}
for index, row in gdf.iterrows():
    line = row['geometry']
    idxs=row['idxs']
    points = list(line.coords)
    start = points[0]
    end = points[len(points)-1]
    coords[index]=[start,end]
    seg_cnt = len(line.coords)
    print("index=%i,length=%i,seg cnt=%i,avg_len=%.1f,start=%s,end=%s,idxs=%i" %(index,line.length,seg_cnt,line.length/seg_cnt,start,end,idxs))

link={}
for key in coords.keys():
    start,end = coords[key]
    for key2 in coords.keys():
        start2,end2 = coords[key2]
        if end==start2:
            if key in link:
                print("index=%i already have start=%i" %(key,link[key]))
            link[key]=key2
print("link=%s" %(link))

if 0:
    for i in range(len(gdf.index)):
        if i in link.keys():
            #print("index=%i PASS" %(i))
            pass
        else:
            print("index=%i,info=%s" %(index,gdf.iloc[i]))

for l in link:
    print("N%i->N%i" %(l,link[l]))
#print(gdf.iloc[22]['geometry'])
#測試用一個接近點，切開一個線段
from shapely.ops import split
from shapely.geometry import *
from shapely.ops import snap
from shapely.ops import nearest_points
if 0: #split
    pt = Point((1.15, 1.1))
    line = LineString([(0,0), (2,2)])
    result = split(line, pt)
    #result.wkt 
if 0: # snap
    square = Polygon([(1,1), (2, 1), (2, 2), (1, 2), (1, 1)])
    line = LineString([(0,0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)])
    result = snap(line, square, 0.5)
    print(result.wkt)
if 0: # 點切開線
    pt = Point((1.1, 1))
    line = LineString([(0,0), (2,2)])
    print(line)
    #print("line within point:%s,dist=%f"% (line.within(pt),line.distance(pt)))
    nearp=nearest_points(line, pt)
    pt_in = nearp[0]
    #print(pt_in)
    result = split(line, pt_in)    
    print(result.wkt)
if 1:# 圳路接入主流
    line = LineString([(0,0), (2,2)])
    line_add=LineString([(1.1,1), (1.5,0.5)])
    pt = Point(list(line_add.coords)[0])
    if 1:
        pt_in = nearest_points(line, pt)[0]
        print(pt_in)

        pts = []
        pts.append(list(pt_in.coords)[0])
        pts.extend(list(line_add.coords))
        line_add_new = LineString(pts)
        print(line_add_new)

        result = split(line, pt_in)  
        print(result.wkt)
    
type(result)
    
#圳路接入主流
# line_ori:要修改的線, line_need:想加入的線, line_append:更新起點的想加入線, line_split:要修改的線被匯流點切開的結果
#取得圳路 linestring
from shapely import wkt
wkt_str="MultiLineString ((255779.34444821099168621 2742184.59869130607694387, 255062.52472444207523949 2741882.12604631343856454, 254328.86074706495855935 2742279.99766481388360262))"
line_need = wkt.loads(wkt_str)
print("line_need=%s" %(line_need))
#經由圳路的起點，找到最接近需要修改的線
start = Point(list(line_need.geoms[0].coords)[0])
print("start:%s" %(start))
dist_min=5000
idx_min=None
for index, row in gdf.iterrows():
    line = row['geometry']
    dist = line.distance(start)
    if dist<dist_min:
        dist_min=dist
        idx_min=index
print("index=%i, minimal distance=%f" %(idx_min,dist_min))
#找到匯流點
line_ori = gdf.loc[idx_min]['geometry']
pt_in = nearest_points(line_ori, start)[0]
#修改圳路起點為匯流點
pts = []
pts.append(list(pt_in.coords)[0])
for i in range(len(line_need.geoms)):
    pts.extend(list(line_need.geoms[i].coords))
line_append = LineString(pts)
print("line_append=%s" %(line_append))
#取得主流切開後的兩線段
line_split = split(line_ori, pt_in) 
print("line_split=%s"%(line_split))
#在 gdf 中刪掉原線段，加入兩新線段,加入修改後的圳路
idx_start = len(gdf.index)

gdf.loc[idx_start]=[line_append,100,False,9] #geometry,idxs,pit,strord
for i in range(len(line_split)):
    gdf.loc[idx_start+i+1]=[line_split[i],100+i+1,False,9]
gdf.drop(idx_min,axis=0,inplace=True)
gdf.to_file('output/river_c1300_mergeline.geojson', driver='GeoJSON')<jupyter_output>line_need=MULTILINESTRING ((255779.344448211 2742184.598691306, 255062.5247244421 2741882.126046313, 254328.860747065 2742279.997664814))
start:POINT (255779.344448211 2742184.598691306)
index=11, minimal distance=0.000000
line_append=LINESTRING (255779.344448211 2742184.598691306, 255779.344448211 2742184.598691306, 255062.5247244421 2741882.126046313, 254328.860747065 2742279.997664814)
line_split=GEOMETRYCOLLECTION (LINESTRING (259119.7276749082 2738084.128264361, 259099.7253801376 2738104.130559126, 259079.7230853669 2738124.132853892, 259059.7207905963 2738124.132853892, 259039.7184958257 2738144.135148657, 259019.716201055 2738144.135148657, 258999.7139062844 2738164.137443423, 258979.7116115137 2738164.137443423, 258959.7093167431 2738164.137443423, 258939.7070219725 2738184.139738188, 258919.7047272018 2738184.139738188, 258899.7024324312 2738184.139738188, 258879.7001376605 2738184.139738188, 258859.6978428899 2738204.142032954, 258839.6955481193 2738204.142032954, 258819.6932[...]<jupyter_text>Debug<jupyter_code>#geopandas drop/insert rows test
if 0:
    line = LineString([(0,0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)])
    #gdf.loc[23]=[line,1234,False,3]
    #gdf.columns
    gdf.drop(23,axis=0,inplace=True)
gdf
type(line_split[0])
#len(gdf.index)
link
#gdf
#transform
#c=gdf.iloc[19]['geometry']
#print(c)
#c
#df_c = gpd.GeoSeries([c])
#df_c.plot()
#crs.to_epsg()
#flwdir
#elevtn
#np.savetxt('output/debug.csv', flwdir, delimiter=',')
#src2.crs

#flwdir[flwdir==-1].dtype
#flwdir[flwdir==-32767]
#flwdir
#flwdir[flwdir==-9999]
#np.savetxt('output/debug.csv', uparea[uparea!=-9999], delimiter=',')

    #pt_in = snap(pt, line, 1)
    #line_in = snap(line,pt, 0.5)
import numpy as np
import matplotlib.pyplot as plt

H = np.array([[1, 2, 3, 4],
              [5, 6, 7, 8],
              [9, 10, 11, 12],
              [13, 14, 15, 16]])  # added some commas and array creation code

fig = plt.figure(figsize=(6, 3.2))

ax = fig.add_subplot(111)
ax.set_title('colorMap')
plt.imshow(uparea)
ax.set_aspect('equal')

cax = fig.add_axes([0.12, 0.1, 0.78, 0.8])
cax.get_xaxis().set_visible(False)
cax.get_yaxis().set_visible(False)
cax.patch.set_alpha(0)
cax.set_frame_on(False)
plt.colorbar(orientation='vertical')
plt.show()
from shapely import wkt
wkt_str="MultiLineString ((255779.34444821099168621 2742184.59869130607694387, 255062.52472444207523949 2741882.12604631343856454, 254328.86074706495855935 2742279.99766481388360262))"
p1 = wkt.loads(wkt_str)
gs = gpd.GeoSeries([p1])
gs<jupyter_output><empty_output><jupyter_text>geopandas 練習<jupyter_code>import geopandas

path_to_data = geopandas.datasets.get_path("nybb")
gdf = geopandas.read_file(path_to_data)

gdf<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15981033791905971, 'fraction_numerical': 0.11312942433551489, 'mean_word_length': 4.741947852760736, 'pattern_counts': {'":': 3, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '388044', 'license': 'MIT', 'n_tokens_mistral': 7526, 'n_tokens_neox': 6634, 'n_words': 1081, 'path': 'pyflwdir.ipynb', 'repo_name': 'senkukogetsu/waterswak', 'size': '14974'}