ラベル GIS の投稿を表示しています。 すべての投稿を表示
ラベル GIS の投稿を表示しています。 すべての投稿を表示

2022年6月14日火曜日

Pythonで海面水温衛星データを可視化する

G-Portalという地球観測衛星データ提供システムがあり、登録するとさまざまな衛星データを取得できる。今回は気候変動観測衛星「しきさい」（GCOM-C）の海面水温プロダクト(SST) を取得して可視化してみる。


環境

WSL2(Ubuntu20.04)とJupyter。

$ lsb_release -dr
Description:    Ubuntu 20.04.4 LTS
Release:        20.04
$ python3 -V
Python 3.8.10
$ jupyter --version
jupyter core     : 4.7.1
jupyter-notebook : 6.2.0
qtconsole        : 5.0.2
ipython          : 7.20.0
ipykernel        : 5.4.3
jupyter client   : 6.1.12
jupyter lab      : 3.0.9
nbconvert        : 6.0.7
ipywidgets       : 7.6.3
nbformat         : 5.1.2
traitlets        : 5.0.5


G-Portalの登録

衛星データをダウンロードするために必要なユーザ登録をしておく。ユーザ登録画面から、アカウント名やメールアドレスなどを登録する。


衛星データのダウンロード

G-Portalトップページのメニューからログインにアクセスして、登録したユーザーアカウントとパスワードでログインする。［衛星からの検索］をクリックして、ダウンロードするデータを選択する画面に遷移する。


1.絞り込み

ダウンロードする衛星データの種類を指定する。

［1.絞り込み］タブの［衛星、センサから選ぶ］を選択する。そして、［GCOM-C/SGLI］ ＞ ［LEVEL2］ ＞ ［海洋圏］ ＞ ［L2-海面水温］のチェックボックスをチェックする。


2.期間指定

取得するデータを検索する期間を指定する。

［2.期間指定］タブで期間を指定する。


3.範囲指定

取得するデータの範囲を指定する。

［3.範囲指定］タブを開いて、右側に表示されている地図で指定するか、緯度・経度を数値で指定する。


検索するデータの指定ができたら［検索する］ボタンを押す。対象データの一覧が表示されるので、使用するデータをダウンロードする。ダウンロードしたファイル名はGC1SG1_L2SG_SSTDQ_3001.h5とする。

ダウンロードしたファイルはHDF形式で、可視化する際にはGeoTIFFが必要になる。データ一覧の［加工］ボタンを押すと、GeoTIFFへの加工要求ができる。加工には少し時間がかかるので、今回は後述するツールでGeoTIFFに変換する。


衛星データをGeoTIFFに変換する

ダウンロードしたHDF形式の衛星データをGeoTIFFに変換する。

まずはG-Portalから変換ツールをダウンロードする。トップページから［ツール・ドキュメント］画面へ遷移し、［GCOM-C］を選択。SGLI 地図投影・GeoTIFF出力ツール(Ver.1.1)をダウンロードする。ここではLinux版を使う。

ファイル名をsgli_geo_map_linuxとしてダウンロードし、先にダウンロードした衛星データを、WSL2上の同じディレクトリに配置する。

変換する前に、ダウンロードした衛星データ内のデータセット名を確認しておく。

! pip3 install h5py

import h5py

def print_dataset_name(name, obj):
    if isinstance(obj, h5py.Dataset):
        print(name)

with h5py.File('./GC1SG1_L2SG_SSTDQ_3001.h5', mode='r') as f:
    f.visititems(print_dataset_name)


上記コードを実行すると、以下のようにHDFファイル内のデータセット一覧が表示される。

Geometry_data/Latitude
Geometry_data/Longitude
Geometry_data/Obs_time
Geometry_data/Sensor_zenith
Geometry_data/Solar_zenith
Image_data/Cloud_probability
Image_data/Line_tai93
Image_data/QA_flag
Image_data/SST

上記データセットのうち、今回使用する海面水温プロダクト(SST)はImage_data/SST。変換時にはこのデータセットを指定する。

以下コマンドでHDF形式からGeoTIFFに変換する。はじめにダウンロードしたツールに実行権限を付与しておく。変換されたGeoTIFFのファイル名はGC1SG1_L2SG_SSTDQ_3001_SST.tif となる。

$ chmod u+x sgli_geo_map_linux
$ ./sgli_geo_map_linux GC1SG1_L2SG_SSTDQ_3001.h5 -d Image_data/SST


