定期横断測量をQGISへ_対象flagの検出とgeopkgeの作成

## １．定期横断測量のDL、データ確認、距離標データの作成

仕事の関係で、堤外地の堤防法尻点の座標点が必要となりました。
定期横断測線に座標点を与えるのは、少々手間なので、その作業をpythonで実施しています。

北海道開発局が[定期横断測量データ](https://www.hkd.mlit.go.jp/ky/kn/kawa_kan/slo5pa000000zpx4.html)を公開しています。

データはzip形式で、以下の３つのファイルが河川ごとに格納されています。
距離標データ（数値データ）
縦断データ（整理表、数値データ）
横断データ（整理表、数値データ）

### 定期横断測量のDL
ここでは、以下のデータを利用します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3956606/16aabe51-32be-dee8-f677-e3539856633e.png)

### 各種データの確認
#### 距離標データ
距離標データは、杭の座標点が格納されてます。
既設（古い）の距離標の座標点と移設の距離標の座標点が2つ載っており、測量年に併せてデータを加工する必要があります。データを利用するよりも表としての見やすさを優先されているようで、改善が必要な箇所ですね。

#### 距離標データの作成
使用する横断データの計測年に併せて、既設か移設のKP座標のファイルを別途作成します。

<details><summary>サンプルcsv</summary>

