1up noob: 2021

2021/09/14

【Python】matplotlibでDEMデータを表示する

国土地理院の基盤地図情報（数値標高モデル）などのDEMデータをその他のデータと組み合わせてmatplotlibでプロットしたい時に使えそうなプログラムが出来たので公開します。
例で使用したDEMデータは、国土地理院からダウンロードしたデータを予め結合したGeoTIFFを使用しています。

また、GeoTIFFの読み込みには、GDALを使用しています。

DEMデータの結合方法やGDALの導入については、色々なサイトで公開されているので、検索してください。本記事では、割愛します。

matplotlibでの標高塗り分けが本記事の中心部分となります。matplotlibのterrainでの塗り分けで満足出来ない場合にご活用いただければ幸いです。

スクリプト例

	import requests
	import numpy as np
	from osgeo import gdal
	from scipy.interpolate import griddata
	import matplotlib.pyplot as plt
	from matplotlib.colors import Normalize
	from matplotlib.colors import ListedColormap

	# 標高モデル
	ds = gdal.Open('./merge.tif', gdal.GA_ReadOnly)
	elevation = ds.GetRasterBand(1).ReadAsArray()

	# 切り出したい領域を設定
	new_x_coord = np.linspace(139.5, 140.5, 910)
	new_y_coord = np.linspace(35, 35.5, 550)
	xx, yy = np.meshgrid(new_x_coord, new_y_coord)
	nrows, ncols = elevation.shape
	x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()
	x1 = x0 + dx * ncols
	y1 = y0 + dy * nrows
	ele_xx, ele_yy = np.meshgrid(np.linspace(x0, x1, ncols), np.linspace(y1, y0, nrows))
	ele_zz = griddata(points=np.stack([ele_xx.reshape(-1,), ele_yy.reshape(-1,)], 1), values=np.flipud(elevation).reshape(-1,), xi=(xx, yy), method='nearest')

	# カラーマップを取得して、設定
	pg = requests.get('http://soliton.vm.bytemark.co.uk/pub/cpt-city/wkp/schwarzwald/wiki-schwarzwald-cont.pg')
	clist = [list(map(int, txt.split())) for txt in pg.text.splitlines()][::-1]
	newcmp = ListedColormap(np.array(clist)[:, 1:]/255)

	fig = plt.figure(figsize=(4, 3), dpi=300)
	plt.rcParams['font.size'] = 5
	ax1 = plt.subplot(111)
	# 標高モデル
	ax1.pcolormesh(xx, yy, np.where(ele_zz > -9999, ele_zz, np.nan), norm=Normalize(vmin=0, vmax=400), cmap=newcmp, shading='auto')
	ax1.contour(xx, yy, ele_zz, levels=[0, 100, 200, 300, 400], colors=['black'], linewidths=0.1)

	# 目盛設定
	ax1.set_xticks(np.linspace(139.5, 140.5, 11))
	xlabels = [str(round(label, 3)) + '°E' for label in np.linspace(139.5, 140.5, 11).tolist()]
	ax1.set_xticklabels(xlabels)
	ax1.set_xlim(139.5, 140.5)
	ax1.set_yticks(np.linspace(35, 35.5, 6))
	ylabels = [str(round(label, 3)) + '°N' for label in np.linspace(35, 35.5, 6).tolist()]
	ax1.set_yticklabels(ylabels)
	ax1.set_ylim(35, 35.5)
	plt.show()

view raw matplotlib_dem.py hosted with ❤ by GitHub

出力例

matplotlibのカラーマップに好みのものがなかったので、ArcGISっぽい?カラーマップを外部サイトから読み込んで使用しています。

カラーマップのファイルは、ローカルから直接読み込んで使用した方が良いと思いますので、上記のコードを利用される際はファイルをダウンロードして、コードを書き換えてください。

	import requests
	import numpy as np
	from matplotlib.colors import ListedColormap
	pg = requests.get('http://soliton.vm.bytemark.co.uk/pub/cpt-city/wkp/schwarzwald/wiki-schwarzwald-cont.pg')
	clist = [list(map(int, txt.split())) for txt in pg.text.splitlines()][::-1]
	newcmp = ListedColormap(np.array(clist)[:, 1:]/255)
	newcmp

view raw create_schwarzwald_cmap.py hosted with ❤ by GitHub

関連記事

2021/08/06

【Python】複数画像をPDFに変換する方法の比較

以前も、複数画像をPDFに変換する方法に関して、img2pdfの紹介記事を書いたのですが、img2pdfやそれ以外の手法での実行時間やファイルサイズなどを比較してみました。

本記事では、img2pdfを用いてpdfにする方法と、reportlabとPillowを組み合わせる方法、Pillowのみを使う方法の三種類を比較した結果を紹介します。

実行時間は、実行環境などで変わるため、あくまで目安です。計測に関しても、Jupyter Lab上で、実行した結果であり、複数回測って平均をとるというような厳密なことはしておりません。

実行環境は、以下の通りです。

MacBook Pro 13-inch (2017)
------
Python 3.8.11
img2pdf 0.4.1
reportlab 3.3.0
Pillow 8.1.0

下記のスクリプトは、いずれも書籍をスキャンしたpngファイル（連番ファイル名）を、pdfに変換することを想定しています。

Method1 : img2pdfを使う

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw convert_pdf_test1.ipynb hosted with ❤ by GitHub

Method2 : reportlabとPillowを組み合わせる

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw convert_pdf_test2.ipynb hosted with ❤ by GitHub

Method3 : Pillowだけを使う

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw convert_pdf_test3.ipynb hosted with ❤ by GitHub