SSTデータの値を摂氏（℃）に変換

SSTの値は一般的に使われる温度の単位（摂氏）ではないので、摂氏に変換する必要がある。しきさい（GCOM-C）のページにSSTデータの仕様とその変換方法が記載されている。トップページから［データを使う］＞［標準プロダクト&アルゴリズム］を開き、Ver.3の海面水温を選択する。開いた画面のAttribute情報にあるSSTを開く。

この仕様によると、有効な値は0~65531で、65532～65535は陸地や雲、その他のエラーなどで水温が取得できなかった場合。SSTの値を摂氏（℃）に変換する式は以下の通り。DNがSSTの値で、Slopeは0.0012、Offsetは-10となっている。

SST[degree]=DN*Slope+Offset


海面水温データを可視化する

ダウンロードした海面水温衛星データを、Jupyterで可視化してみる。地理空間データ用のPythonライブラリであるrasterioとEarthPyを使う。

! pip3 install rasterio # ラスター形式の地理空間データを扱うためのPythonライブラリ
! pip3 install earthpy # ラスターおよびベクター形式空間データを扱うためのPythonライブラリ

import rasterio
from earthpy.mask import mask_pixels
from earthpy.plot  import plot_bands

file_path = './GC1SG1_202206100147B05610_L2SG_SSTDQ_3001_SST.tif'
sst = rasterio.open(file_path)

# バンド番号（1からはじまる）を指定してデータ読み込み
# https://rasterio.readthedocs.io/en/latest/topics/reading.html
sst_array = sst.read(1)
print(sst_array.shape)

# SSTの非有効値をマスクする
# https://earthpy.readthedocs.io/en/latest/api/earthpy.mask.html
mask_values = [65532, 65533, 65534, 65535] # マスクする値
sst_masked = mask_pixels(sst_array, sst_array, vals=mask_values)

# SSTの値を摂氏（℃）に変換
slope = 0.0012
offset = -10.0
sst_deg = sst_masked * slope + offset
sst_deg_min = sst_deg.min() # 最低水温
sst_deg_max = sst_deg.max() # 最高水温

# 海面水温を可視化
plot_bands(sst_deg, vmin=sst_deg_min, vmax=sst_deg_max, cmap='rainbow')

以下のような画像が表示される。陸地や雲などで有効値がない箇所は白色になっている。


投稿者 時刻: 0 件のコメント:
ラベル: , , ,

2022年4月22日金曜日

PythonとCartopyで地図データを可視化する

PythonでGIS用シェープファイルを可視化するでは、日本地図とフェリー航路のシェープファイルをPythonライブラリのGeoPandasで可視化した。日本地図のシェープファイルは国土地理院のホームページからダウンロードしたが、地理空間データ生成用のPythonパッケージであるCartopyを使うと、自分で地図のシェープファイルを用意する必要がない。CartopyにはNatural Earthのシェープファイルをロードするメソッドが用意されている。今回はCartopyを使って地図データを可視化してみる。


環境

WSL2(Ubuntu20.04)とJupyter。

$ lsb_release -dr
Description:    Ubuntu 20.04.4 LTS
Release:        20.04
$ python3 -V
Python 3.8.10


Pythonライブラリのインストール

必要なPythonライブラリをインストールする。Cartopyのインストールは少し手間がかかるので後述。

$ pip3 install matplotlib geopandas

インストールされたバージョン。

$ pip show matplotlib
Name: matplotlib
Version: 3.5.1
...
$ pip show geopandas
Name: geopandas
Version: 0.10.2
...


Cartopyのインストール

Cartopyはpipでインストールというわけにはいかず、少し手間がかかる。Cartopyをインストールする前に、proj、geos、pyshp、shapelyをインストールしておく必要がある。

まずはprojのインストール。projはaptコマンドでもインストールできるが、バージョンが古いのでソースからインストールする。ビルドの前に必要なパッケージをインストールする。

$ sudo apt install cmake sqlite3 libsqlite3-dev 

続いてソースファイルをダウンロードしてビルドする。今回はバージョン8.2.1を使う。

$ curl -O https://download.osgeo.org/proj/proj-8.2.1.tar.gz
$ tar xzvf proj-8.2.1.tar.gz
$ cd proj-8.2.1
$ mkdir build
$ cd build
$ cmake ..
$ cmake --build .
$ sudo cmake --build . --target install

次にgoesのインストール。

$ sudo apt install libgeos-dev
$ pip3 install geos

pyshpとshapelyのインストール。