```
以下をyuubetsu.csvで保存してください。

KP,LX,LY,RX,RY,sokuten
0,"26,700.580 ","109,456.805 ","26,608.004 ","109,526.703 ",116.000 
0.2,"26,717.296 ","109,368.164 ","26,550.595 ","109,358.167 ",169.000 
0.4,"26,742.624 ","109,307.594 ","26,557.055 ","109,233.001 ",200.000 
0.5,"26,786.367 ","109,190.553 ","26,574.120 ","109,127.200 ",221.500 
0.6,"26,853.991 ","109,026.201 ","26,587.064 ","109,001.539 ",228.064 
0.8,"26,894.888 ","108,887.798 ","26,582.502 ","108,847.016 ",275.036 
1,"26,865.157 ","108,518.789 ","26,525.207 ","108,749.049 ",410.592 
1.2,"26,755.704 ","108,384.165 ","26,409.101 ","108,619.067 ",418.703 
1.4,"26,617.376 ","108,212.772 ","26,308.052 ","108,478.444 ",407.754 
1.6,"26,472.648 ","108,046.080 ","26,179.328 ","108,330.500 ",408.572 
1.8,"26,294.429 ","107,885.497 ","26,027.567 ","108,192.131 ",406.497 
2,"26,107.890 ","107,748.564 ","25,901.169 ","108,099.636 ",407.413 
2.2,"25,945.519 ","107,652.536 ","25,788.879 ","108,031.375 ",409.945 
2.4,"25,786.529 ","107,574.742 ","25,640.963 ","107,955.991 ",408.093 
2.6,"25,585.687 ","107,498.709 ","25,470.934 ","107,889.390 ",407.185 
2.8,"25,423.886 ","107,453.027 ","25,338.836 ","107,852.466 ",408.393 
3,"25,273.554 ","107,423.833 ","25,202.898 ","107,825.494 ",407.828 
3.2,"25,078.097 ","107,403.076 ","25,045.821 ","107,808.646 ",406.852 
3.4,"24,891.559 ","107,397.103 ","24,891.001 ","107,805.154 ",408.051 
3.6,"24,712.446 ","107,394.072 ","24,706.799 ","107,801.646 ",407.613 
3.8,"24,557.284 ","107,385.711 ","24,512.046 ","107,791.884 ",408.684 
4,"24,399.582 ","107,369.810 ","24,325.354 ","107,774.303 ",411.247 
4.2,"24,236.140 ","107,347.859 ","24,169.278 ","107,757.436 ",414.999 
4.4,"24,002.847 ","107,302.477 ","23,914.929 ","107,701.190 ",408.291 
4.6,"23,747.054 ","107,234.104 ","23,698.226 ","107,661.698 ",430.373 
4.8,"23,627.863 ","107,197.333 ","23,456.893 ","107,604.406 ",441.519 
5,"23,493.592 ","107,148.412 ","23,342.287 ","107,496.470 ",379.523 
5.2,"23,320.131 ","107,077.300 ","23,170.845 ","107,434.461 ",387.105 
5.4,"23,123.207 ","106,986.404 ","22,973.261 ","107,356.249 ",399.085 
5.6,"22,958.679 ","106,909.496 ","22,791.587 ","107,280.481 ",406.878 
5.8,"22,797.868 ","106,830.734 ","22,623.395 ","107,204.992 ",412.928 
6,"22,637.958 ","106,707.542 ","22,441.505 ","107,118.430 ",455.437 
6.2,"22,506.037 ","106,650.711 ","22,233.096 ","107,073.793 ",503.483 
6.4,"22,259.335 ","106,591.929 ","22,094.142 ","106,947.845 ",392.384 
6.6,"22,149.378 ","106,324.056 ","21,925.154 ","106,854.571 ",575.954 
6.8,"22,009.118 ","106,227.354 ","21,712.441 ","106,720.380 ",575.406 
7,"21,796.726 ","106,224.353 ","21,546.586 ","106,602.767 ",453.616 
7.2,"21,613.057 ","106,069.973 ","21,320.838 ","106,421.863 ",457.404 
7.4,"21,415.835 ","105,984.201 ","21,175.238 ","106,297.228 ",394.807 
7.6,"21,404.554 ","105,746.338 ","21,041.919 ","106,175.584 ",561.922 
7.8,"21,315.410 ","105,571.696 ","20,887.152 ","106,030.280 ",627.459 
8,"21,231.918 ","105,479.250 ","20,742.137 ","105,893.606 ",641.542 
8.2,"20,995.135 ","105,343.031 ","20,581.785 ","105,741.753 ",574.315 
8.4,"20,830.425 ","105,237.474 ","20,442.428 ","105,610.571 ",538.278 
8.6,"20,758.495 ","105,092.449 ","20,344.845 ","105,522.262 ",596.528 
8.8,"20,742.887 ","104,917.205 ","20,264.363 ","105,451.552 ",717.295 
9,"20,535.580 ","104,820.529 ","20,177.317 ","105,376.800 ",661.657 
9.2,"20,390.861 ","104,824.311 ","20,084.403 ","105,298.962 ",564.987 
9.4,"20,221.921 ","104,802.800 ","19,961.709 ","105,202.397 ",476.852 
9.6,"20,041.100 ","104,717.728 ","19,798.180 ","105,081.061 ",437.059 
9.8,"19,847.963 ","104,626.810 ","19,635.338 ","104,965.922 ",400.258 
10,"19,718.491 ","104,549.907 ","19,450.980 ","104,842.500 ",395.450 
10.2,"19,552.021 ","104,443.731 ","19,330.170 ","104,765.717 ",391.015 
10.4,"19,400.488 ","104,348.027 ","19,195.808 ","104,683.644 ",393.106 
10.6,"19,216.422 ","104,231.848 ","19,004.591 ","104,576.074 ",404.183 
10.8,"19,033.810 ","104,139.015 ","18,846.822 ","104,491.251 ",397.792 
11,"18,877.987 ","104,061.177 ","18,680.959 ","104,403.841 ",394.271 
11.2,"18,687.581 ","103,966.383 ","18,514.060 ","104,322.114 ",394.795 
11.4,"18,513.194 ","103,867.593 ","18,363.019 ","104,256.742 ",416.120 
11.6,"18,373.295 ","103,722.518 ","18,191.681 ","104,182.079 ",494.146 
11.8,"18,265.417 ","103,591.817 ","18,049.356 ","104,121.992 ",572.510 
12,"18,095.572 ","103,480.378 ","17,883.848 ","104,058.444 ",615.619 
12.2,"17,915.961 ","103,409.663 ","17,709.257 ","104,004.598 ",629.821 
12.4,"17,739.037 ","103,341.059 ","17,532.011 ","103,947.845 ",641.131 
12.6,"17,565.904 ","103,277.816 ","17,396.872 ","103,905.310 ",649.862 
12.8,"17,339.372 ","103,229.009 ","17,222.260 ","103,856.360 ",638.188 
13,"17,162.031 ","103,191.543 ","17,033.838 ","103,816.383 ",637.855 
13.2,"16,963.918 ","103,100.202 ","16,846.039 ","103,771.712 ",681.778 
13.4,"16,760.075 ","103,025.724 ","16,667.497 ","103,743.593 ",723.814 
13.6,"16,550.625 ","103,012.941 ","16,481.406 ","103,717.122 ",707.575 
13.8,"16,288.814 ","103,158.900 ","16,285.273 ","103,700.095 ",541.207 
14,"16,128.708 ","103,245.769 ","16,079.625 ","103,686.200 ",443.158 
14.2,"15,902.368 ","103,249.559 ","15,867.918 ","103,672.030 ",423.873 
14.4,"15,731.328 ","103,196.838 ","15,638.184 ","103,656.173 ",469.084 
14.6,"15,498.684 ","103,059.806 ","15,432.692 ","103,642.531 ",586.450 
15,"15,214.686 ","103,041.076 ","15,188.890 ","103,625.049 ",584.542 
15.2,"15,071.395 ","103,070.200 ","14,979.523 ","103,621.896 ",559.293 
15.4,"14,763.596 ","103,194.506 ","14,705.963 ","103,592.970 ",403.010 
15.6,"14,483.551 ","103,173.920 ","14,508.344 ","103,578.675 ",405.514 
15.8,"14,253.487 ","103,156.999 ","14,333.523 ","103,566.683 ",417.429 
16,"14,121.997 ","103,126.070 ","14,158.450 ","103,567.868 ",443.699 
16.2,"13,902.504 ","103,124.478 ","14,016.304 ","103,594.630 ",483.729 
16.4,"13,721.649 ","103,147.075 ","13,839.226 ","103,630.651 ",497.965 
16.6,"13,440.581 ","103,189.910 ","13,679.346 ","103,677.932 ",526.299 
17,"13,281.440 ","103,372.115 ","13,534.654 ","103,731.312 ",439.477 
17.2,"13,080.452 ","103,560.488 ","13,379.338 ","103,855.559 ",420.000 
17.4,"12,969.605 ","103,626.601 ","13,286.952 ","103,978.556 ",468.400 
17.6,"12,817.147 ","103,678.650 ","12,853.116 ","104,236.611 ",552.619 
18,"12,625.985 ","103,689.956 ","12,579.640 ","104,175.683 ",487.933 
18.2,"12,428.074 ","103,669.603 ","12,385.042 ","104,163.687 ",495.954 
18.4,"12,243.636 ","103,651.909 ","12,227.727 ","103,939.608 ",288.139 
18.6,"12,095.462 ","103,619.012 ","12,063.241 ","103,920.210 ",298.917 
18.8,"11,923.217 ","103,580.169 ","11,880.158 ","103,881.487 ",304.379 
19,"11,799.292 ","103,516.556 ","11,651.349 ","103,753.646 ",279.462 
19.2,"11,706.492 ","103,397.329 ","11,434.578 ","103,595.624 ",336.538 
19.4,"11,601.620 ","103,241.414 ","11,328.268 ","103,455.749 ",347.363 
19.6,"11,360.357 ","103,068.631 ","11,155.286 ","103,357.336 ",354.125 
19.8,"11,137.248 ","102,974.480 ","11,003.251 ","103,302.192 ",355.049 
20,"10,991.852 ","102,912.539 ","10,954.860 ","103,261.016 ",351.435 
20.2,"10,803.846 ","102,832.573 ","10,691.133 ","103,214.451 ",398.165 
20.4,"10,639.843 ","102,765.367 ","10,565.593 ","103,153.905 ",394.569 
20.6,"10,500.617 ","102,704.629 ","10,382.594 ","103,107.755 ",420.048 
20.8,"10,294.794 ","102,661.468 ","10,251.939 ","103,069.315 ",410.092 
21,"10,051.581 ","102,644.188 ","10,128.298 ","103,119.968 ",481.925 
21.2,"9,901.248 ","102,650.157 ","9,923.507 ","103,200.849 ",551.142 
21.4,"9,657.805 ","102,683.107 ","9,765.980 ","103,204.002 ",532.009 
21.7,"9,447.952 ","102,713.037 ","9,494.084 ","103,247.760 ",536.709 
22,"9,220.075 ","102,763.156 ","9,263.475 ","103,105.398 ",344.983 
22.2,"8,982.377 ","102,820.375 ","9,065.935 ","103,151.833 ",341.828 
22.4,"8,791.737 ","102,880.988 ","8,913.913 ","103,114.200 ",263.277 
22.6,"8,662.858 ","102,921.847 ","8,735.274 ","103,160.428 ",249.329 
22.8,"8,467.144 ","102,984.240 ","8,551.817 ","103,267.801 ",295.933 
23,"8,257.983 ","103,051.101 ","8,505.357 ","103,416.554 ",441.305 
23.2,"8,172.252 ","103,076.370 ","8,171.204 ","103,496.890 ",420.521 
23.4,"8,048.327 ","103,072.080 ","7,748.145 ","103,586.007 ",595.172 
23.7,"7,870.657 ","102,981.604 ","7,619.957 ","103,342.093 ",439.093 
24,"7,698.027 ","102,880.022 ","7,494.336 ","103,247.692 ",420.323 
24.2,"7,548.238 ","102,786.894 ","7,223.354 ","103,044.228 ",414.452 
24.4,"7,303.463 ","102,621.495 ","7,131.803 ","102,975.363 ",393.306 
24.6,"7,117.375 ","102,542.931 ","6,980.587 ","102,859.716 ",345.056 

```
</details>


