Pythonで世界地図-20(地図画像に位置情報を付与する)

mpl_toolkits.basemapのサンプルプログラムを探していたら、海外のサイトで*.tifをロードするプログラムがよくあります。
プログラム自体は動作するのですが、tifがよくわからないのでサンプルのtifを開こうとしても開くことが出来ません？？
昔はtifファイルを開いたことがあるのですが？？
（tifは開くことができました、tiffが開けなかった。）

調べていたら、GEOTIFFと言う位置情報が付与されたtifファイルのようです。

QGISと言うアプリで作成出来るようなのでインストールしました。

#####インストール
Ubuntu Studio 18.04の場合
qgis/bionic 2.18.17+dfsg-1 amd64  地理情報システム (GIS)

```shell-session
$ sudo apt install qgis
```

最新バージョンをインストールする場合
Windows・Mac・Linux(自分の環境にあったQGISのダウンロード)
最新バージョン QGIS 3.2.0 'Bonn'  22.06.2018
https://qgis.org/ja/site/forusers/download.html

#####地図画像を作成
官公庁のサイトからダウンロード、またはGoogleMap等で作成します。
私は、mpl_toolkits.basemapで作成しました。（これも、ポイントを決める時にズームすると、緯度・経度線が太くてなかなか精密にポイントを決めるのが難しいです。）

一度、ネットにあった精度の悪そうな日本全土の地図でやって見ましたが、実際の位置と100km単位でずれて、使い物になりませんでした。

![geotiff.jpg](https://qiita-image-store.s3.amazonaws.com/0/238800/d8e58bde-9092-c666-3d33-c327c5b98d84.jpeg)

```python:
#!/usr/bin/python3
# coding: UTF-8

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(10,8))
map = Basemap(projection='cyl',
                 resolution='f',
                 llcrnrlon=123,
                 llcrnrlat=24,
                 urcrnrlon=147,
                 urcrnrlat=46)

map.drawcoastlines(color='lightgray')
map.drawcountries(color='lightgray')
map.fillcontinents(color='white', lake_color='#eeeeee');
map.drawmapboundary(fill_color='#eeeeee')
map.drawparallels(np.arange(24.0, 46, 1.0), labels = [1,0,0,0], fontsize=12)
map.drawmeridians(np.arange(123.0, 147, 2.0), labels = [0,0,0,1], fontsize=12)

plt.show()
```

#####地図画像に位置情報を付与する。
参考
ジオレファレンサプラグイン
https://docs.qgis.org/2.18/ja/docs/user_manual/plugins/plugins_georeferencer.html

１．地図画像を作成する。
GoogleMap等の地図を利用する時は、目印になる道路と道路の交差点等を、右クリックの「この場所について」等で東経・北緯を先に調べておく。（４箇所）

mpl_toolkits.basemapで作成する時は、東経・北緯線を表示しておく。
地図の図法は何でも良いわけではなく、GISに向かないものがあるそうです。
GoogleMap等はwebメルカトル図法らしいので、webメルカトル図法はなさそうなので今回はメルカトル図法projection='merc'で作成しました。（projection='merc'　で作成すると非常に誤差が大きいです。）
####注意
必ず、projection='cyl'で作成する必要があります。
'cyl'で作成するとほとんど誤差がありませんでした。

２．QGIS　Desktopを起動する。

地図の保存を行う。
プロジェクト（J)　ー＞　名前を付けて保存(A)　ー＞　フォルダを指定して名前を入力する。
****.qgs

３．位置情報を付加するにはジオリファレンサーを有効にする必要がある。
プラグイン(P)　ー＞　プラグインの管理とインストール...　ー＞　「GDALジオリファレンサー」にチェックを入れ、閉じる。
４．ラスタ(R)　ー＞　ジオリファレンサー(G)　ー＞　ジオリファレンサー(G)　ー＞　 ジオリファレンサーのウインドウが表示される。
ファイル　ー＞　ラスタを開く...　ー＞　位置情報を付与する画像を開く。
空間参照システム選択　ダイアログが開く。
WGS84を選択してOKをクリックする。　ー＞　地図が表示される。

５．地図に位置情報を付与する。
ジオリファレンサー　ウインドウの編集　ー＞　ポイントの追加　ー＞　ポイントをクリックすると、地図座標を入力　ダイアログが表示されるので、クリックした地点の東経と北緯を入力しOKをクリックする。　ー＞　GCPテーブルに変換元と変換先の座標が表示される。

注意）
この作業は、正確に行わないと地図がずれてしまいます。
特に上記の地図のように範囲の大きな地図の場合は、特にずれてしまいます。
地図をズームして正確にポイントを決めなくてはいけません。

この作業を4箇所行う。（間違っていたら、ダブルクリックで修正する。）

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/62302c8d-dca2-716d-fdf0-7482c82d7546.png)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/b4099010-0fbf-5886-3a13-84a2b7470199.png)

６．左上の右三角形の「ジオリファレンスの開始」をクリックする。　ー＞　情報　ダイアログが表示される。　ー＞　　
　OK　ー＞　変換の設定　ダイアログが表示される。　ー＞　
変換タイプ　：　多項式１
リサンプリング方法　：　最近傍
変換先SRS　：　デフォルトCRS(EPSG:4326-WGS 84)
出力ラスタ　：　フォルダとファイル名を指定
OK
もう一度左上の右三角形の「ジオリファレンスの開始」をクリックするとGEOTIFFファイルが作成される。

計算して地図が歪んでいたら、その地点に赤い線が表示されます。（ポイントの追加した地点）
この線が長いと使い物になりません。（長ければ長いほど誤差が大きい）

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/cc94b3d5-60c2-8e69-7fdb-65b27a1eabad.png)

#####QGIS Desktopで作成したtifファイルを開く
ラスタレイヤの追加で、作成したtifファイルを開く。

マウスカーソルが手のひらになっているので、地物情報表示をクリックするとマウスカーソルが矢印になりました。

これで、マウスカーソルを緯度線と経度線の交点に持っていくと、ウインドウ下側の座標がだいたい良い値になっています。
ちょっとはずれている気がする。

確認
上記の地図で、四国のシェープファイルを作成し地図に描画しましたが、ポイントを決める時に適当にやったのが悪いのか、少しずれています。
上記の地図自体が範囲が広いせいか、resolution =の精度が悪いのか、ポリゴンが荒いです。

![image.png](https://qiita-image-store.s3.amazonaws.com/0/238800/b1d1e840-9f00-e308-267e-4f0a1f15e902.png)