10個の異なるpng画像ファイル群での実行時間と出力ファイルサイズの比較は、以下の通りです。

img2pdfを用いた方法は、実行時間（変換にかかる時間）が特に短いです。出力されたpdfファイルサイズは、pngファイルの合計サイズよりとほぼ同じサイズになっており、pngファイルを複雑に処理していないことが伺われます。 img2pdfに関して、残念なのはアルファチャンネルのある透過pngに対応していない点です。事前に処理しておいてね、という仕様となっています。

reportlabとPillowを組み合わせる方法ですが、実行時間が特に長く、画像のみのpdfを作るのであれば、reportlabを組み合わせる必要はないかなと感じました。

Pillowだけを用いる方法は、入力ファイルによって結果がばらついていますが、基本的にはファイルサイズが圧縮されています。リサンプリングしているため実行時間が伸びますが、この方法では、qualityの値を調整することでファイルサイズをさらに圧縮することも可能なので、ファイルサイズを調整したい時には良い方法だと思います。
ただし、img2pdfのように、見開きの設定など複雑なことはできないようなので、割り切って使うか他のツールと組み合わせる必要があります。

追記:2022/12/09
PyMuPDFを使って、変換する方法が分かったので追加しました。実行速度などは測っていません。

Method4 : PyMuPDFを使う

	import os
	import fitz

	file_dir = '(画像ファイルが入っているフォルダのパス)'
	file_ext = '.png'
	img_files = sorted([os.path.join(file_dir, _) for _ in os.listdir(file_dir) if _.endswith(file_ext)])

	# PDFの新規作成
	doc = fitz.open()
	# PDFにPNG画像を書き込む
	for file in img_files:
	# 画像を読み込む
	with fitz.open(file) as img:
	rect = img[0].rect
	pdfbytes = img.convert_to_pdf()
	# 新しいページを追加
	page = doc.new_page(width=rect.width, height=rect.height)
	page.show_pdf_page(rect, fitz.open('pdf', pdfbytes), 0)

	# PDFを保存
	doc.save('output4.pdf')

view raw PyMuPDF_png_pdf.py hosted with ❤ by GitHub

参考：

Method5 : borbを使う

borbを使って変換することも出来ましたが、v2.1.7では画像をPDFにする際に自動でOCR処理をするようで変換処理が遅いです。OCR処理をスキップ指定できる方法が分かれば更新します。基本的には余白が自動で入る設計で、このような変換は想定していないと思われます。

	import os
	import borb
	from PIL import Image

	file_dir = '(画像ファイルが入っているフォルダのパス)'
	file_ext = '.png'
	img_files = sorted([os.path.join(file_dir, _) for _ in os.listdir(file_dir) if _.endswith(file_ext)])

	doc = borb.pdf.Document()
	for file in img_files:
	img = Image.open(file)
	page = borb.pdf.Page(width=img.width, height=img.height)
	doc.add_page(page)
	page_layout = borb.pdf.SingleColumnLayoutWithOverflow(page, horizontal_margin=0, vertical_margin=0)
	page_layout.add(borb.pdf.Image(img, width=img.width, height=img.height, margin_bottom=0, margin_left=0, margin_right=0, margin_top=0, border_width=0))

	with open('output5.pdf', 'wb') as pdf_file_handle:
	borb.pdf.PDF.dumps(pdf_file_handle, doc)

view raw borb_png_pdf.py hosted with ❤ by GitHub

続き

2021/05/10

【Python】foliumで地質図を表示する方法

先日「【Python】foliumで活断層を表示する方法」という記事を書いたばかりですが、今回もfoliumネタです。

foliumは、LeafletというJavaScriptライブラリをPythonから使える様にしたものですので、地図タイルとして様々な地図が利用できます。有名どころでは、国土地理院の地理院地図の地理院タイルなどがあります。そのほかにも、foliumにビルドインされている地図もあります（参考）。

本記事では、日本シームレス地質図を地図タイルとして使う際のサンプルコードを掲載します。

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw folium-seamless-geological-map.ipynb hosted with ❤ by GitHub

出力参考：

地理院地図を使う際のサンプルコードも併せて掲載します。

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw folium-GSI-tiles.ipynb hosted with ❤ by GitHub

出力参考：

地理院タイルでは、ズームレベルの制限などもありますが、英語版や土地利用図、火山基本図、衛星画像など様々なタイルが公開されています。

利用にあたっては、承認などが必要な場合もある様ですのでご注意下さい。

参考：「国土地理院の地図の利用手続」

【Python】foliumで活断層を表示する方法

昔書いた記事ですが、「Google Earth上に活断層を表示する方法」にそこそこアクセスがあるので、同様に、産総研が公開している活断層データベースのkmzファイルをfoliumを使って表示するプログラムを書いたので公開します。

kmzを開いて、kmlの中身をPolylineとして表示するプログラムとなっています。

kmzファイルの仕様を確認して作成したわけではないので、kmzファイルによってはそのまま使用できないと思いますが、参考になれば幸いです。

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw unzip_kmz.ipynb hosted with ❤ by GitHub

出力参考：

登録: 投稿 (Atom)

2021/09/14

【Python】matplotlibでDEMデータを表示する

スクリプト例

出力例

2021/08/06

【Python】複数画像をPDFに変換する方法の比較

Method1 : img2pdfを使う

Method2 : reportlabとPillowを組み合わせる

Method3 : Pillowだけを使う

Method4 : PyMuPDFを使う

Method5 : borbを使う

2021/05/10

【Python】foliumで地質図を表示する方法

2021/05/08

【Python】foliumで活断層を表示する方法