#### 横断データ
WZAAA4001.CSVの中身をみると、数値データが羅列されています。
1行目は、以下に示す基本的な内容が記載されてます。
なお、堤防法尻は、堤内地側の法尻を示しています。
詳細は、[河川定期縦横断データ作成ガイドライン（H20)](
https://www.mlit.go.jp/river/shishin_guideline/kasen/gis/pdf_docs/juoudan/guideline0805.pdf) を参考にしてください。

| 距離標 | 流心間距離 | 左岸距離杭高 | 右岸距離杭高 |  左岸水際杭高 | 右岸水際杭高 |  測点数 | 左岸堤防法尻高 | 左岸堤内地平均 | 右岸堤防法尻高 | 右岸堤内地平均 | 河川 構造フラグ | 測量年月日 | 河川番号 | 水系名      | 河川名称 |
|--------|------------|--------------|--------------|---------------|--------------|---------|----------------|----------------|----------------|----------------|-----------------|------------|----------|-------------|----------|
| 2.2    | 170        | 7.249        | 6.732        | 3.608         | 3.94         | 110     | 3.54           | 2.55           | 3.69           | 3.48           | 0               | 20190122   | 8.1E+09  | 湧別川水系  | 湧別川   |

測線データは2行目以下になり、以下の順に並んでいます。
標高は平面直角座標系（12系）での値を示しています。

| flag | 左岸杭からの距離 | 標高 |
|------|------------------|------|
| 4    | -28.5            | 2.37 |
| 4    | -21.53           | 2.4  |
| 28   | -14.99           | 2.55 |

### flagの意味
| 記号 | 座標点の特性                         |
|------|--------------------------------------|
| 1    | 距離杭高（左岸）                     |
| 2    | 距離杭下                             |
| 3    | 堤防法肩位置（始点，終点）           |
| 4    | 道路（始点，終点）                   |
| 5    | 水路（始点，終点）                   |
| 6    | 河川管理境界                         |
| 7    | 官民界                               |
| 8    | コンクリート張（始点，終点）         |
| 9    | 石張（始点， 終点）                  |
| 10   | コンクリートブロック張（始点，終点） |
| 11   | じゃかご張（始点，終点）             |
| 12   | 水際杭高（左岸，右岸）               |
| 13   | 低水路肩（右岸，左岸）               |
| 14   | アスファルト（始点，終点）           |
| 15   | 水制（始点，終点）                   |
| 16   | 水面（右岸，左岸）                   |
| 17   | 根固ブロック（始点，終点）           |
| 18   | 距離杭高（右岸）                     |
| 19   | 鉄柵（地盤高，鉄柵の上端）           |
| 20   | 生垣（地盤高，生垣の上端）           |
| 21   | ブロック塀（始点，終点）             |
| 22   | ガードレール（地盤高，生垣の上端）   |
| 23   | 家屋（地盤高）                       |
| 24   | 左岸本堤防法尻（表・裏2点）          |
| 25   | 右岸本堤防法尻（表・裏2点）          |
| 26   | その他の堤防法尻                     |
| 27   | ガンバン（始点，終点）               |
| 28   | 堤内地盤（始点，終点）               |

この表は、[Table Generator](https://www.tablesgenerator.com/markdown_tables)を用いて作成されてます。

## ２．座標変換からgeopkageの作成
　横断側線上の各点に座標を与える処理をします。コードは以下のサンプルコードを参照してください。
　左岸杭はleft_point、同様に左岸の堤外法尻はleft_toe_pointとしています。
 
| kp  | flag | distance | elevation | x        | y        | type            |
|-----|------|----------|-----------|----------|----------|-----------------|
| 2.2 | 4    | -28.5    | 2.37      | 25956.41 | 107626.2 |                 |
| 2.2 | 1    | 0        | 7.249     | 25945.52 | 107652.5 | left_point      |
| 2.2 | 2    | 0        | 7.06      | 25945.52 | 107652.5 |                 |
| 2.2 | 24   | 14.02    | 3.66      | 25940.16 | 107665.5 | left_toe_point |　
<details><summary>座標変換</summary>

```python
#zahyou_henkan.pyとして保存
import math
import pandas as pd
import numpy as np
import os
import glob

def process_files(coordinate_file_path, input_pattern, output_dir):
    """
    ファイルを処理する主要な関数
    
    Args:
    - coordinate_file_path: 座標ファイル（yuubetsu.csv）のパス
    - input_pattern: 入力ファイルのパターン（出力フォルダ内のW*.CSV）
    - output_dir: 出力フォルダのパス
    """
    try:
        # 1. yuubetsu.csvを読み込む
        df_coordinates = pd.read_csv(coordinate_file_path, encoding='shift_jis', thousands=',')
        
        # 座標列の型変換
        coordinate_columns = ['KP', 'LX', 'LY', 'RX', 'RY', 'sokuten']
        for col in coordinate_columns:
            df_coordinates[col] = pd.to_numeric(df_coordinates[col], errors='coerce')
        
        print("座標データの最初の行:")
        print(df_coordinates.head(1))
    except Exception as e:
        print(f"座標ファイル読み込みエラー: {str(e)}")
        return
    
    # 2. outputフォルダのcsvファイルを読み込む（W始まりのファイル）
    input_files = glob.glob(input_pattern)

    for file in input_files:
        try:
            # 出力ファイルのキロポストを読み取る
            result_df = pd.read_csv(file, encoding='shift_jis', skiprows=1, header=None)
            result_df = result_df.iloc[:, :3]
            result_df.columns = ['flag', 'distance', 'elevation']
            
            # ファイルのキロポストを取得
            kilopost = pd.read_csv(file, encoding='shift_jis', nrows=1, header=None).iloc[0, 0]
            print(f"\n処理中のファイル: {file}")
            print(f"ファイルのキロポスト: {kilopost}")
            
            # 3. yuubetsu.csvのKPと完全一致するかを確認
            matching_coords = df_coordinates[df_coordinates['KP'] == float(kilopost)]
            
            if matching_coords.empty:
                print(f"警告: キロポスト {kilopost} に一致する座標が見つかりません。スキップします。")
                continue
            
            # 4. マッチングした座標を取得
            matching_row = matching_coords.iloc[0]
            
            # flag=1の箇所にLX, LYを転記
            left_point_mask = result_df['flag'] == 1
            if left_point_mask.any():
                left_point_index = left_point_mask[left_point_mask].index[0]
                result_df.loc[left_point_index, 'x'] = matching_row['LX']
                result_df.loc[left_point_index, 'y'] = matching_row['LY']
                result_df.loc[left_point_index, 'type'] = 'left_point'
            
            # flag=18の箇所にRX, RYを転記
            right_point_mask = result_df['flag'] == 18
            if right_point_mask.any():
                right_point_index = right_point_mask[right_point_mask].index[0]
                total_line_length = result_df.loc[right_point_index, 'distance']
            else:
                # flag=18がない場合、RX, RYを新しい行として追加
                new_row = pd.DataFrame({
                    'flag': [18],
                    'distance': [matching_row['sokuten']],
                    'elevation': [0],  # デフォルト値
                    'x': [matching_row['RX']],
                    'y': [matching_row['RY']],
                    'type': ['right_point_added']
                })
                result_df = pd.concat([result_df, new_row], ignore_index=True)
                total_line_length = matching_row['sokuten']
                print(f"注意: キロポスト {kilopost} にflag=18の点がないため、sokuten値を使用して追加しました。")
            
            # 左右の杭座標
            left_bank_x, left_bank_y = matching_row['LX'], matching_row['LY']
            right_bank_x, right_bank_y = matching_row['RX'], matching_row['RY']
            
            # 線分の方向ベクトルを計算
            line_dx = right_bank_x - left_bank_x
            line_dy = right_bank_y - left_bank_y
            
            # まだx,yが設定されていない行の座標を計算
            mask_without_coords = result_df['x'].isna()
            result_df.loc[mask_without_coords, 'x'] = (
                left_bank_x + 
                (line_dx / total_line_length) * result_df.loc[mask_without_coords, 'distance']
            )
            result_df.loc[mask_without_coords, 'y'] = (
                left_bank_y + 
                (line_dy / total_line_length) * result_df.loc[mask_without_coords, 'distance']
            )
            
            # KPの列を追加
            result_df.insert(0, 'kp', kilopost)
            
            # 6. 法尻高情報の追加, 法尻高の抽出
            result_df['left_toe_point'] = None
            result_df['right_toe_point'] = None
            
            # 法尻高の抽出
            left_toe_flags = result_df[result_df['flag'] == 24]
            if len(left_toe_flags) >= 2:
                left_toe_point = left_toe_flags.iloc[1]['elevation']
                result_df.loc[left_toe_flags.index[1], 'type'] = 'left_toe_point'
                result_df.loc[left_toe_flags.index[1], 'left_toe_point'] = left_toe_point
            
            right_toe_flags = result_df[result_df['flag'] == 25]
            if len(right_toe_flags) >= 1:
                right_toe_point = right_toe_flags.iloc[0]['elevation']
                result_df.loc[right_toe_flags.index[0], 'type'] = 'right_toe_point'
                result_df.loc[right_toe_flags.index[0], 'right_toe_point'] = right_toe_point
            
            # 丸め処理
            columns_to_round = ['x', 'y', 'distance', 'elevation']
            result_df[columns_to_round] = result_df[columns_to_round].round(6)
            
            # 結果の保存
            base_name = os.path.basename(file)
            output_filename = os.path.join(output_dir, 
                os.path.splitext(base_name)[0] + '_points.csv')
            
            result_df.to_csv(output_filename, index=False)
            print(f"処理完了: {base_name} → {os.path.basename(output_filename)}")
            
        except Exception as e:
            print(f"エラー: {os.path.basename(file)} の処理中にエラーが発生しました: {str(e)}")

def main():
    # スクリプトがあるディレクトリを取得
    current_dir = os.path.dirname(os.path.abspath(__file__))
    
    # 設定パラメータ
    input_dir = os.path.join(current_dir, 'oudanseika_03')
    coordinate_file_path = os.path.join(current_dir, 'yuubetsu.csv')
    input_pattern = os.path.join(input_dir, 'W*.CSV')
    output_dir = os.path.join(current_dir, 'output')
    
    # 出力ディレクトリ作成
    os.makedirs(output_dir, exist_ok=True)
    print(f"出力ディレクトリ: {output_dir}")
    
    # ファイル処理の実行
    process_files(coordinate_file_path, input_pattern, output_dir)
    
    print("\n注意: このデータの座標系は JGD2000 / 平面直角座標系 12系 (EPSG:2451) です")

if __name__ == "__main__":
    main()
```
</details>

<details><summary>geopkgeへの変換</summary>

```python
#geopkge_hankan.pyで保存

import os
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

def csv_to_geopackage(csv_path, output_path, epsg_code):
    # CSVファイルを読み込む
    df = pd.read_csv(csv_path)

    # 必要な列だけを選択し、データ型を最適化する
    df = df[['x', 'y', 'type', 'left_toe_point', 'right_toe_point']].copy()
    
    # Toe heightが設定されているかどうかを確認する列を追加
    df['marker'] = df.apply(lambda row: 'toe_point' if row.get('left_toe_point', 0) > 0 or row.get('right_toe_point', 0) > 0 else '', axis=1)
    
    # GeoDataFrameを作成するためにポイントを作成
    geometry = [Point(xy) for xy in zip(df['y'], df['x'])]
    gdf = gpd.GeoDataFrame(df, geometry=geometry)

    # CRSを設定（引数で指定されたEPSGコードを使用）
    gdf.set_crs(epsg=epsg_code, inplace=True)

    # GeoPackageに保存
    gdf.to_file(output_path, driver='GPKG')

    print(f"Converted {csv_path} to {output_path}")

def main():
    # ファイルの箇所を指定
    input_dir = r'E:\python_zahyou\output'
    output_dir = os.path.join(input_dir, 'geopackage_output')
    os.makedirs(output_dir, exist_ok=True)

    # EPSGコードを指定
    epsg_code = 2454  # JGD2000の場合

    # ファイル数をカウント
    csv_files = [f for f in os.listdir(input_dir) if f.endswith('.csv')]
    total_files = len(csv_files)
    
    for i, filename in enumerate(csv_files, start=1):
        input_path = os.path.join(input_dir, filename)
        output_path = os.path.join(output_dir, f'{os.path.splitext(filename)[0]}.gpkg')
        try:
            csv_to_geopackage(input_path, output_path, epsg_code)
            print(f"Progress: {i}/{total_files} files processed.")
        except Exception as e:
            print(f"Error processing {filename}: {e}")

if __name__ == '__main__':
    main()

```
</details>

### 任意のflag（堤外地の法尻抽出）データを抽出する
　法尻点のみ抽出したpkgeとcsvファイルの作成をします。

<details><summary>法尻点の抽出</summary>

```python
#pont.pyとして保存
import os
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

def csv_to_geopackage(csv_path, output_path, epsg_code):
    # CSVファイルを読み込む
    df = pd.read_csv(csv_path)

    # 必要な列だけを選択し、データ型を最適化する
    df = df[['x', 'y', 'type', 'left_toe_point', 'right_toe_point']].copy()
    
    # Toe heightが設定されているかどうかを確認する列を追加
    df['marker'] = df.apply(lambda row: 'toe_point' if row.get('left_toe_point', 0) > 0 or row.get('right_toe_point', 0) > 0 else '', axis=1)
    
    # GeoDataFrameを作成するためにポイントを作成
    geometry = [Point(xy) for xy in zip(df['y'], df['x'])]
    gdf = gpd.GeoDataFrame(df, geometry=geometry)

    # CRSを設定（引数で指定されたEPSGコードを使用）
    gdf.set_crs(epsg=epsg_code, inplace=True)

    # GeoPackageに保存
    gdf.to_file(output_path, driver='GPKG')

    print(f"Converted {csv_path} to {output_path}")

def main():
    # ファイルの箇所を指定
    input_dir = r'E:\python_zahyou\output'
    output_dir = os.path.join(input_dir, 'geopackage_output')
    os.makedirs(output_dir, exist_ok=True)

    # EPSGコードを指定
    epsg_code = 2454  # JGD2000の場合

    # ファイル数をカウント
    csv_files = [f for f in os.listdir(input_dir) if f.endswith('.csv')]
    total_files = len(csv_files)
    
    for i, filename in enumerate(csv_files, start=1):
        input_path = os.path.join(input_dir, filename)
        output_path = os.path.join(output_dir, f'{os.path.splitext(filename)[0]}.gpkg')
        try:
            csv_to_geopackage(input_path, output_path, epsg_code)
            print(f"Progress: {i}/{total_files} files processed.")
        except Exception as e:
            print(f"Error processing {filename}: {e}")

if __name__ == '__main__':
    main()
```
</details>

<details><summary>一括処理(おまけ)</summary>

```python
'''
 上記のすべてのファイルを一括で実行するコードです。
 各ファイル中のmain()で設定されているdirectory等も
 ここで修正すると便利だと思います。
'''
import os
import sys
import subprocess

def run_script(script_path):
    """
    スクリプトを実行し、エラーがあれば表示する関数
    """
    try:
        print(f"\n----- {os.path.basename(script_path)} の実行を開始 -----")
        result = subprocess.run([sys.executable, script_path], 
                                capture_output=True, 
                                text=True, 
                                check=True)
        print(result.stdout)
        print(f"----- {os.path.basename(script_path)} の実行が完了 -----\n")
    except subprocess.CalledProcessError as e:
        print(f"エラー: {script_path} の実行中にエラーが発生しました")
        print("標準出力:", e.stdout)
        print("標準エラー:", e.stderr)

def main():
    # スクリプトの相対パスまたは絶対パスを指定
    scripts = [
        'zahyou_henkan.py',     # 座標変換スクリプト
        'geopkge_henkan.py',    # GeoPackage変換スクリプト
        'point.py'               # ポイント処理スクリプト
    ]

    # カレントディレクトリを設定（このmain.pyがある場所）
    script_dir = os.path.dirname(os.path.abspath(__file__))
    
    # 各スクリプトを順番に実行
    for script in scripts:
        script_path = os.path.join(script_dir, script)
        if os.path.exists(script_path):
            run_script(script_path)
        else:
            print(f"警告: スクリプト {script} が見つかりません")

if __name__ == "__main__":
    main()
```
</details>

## ３．QGISへ図化する
### 背景画像の表示
[国土地理院QGISマニュアル](https://nlftp.mlit.go.jp/ksj/other/QGIS_manual.pdf)を参考に、背景に写真を追記します。
　XYZ Tiles－新規接続を選び、接続の詳細にURL、最小・最大ズームを入れてください。
 
### 測点を表示
geopacage_outputのファイルをQGIS上にドラック＆ドロップすると側線のポイントデータが表示されます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3956606/98950bba-347d-2a63-32b6-ed1888c1c4d2.png)

### 任意のflag点（法尻のpointデータ）を表示
任意のflagのデータについても同様に表示できます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3956606/607a5538-4f23-220c-286d-75596e9e6793.png)