$ pip3 install pyshp shapely

インストールした各ライブラリのバージョン。

$ pip3 show geos
Name: geos
Version: 0.2.3
...
$ pip3 show pyshp
Name: pyshp
Version: 2.2.0
...
$ pip3 show shapely
Name: Shapely
Version: 1.8.1.post1
...

最後にCartopyのインストール。

$ pip3 install cartopy

インストールされたバージョン

$ pip3 show cartopy
Name: Cartopy
Version: 0.20.2
...


CartopyとGeoPandasで地図とフェリー航路を可視化する

地図の他に、PythonでGIS用シェープファイルを可視化すると同様にフェリー航路も可視化する。フェリー航路のシェープファイルは、国土地理院のホームページにある地球地図日本の第2.2版ベクタの全レイヤを使う。ファイルをダウンロードして解凍しておく。フェリー航路のシェープファイルはferryl_jpn.shp。

地図とフェリー航路を可視化するコードは以下の通り。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import geopandas as gpd

# 解凍したフェリー航路のシェープファイルのパス
path_ferry =  "gm-jpn-all_u_2_2/ferryl_jpn.shp"
df_ferry = gpd.read_file(path_ferry)

plt.figure(figsize=(20, 16))

# 地図投影法を指定
ax = plt.axes(projection=ccrs.PlateCarree())

# Natural Earthの地図シェープファイル（解像度10m）を読み込む
ax.coastlines(resolution='10m')

# フェリー航路をプロット
df_ferry.plot(ax=ax, color='blue', linewidth=1.0, label='Ferry routes')

以下のように地図とフェリー航路が表示される。クリッピングされる地図の位置は、読みこんだシェープファイル（今回はフェリー航路）の位置により自動調整される。


北海道が途中で切れているので北海道全体がクリッピング範囲に入るようにしてみる。Cartopyにはset_extentというメソッドが用意されていて、任意の場所でクリッピングできるらしいのだが、実行しても処理が途中で落ちてしまう。ここでは回避策として、透明なポリゴンの境界を配置して、クリッピング位置が自動調整されたときに北海道が切り取られないようにする。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import geopandas as gpd
from shapely.geometry import Polygon

# 解凍したフェリー航路のシェープファイルのパス
path_ferry =  "gm-jpn-all_u_2_2/ferryl_jpn.shp"
df_ferry = gpd.read_file(path_ferry)

plt.figure(figsize=(20, 16))

# 地図投影法を指定
ax = plt.axes(projection=ccrs.PlateCarree())

# Natural Earthの地図シェープファイル（解像度10m）を読み込む
ax.coastlines(resolution='10m')

# フェリー航路をプロット
df_ferry.plot(ax=ax, color='blue', linewidth=1.0, label='Ferry routes')

# 指定範囲で地図を切り取る
# set_extentがエラーになるので透明なポリゴンの境界で対処
# ax.set_extent([130, 150, 30, 40], crs=ccrs.PlateCarree())
polygon = Polygon([(130, 40), (150, 40), (150, 30), (130, 30), (130, 40)])
poly_gdf = gpd.GeoDataFrame([1], geometry=[polygon], crs=df_ferry.crs)
poly_gdf.boundary.plot(ax=ax, alpha=0)

以下のように北海道が途中で切れることなくプロットされた。


投稿者 時刻: 0 件のコメント:
ラベル: , , ,

2022年1月24日月曜日

Raspberry PiにGeoPandasをインストールする

PythonでGISデータを扱うときにはGeoPandasが使える。例えばPythonでGIS用シェープファイルを可視化するやPythonでシェープファイルから都道府県境界図を作成するなど。名前の通りGISデータに対応したPandasのようなもので、GISデータをPandasのデータフレームと同じように操作できる。

しかしRasoberry Piでpipを使ってインストールしようとしても失敗する。どうやら依存関係があるpyprojのインストールでエラーになっている。そこでRaspberry PiにGeoPandasをインストールする手順をまとめておく。


環境

Raspberry Pi OS(bullseye)。

$ lsb_release -dr
Description:    Raspbian GNU/Linux 11 (bullseye)
Release:        11
$ python3 -V
Python 3.9.2


GeoPandasをインストールする方法

pipによるインストールで失敗するのでaptコマンドを使う。aptだとすんなりインストールできる。

$ sudo apt install python3-geopandas
$ pip3 show geopandas
Name: geopandas
Version: 0.8.2

