UAV写真測量の撮影枚数(測量士試験問題)をpython で解く

# はじめに

測量士の過去問を見ていると、最近はUAV写真測量の問題が出題されています。
同じような問題を出題しているのですが、それはきっと大切な内容なのでしょう。

今回、それをプログラミングしてみました。

問題はこちら。直接ダウンロードもできます。
https://www.gsi.go.jp/LAW/SHIKEN/past.html

令和元年度、午後（記述式）のNo.２-C を解きます。



# UAV写真測量の撮影枚数をpython で計算する

## 試験問題の設定条件

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294040/6d6bd648-cdaa-b92b-6a08-f5f22ee52609.png)

設定は以下の通りです。

(1) カメラパラメータ

要素名 | 記号 | 数値 | 内容
--- | --- | --- | ---
画面距離(焦点距離) | f | 16mm | イメージセンサへの入力が像を結ぶ点
撮像素子 | s | 4um | イメージセンサの１素子の一片の大きさ
画面の大きさ（画素数) | | 6000 x 4000 | センサ素子の数

- (2) 対象区域

　要素名 | 記号 | 大きさ
  --- | --- | ---
  東西 | Y | 300 m
  南北 | X | 200 m

- (3) 撮影条件

要素名 | 記号 | 値
--- | --- | ---
地上画素寸法(撮影基準面における) | S | 1 cm
コース内の隣接写真との重複度 | $w_o$ | 80%
隣接撮影コースの空中写真との重複度 | $w_s$ | 60 %
UAV 飛行速度 | v | 4 m/s
南北両端の撮影コース | $e_1$ | 撮影される範囲のすべてが計測対象範囲の外とする
各撮影コースの両端 | $e_2$ | 撮影される範囲のすべてが計測対象範囲

問題は以下の通り。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294040/5a593656-d495-820b-a5cc-fba8b56a4ec4.png)

では解いていきましょう。

## Python による計算

### 問C-1 撮影高度

撮影基準面からの高度は以下で求まります。

$H = f\frac{S}{s}$

撮影基準面の標高が与えられているので、「海面からの撮影高度」を計算できます。

```Python
# Question C-1: flight altitude
f = 16E-3 # メートル
s = 4E-6 # 画素サイズ
gsd = 0.010 # 地上画素分解能 meter
elevation = 100 # 標高 meter
h = height_from_camera_and_gsd(f, s, gsd)
print(f"height on ground={h:.3f} [m]")
print("elevation={:.3f} [m]".format(h + elevation))
```

こたえはずばり、、

```
height on ground=40.000 [m]
elevation=140.000 [m]
```

### 問C-2 シャッター間隔

撮影間隔(基線長というらしいです）を求めます。進行方向の画像の長さを$s$とします。撮影画像のoverlap だけ進むので、これが画像間の距離になります。

$ B= S * w_0$

飛行速度$v$が与えられているので、この基線長を移動するのにかかる時間は$B/v$ です。

```Python
# Question C-2: shutter interval.
print("---")
img_size=(6000, 4000)
flight_speed=4.0 # meter/sec
overlap = 0.80

def image_baseline(overlap: float, gnd_img_len: float):
    return (1-overlap) * gnd_img_len

B = image_baseline(overlap, img_size[1] * gsd)
print("ground image size: ({:.2f}, {:.2f}) [m]".format(img_size[0]*gsd, img_size[1]*gsd))
print("baseline: {:.2f} [m]".format(B))
print("shutter interval {:.2f} [sec]".format(B / flight_speed))
```

答えはずばり、

```
ground image size: (60.00, 40.00) [m]
baseline: 8.00 [m]
shutter interval 2.00 [sec]
```

### 問C-3 コース数

ずばりコース数をn として以下を解きます。

$ Y + L\cdot e_1 \cdot 2 < L_w [1 + (n-1)(1-w_c)]$

これを満たす整数ｎを求める計算を、そのまま実装します。

```Python
land_size = (200, 300)
cource_edge_overlap = 1.0
sidelap = 0.60

from numpy import ceil
def course_number(land_len: float, img_len: float, sidelap: float, edge_overlap: float) -> int:
    min_n = (land_len/img_len + 2 * edge_overlap - sidelap) / (1.0 - sidelap)
    return int(ceil(min_n))

N1 = course_number(land_size[0], img_size[0]*gsd, sidelap, cource_edge_overlap)
print(f"{N1} cources")
```

こたえはずばり、、、

```
12 cources
```

### 問C-4 撮影枚数

同じように、１コース内の撮影枚数を求めます。コース数はさきほど求めたので合計枚数を計算できます。

```Python
overlap = 0.80

def photo_in_cource(land_len: float, img_len: float, overlap, edge_overlap: float) -> int:
    return course_number(land_len, img_len, overlap, edge_overlap)

N2 = photo_in_cource(land_size[1], img_size[1]*gsd, overlap, cource_edge_overlap)
print(f"{N2} photo in a cource")
print("total: {} photos".format(N1*N2))
```

答えはずばり

```
44 photo in a cource
total: 528 photos
```

答え合わせ。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294040/97a6e35a-abf9-57a5-ca21-e76564736fc9.png)


# まとめ

UAV写真測量の写真枚数の計算を、測量士の試験問題を例にpython で実装した。答え合わせをして計算結果が正しいことを確認できた。

今後、改良するとしたら、、、
- class にする。単純なrefactoring
- 図で表示したい
- ブラウザで動作するようにしたい

とは思うが、今日はこのあたりで。疲れた。測量士の試験の準備、中途半端だったなー(^^;
(2021/09/14)