ただaptコマンドでインストールできるバージョンは0.8.2で最新の0.10.2と比べるとけっこう古い。そこでpipコマンドでアップグレードする。これで無事GeoPandasの最新バージョンがインストールできた。

$ pip3 install geopandas --upgrade
$ pip3 show geopandas
Name: geopandas
Version: 0.10.2


投稿者 時刻: 0 件のコメント:
ラベル: , , ,

2020年1月28日火曜日

Pythonでシェープファイルから都道府県境界図を作成する

PythonでGIS用シェープファイルを可視化するでは、PythonライブラリのGeoPandasとMatplotlibでシェープファイルを読みこんで日本地図を表示した。今回は、シェープファイルをダウンロードして、何かと使うことのありそうな都道府県境界図をPythonで作成してみる。


環境


Windows10(1903)のWSL(Ubuntu 18.04)とJupyter Notebookを使用。



境界図のシェープファイル


都道府県境界図を作成するためにシェープファイルをダウンロードする。国土数値情報　ダウンロードサービスや地図で見る統計(統計GIS)からダウンロードできる。国土数値情報からダウンロードできるシェープファイルだと全国のデータが1ファイルにまとまっているので作成が簡単だと思ったのだが、都道府県境図にしようとしても市町村区境の図になってしまった。そこで、今回は地図で見る統計(統計GIS)から都道府県ごとのシェープファイルをダウンロードして、それを連結して都道府県境図を作成する。


シェープファイルのダウンロード


ダウンロードするシェープファイルは地図で見る統計(統計GIS)の境界データ。今回は以下の場所にあるファイルを使う。
境界データダウンロード＞小地域＞国勢調査＞2015年＞小地域
＞世界測地系緯度経度・Shapefile
都道府県ごとに全域のファイルがあるのでそれをダウンロードする。全47ファイルをダウンロードしたら、zip圧縮されているのですべて解凍して同じディレクトリに保存しておく。

また、以下の場所から定義書がダウンロードできるので、これをみるとどんな項目があるか確認できる。
境界データダウンロード＞小地域＞国勢調査＞2015年＞小地域
今回使うのは次の3項目。
  • PREF_NAME 都道府県名
  • AREA 面積(㎡)
  • JINKO 人口

ちなみに、北海道全域のファイルをそのまま表示すると以下のような国勢調査の小地域の境界図となる。


GeoPandasで都道府県境図を作成


まずはGeoPandasとそのほか必要なパッケージなどをインストールする。


インストールされたバージョン。


準備ができたので、GeoPnadasで都道府県境界図を作成する。手順としては、ダウンロードした都道府県ごとのシェープファイルは小地域境のデータなので、これを都道府県境のデータにする。GISではこういう処理をディゾルブと言うらしく、GeoPandasではそのためのメソッドが用意されている（Aggregation with dissolve）。ここでは都道府県境界図を作成するのでPREF_NAME（都道府県名）でディゾルブする。また、人口などの数値はaggfuncで指定した集計がされる。市区町村境界図を作成する場合はCITY_NAME（区町村名）でディゾルブすればよい。ディゾルブしたデータを連結していき、最終的に日本地図とする。これをシェープファイルとして保存。以下はそのPythonコード。

%matplotlib inline
import matplotlib.pyplot as plt

import os
from pathlib import Path
import glob

import pandas as pd
import geopandas as gpd

# シェープファイル保存ディレクトリ
SHAPE_DIR = os.path.join(str(Path.home()), 'estat/sgis_shape')

def create_outline_map():
    print('*** Creating outline map ***')
    print('Loading from {}'.format(SHAPE_DIR))

    jpn = None

    shpfiles = glob.glob(os.path.join(SHAPE_DIR, '*.shp'))
    for shp_path in shpfiles:
        gdf = gpd.read_file(shp_path)

        # byでディゾルブする単位を指定する
        # byで指定したカラムがindexになる
        # AREAとJINKOは後でコロプレス図で使用
        pref = gdf[['PREF_NAME', 'AREA', 'JINKO', 'geometry']].dissolve(by='PREF_NAME', aggfunc='sum')

        if jpn is None:
            jpn = pref
        else:
            # GeoDataFrameの結合
            jpn = gpd.GeoDataFrame(pd.concat([jpn, pref]), crs=jpn.crs)

    # シェープファイルにはindex情報は保存されないので都道府県名をカラムにしておく
    jpn['PREF_NAME'] = jpn.index

    print()
    print(jpn.info())
    print(jpn.head())
    
    # shapeファイル保存
    jpn.to_file('jpnmap.shp', encoding = 'shift-jis')

    # 境界図のプロット
    fig, ax = plt.subplots(figsize = (20,16)) 
    jpn.plot(ax=ax, facecolor='white', edgecolor='black', linewidth=0.5)
    ax.set_axis_off()
    plt.axis('equal')
    plt.show()

if __name__ == '__main__':
    create_outline_map()


以下のような都道府県境界図が作成された。


コロプレス図を作成する


ダウンロードしたファイルには面積と人口のデータも含まれている。せっかく境界図のシェープファイルを作成したので、人口密度を計算してコロプレス図を作成する。

コードは以下の通り。
%matplotlib inline
import matplotlib.pyplot as plt

import geopandas as gpd

def plot_choropleth():
    # 作成した都道府県境界図のシェープファイルの読み込み
    jpn = gpd.read_file('jpnmap.shp', encoding = 'shift-jis')

    # AREA（面積）の単位は平方メートルなので平方キロメートルに直して人口密度を計算
    jpn['DENSITY'] = jpn['JINKO'] / (jpn['AREA'] / (10**6))

    print()
    print(jpn.info())
    print(jpn.head())

    fig, ax = plt.subplots(figsize = (12,8))

    # 人口密度のコロプレス図
    jpn.plot(column = 'DENSITY', edgecolor = "black",
                    scheme='quantiles', cmap='YlOrRd', ax=ax, 
                    legend = True)

    ax.set_axis_off()
    plt.axis('equal');
    plt.show()

if __name__ == '__main__':
    plot_choropleth()

次のようなコロプレス図が作成された。

投稿者 時刻: 0 件のコメント:
ラベル: , , ,

2019年5月22日水曜日

PythonでGIS用シェープファイルを可視化する

GISで使われる地図データなどをPythonで利用する機会があったので、調べたことをまとめておく。最終的には、シェープファイル形式のデータをPythonで読み込んで、matplotlibでプロットしてみる。


環境


Windows10のWSL(Ubuntu 18.04)。



GISで使われるデータ


GISで使われるデータフォーマットにはいくつか種類があって、まずはラスターデータとベクターデータに分けられる。ラスターデータはざっくり言うとjpgやpngなどの画像で、ベクターデータは座標や線などの情報からなるデータ。PDFで拡大しても画像のように線などが荒くならないのは、PDFではベクターデータを扱えるから。GISで使われるベクターデータにもいくつかフォーマットがあって、そのひとつがシェープファイル。シェープファイルは本体の拡張子shpのファイルとその他の複数のファイルから構成される。


可視化するデータ


国土地理院のホームページにある地球地図日本の第2.2版ベクタの全レイヤをダウンロード。ダウンロードしたファイルを解凍すると、いくつかのシェープファイルがある。今回は日本の海岸線データ（coastl_jpn.shp）とフェリー航路データ（ferryl_jpn.shp）を可視化してみる。




GIS用データを扱えるPythonライブラリ


調べた限りGIS用データを扱えるPythonライブラリはいくつかあるが、GeoPandasを使うことにした。Pandasという名前がついているように、GeoPandasではPandasと同じデータ形式が使える。Pandasに慣れていれば他のライブラリよりも扱いやすいかもしれない。

GeoPandasはpipでインストールできる。



シェープファイルの可視化


ダウンロードしたシェープファイルをmatplotlibで可視化する。

import geopandas as gpd

# matplotlibのターミナル対応
import matplotlib
matplotlib.use('Agg')

import matplotlib.pyplot as plt

def main():
    # シェープファイルのあるディレクトリ
    shape_dir = 'gm-jpn-all_u_2_2/'
    # シェープファイル名
    shape_fname_l = ['coastl_jpn.shp', 'ferryl_jpn.shp']
    # プロットする色
    color_l = ['black', 'blue']

    fig, ax = plt.subplots(figsize = (20,16)) 
    for shape_fname, color in zip(shape_fname_l, color_l):
        # シェープファイルの読み込み
        data = gpd.read_file(shape_dir + shape_fname)
    
        data.plot(ax=ax, color=color)

    plt.savefig('gis_coast_ferry.pdf')
    plt.clf()

if __name__ == '__main__':
    main()

結果は以下の通りで、日本の海岸線は黒、フェリー航路は青でプロットした。


投稿者 時刻: 0 件のコメント:
ラベル: , , ,
登録: 投稿 (Atom)