仮想日本を作ってみた

気候区分の説明をするのに使われる仮想大陸っていうものがあるのをご存知でしょうか。

あの、魚みたいな形をしたやつです。

仮想大陸というのは、緯度帯ごとの面積占有率が実際の大陸分布と等しいんですね。

どういうことかというと、すべての大陸・島を、緯度を変えずにぎゅっと集めたものなんです。

記事を書こうと思って仮想大陸のことを考えていた時、ふとこれを日本列島で作ってみたらどうなるんだろう、と思ったんです。

1 構想
2 陸の幅を取得
3 エクセルで処理
4 シェープファイル作成
5 完成
6 あとがき
- 6.1 共有:

構想

まずはどうやって作るかを考えなければなりません。

まず、仮想大陸とは緯度帯ごとの面積占有率が変わらないように等積変形し一カ所に集めたものと定義します。

ただ日本列島の形を変形しただけではつまらないので、各都道府県も変形したいですね。

QGISのなかに、「線長の合計」というベクタツールがあるのでこれを使うことにしました。

陸の幅を取得

緯度帯ごとにどれだけの長さ陸があるかを調べればよいわけです。

ここで、どれぐらいの解像度で求めるか決定しなくてはなりませんが、今回は0.01°刻みにしました。

日本列島の南北方向の広がりがおよそ25°程度なのでラスタデータを想定すると縦が2500pxあることになり、十分だと考えました。

「線長の合計」はポリゴンレイヤとラインレイヤの共通部分の長さを属性として持つ、ポリゴンレイヤの複製を返すツールなのでまず都道府県ごとのポリゴンレイヤを用意します。

国土数値情報を使って各都道府県の行政区分のshpファイルをダウンロードします。

今回は行政界を使わないので融合します。47都道府県もやるのは大変なのでpythonスクリプトにやらせましょう。

import processing

for n in range(47):
    prefn = str(n+1).zfill(2)
    #nは0始まりなのに対し都道府県コードは01始まりなので、nに1を加え2桁で0埋めする
    
    input = '~/N03-20_'+prefn+'_200101.shp'
    output = '~/dissolved_'+prefn+'.shp'
    #最新の行政区域のシェープファイル名は「N03-20_(都道府県コード)_200101.shp」だった

    processing.run("native:dissolve",{'INPUT': input, 'OUTPUT': output})

各都道府県の北端と南端の座標をここから取得し、0.01°単位で北端を切り上げ、南端を切り下げて各都道府県の処理開始の緯度を記したリストと処理停止の緯度を記したリストを作ります。

開始の緯度から順に緯線の一時ファイルを出力していき、読み込んだポリゴンファイルと合わせて線長の合計を求め、新しいレイヤの属性である線長をcsvファイルに吐き出させました。

import csv
from qgis.core import *
import processing

latitudemin=[41.35,40.21,38.73,37.76,38.86,37.73,36.78,35.73,36.18,35.98,35.75,34.88,20.41,35.11,36.73,36.26,36.06,35.33,35.16,35.18,35.13,34.56,34.56,33.71,34.78,34.7,34.26,34.15,33.85,33.41,35.05,34.3,34.28,34.03,33.7,33.53,34,32.88,32.7,33,32.95,31.98,32.08,32.7,31.35,27.01,24.03]
latitudemax=[45.57,41.57,40.47,39.02,40.52,39.22,37.99,36.95,37.17,37.07,36.3,36.12,35.9,35.69,38.57,36.99,37.87,36.3,35.99,37.04,36.47,35.65,35.44,35.27,35.72,35.79,35.07,35.69,34.77,34.4,35.62,37.25,35.37,35.12,34.8,34.27,34.57,34.32,33.9,34.27,33.64,34.74,33.2,33.75,32.85,32.32,27.9]
#各都道府県の緯度の下限・上限を格納したリスト

for n in range(47):
    prefn = str(n+1).zfill(2)
    #融合の時と同じ
    input = '~/dissolved_'+prefn+'.shp'
    dlayer = QgsVectorLayer(input, 'polygonlayer', 'ogr')
    #2つめの引数はレイヤ名（なんでもいい）、3つめはプロバイダの指定（固定）
    long = round((latitudemax[n]-latitudemin[n])*100+1)
    #計算を繰り返す数、つまり下限から上限まで0.01刻みで何ステップあるかを計算している
    #何故か小数になったので整数に丸めた
    list = []
    for a in range(long):
        lat = a*0.01 + latitudemin[n]
        #aの初期値は0、繰り返すごとに0.01ずつ増える
        line_start = QgsPoint(100,lat)
        line_end = QgsPoint(150,lat)
        #線分の2端を設定、日本が含まれればいいので経度は適当に100°Eと150°Eにした
        line = QgsGeometry.fromPolyline([line_start,line_end])
        #2端を結ぶ線分（ジオメトリ）を作成
        linelayer = QgsVectorLayer('LineString', 'line', 'memory')
        #新規一時レイヤ（3つめのoutputの引数を'memory'にする）を作成
        dp= linelayer.dataProvider()
        #dataproviderというものを介してレイヤはジオメトリなどと接続する
        feat = QgsFeature()
        #新規フィーチャーの作成
        feat.setGeometry(line)
        #featフィーチャーにlineジオメトリを追加
        dp.addFeatures( [feat] )
        #dataproviderによりレイヤにフィーチャーを追加        
        output = '~/'+prefn+'.gpkg'
        #こちらも一時レイヤにするとなぜかエラーが出てしまったので適当にレイヤを出力させる
        #シェープファイルでもジオパッケージでもいい
        len = processing.run('qgis:sumlinelengths', { 'COUNT_FIELD' : 'COUNT', 'LEN_FIELD' : 'LENGTH', 'LINES' : linelayer, 'OUTPUT' : output, 'POLYGONS' : dlayer })
        #「線長の合計」ツールを使う。1つめ・2つめのレイヤは属性の名前を付ける。COUNTにはポリゴンと交差した線分の数、LENGTHには長さが格納される
        result_layer = QgsVectorLayer(len['OUTPUT'], "linelayer", "ogr")
        #できたレイヤを取得する
        for feature in result_layer.getFeatures():
            list.append(feature["LENGTH"])
            #feature名が不明なのでfor文を使った（ちゃんと調べればいい話）。listに追記する
    with open('C:/Users/shun6/Desktop/GIS/result.csv', 'a') as f:
            writer = csv.writer(f)
            writer.writerow(list)
            #csvファイルに追記

6時間くらいかかって全都道府県の処理が完了。

空間インデックスを作成すれば早く済むらしいですがコードを書くのが面倒なのでやめました。

エクセルで処理

日本全体での緯度の範囲とおなじようにエクセルの行に0.1°刻みで緯度を対応させ、各都道府県のデータを並べていきます。

データのない部分を0埋めしたものが以下の図です。

B列が北海道、C列が青森県…というようにA列の緯度に対応する線長が記されていますね。

得られた値が度数表示ではなく、距離表示だったのでそれぞれの緯度帯における1度あたりの平行圏弧長を求め、逆数を掛けて度数表示に改めます。

私が使ったのはWGS84楕円体でしたが、1度あたりの平行圏弧長が記された表などを見つけられなかったのでこれもPyQGISで調達しました。

やっていることは先ほどと同じで、経度差1°の幅を持つ長方形を作図し、線長を求めるだけです。

形が単純なだけあって瞬殺でした。

どのように都道府県を配置するかがここで決まるわけですが、面倒くさくなって都道府県コード順に東から並べることにしました。

兵庫→奈良とか明らかに逆行しているんですが面倒なので目をつむりましょう…

中心は標準時子午線の135°にしましょう。

ここで、エクセルの列を縦に見ると、48本の線が描かれて、一番東の間が北海道、次が青森県という風になるようにしました。

例えば緯度nの線上には北海道が0.5°、青森県が0.3°であったとすると、一番東の線が135.4°、次の線が134.9°、のこりの46本が134.6°というふうにするわけです。

それぞれの線ごとに、日本全体での緯度下限から緯度上限を0.01で区分けしたものに対応する経度を記したリストを作り、経度の二重リストとしてPythonに取り込みます。

シェープファイル作成

北海道の場合、1本目と2本目の間がポリゴンになるわけですが、1本目と2本目が重なっているところ、つまりポリゴンがない場所も併せてレイヤを作ってしまうとエラーが出ることが分かったので、1本目と2本目が重なる場所を区切りとして、ポリゴンの存在する場所をマルチパートとして読み込ませました。

import csv
from qgis.core import *
import processing
from qgis.PyQt.QtCore import QVariant
from qgis.utils import iface

with open('/~.csv', encoding='utf-8-sig') as f:
    reader = csv.reader(f)
    l = 
#csvファイルを読み込み、行ごとのリストをリストとします。
#1行目には緯度、2~49行目には対応する経度長が入っているcsvです
code = ["北海道","青森県","岩手県","宮城県","秋田県","山形県","福島県","茨城県","栃木県","群馬県","埼玉県","千葉県","東京都","神奈川県","新潟県","富山県","石川県","福井県","山梨県","長野県","岐阜県","静岡県","愛知県","三重県","滋賀県","京都府","大阪府","兵庫県","奈良県","和歌山県","鳥取県","島根県","岡山県","広島県","山口県","徳島県","香川県","愛媛県","高知県","福岡県","佐賀県","長崎県","熊本県","大分県","宮崎県","鹿児島県","沖縄県"]
for pref in range(47):
    pointlist=[]
    pointlist1=[]
    pointlist2=[]
    for d in range(2150):
        #2150は今回の下限～上限緯度のステップ数です
        lat = float(l[0][d])
        long = float(l[pref+1][d])
        pointlist1.append(QgsPointXY(long, lat))
        #それぞれの都道府県ポリゴンの東端の座標をpointlist1に格納します。
    for d in range(2150):
        lat = float(l[0][d])
        long = float(l[pref+2][d])
        pointlist2.append(QgsPointXY(long, lat))
        #それぞれの都道府県ポリゴンの西端の座標をpointlist2に格納します。
    p = 0
    #pはステップ数を表します。2149（0始まりなので）になれば処理を終了します。
    while p <= 2149:
        while p <= 2149:
            plist= []
            templist=[]
            ori= 0
            if pointlist1[p]==pointlist2[p]:
                p += 1
                #東端と西端の線が重なっていた場合の処理です。何もせず次の緯度へ進みます。
            elif p == 0:
                plist.append(pointlist1[p])
                plist.append(pointlist2[p])
                p += 1
                break
                #東端と西端の線が重ならず、かつp=0（つまり緯度下限の時）の処理です
                #今回は最南端の緯度帯に沖縄県が0.05°程度の経度長を持っていたためこの条件を付けました（日本の最南端は東京都だがうまく0.1°刻みの緯線に乗らなかったため計上されなかった）
                #緯度帯を増やして、最も南の都道府県でも0になるような緯度が最南端であれば必要ありません
            else:
                plist.append(pointlist1[p-1])
                ori = pointlist1[p]
                plist.append(pointlist2[p])
                p += 1
                break
                #2ステップ目以降に初めて重ならなくなった時の処理です。一致した最後の座標、西端の座標を順番にplistに格納し、東端の座標は一時的にoriに格納します
        while p <= 2149:
            if pointlist1[p]==pointlist2[p]:
                plist.append(pointlist1[p])
                templist.reverse()
                #西端と東端が一致したときは一致した座標をplistに格納します
                if ori==0:
                    p2list=plist+templist
                    pointlist.append(p2list)
                    p += 1
                    break
                    #oriが0、つまりポリゴン南端の緯度で座標が二つある場合(一つ上のループの"elif")はtemplistを逆転させてplistに付け加えます
                    #逆転させる理由は、ポリゴンを作成する手順で、リストの順に点を結んでいくためで、西端を南から北に付け加えたなら東端は北から南に付け加える必要があります
                else:
                    templist.append(ori)
                    p2list=plist+templist
                    pointlist.append(p2list)
                    p += 1
                    break
                    #oriが0以外、つまりポリゴンの南端で座標が1つしかない場合です。oriの座標を逆にしたtemplistの最後尾に付け加えplistに付加します
            elif p == 2149:
                plist.append(pointlist2[p])
                plist.append(pointlist1[p])
                templist.reverse()
                #東端と西端が一致しないまま緯度の上限に達してしまった場合の処理です。plistに西端、東端の座標を順に付加し、以下の分岐ではうえの"if"の際と同じ操作をし、roopから抜け出します
                if ori==0:
                    p2list=plist+templist
                    pointlist.append(p2list)
                    break
                else:
                    templist.append(ori)
                    p2list=plist+templist
                    pointlist.append(p2list)
                    break        
            else:
                templist.append(pointlist1[p])
                plist.append(pointlist2[p])
                p　+= 1
                #西端と東端が一致しないときの処理です。plistに西端、templistに東端の座標を付加します。
    output ='/~+str(pref)+'.shp'
    #完成したレイヤを描きだすパスです
    layerfields = QgsFields()
    #レイヤの属性を作っていきます
    layerfields.append(QgsField('id', QVariant.Int))
    #idフィールド（整数）を作成
    layerfields.append(QgsField('name', QVariant.String))
    #nameフィールド（文字列）を作成
    writer = QgsVectorFileWriter(output, 'UTF-8', layerfields, QgsWkbTypes.Polygon, QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')
    feat = QgsFeature()
    for pl in pointlist:
        feat.setGeometry(QgsGeometry.fromPolygonXY())
        #pointlistには2重リストで、ポリゴンの頂点が格納されているのでこれを使ってマルチパートでジオメトリを作成します。
        feat.setAttributes([pref, code[pref]])
        #属性のidフィールドに都道府県コード、nameフィールドに都道府県名を入れます。
        writer.addFeature(feat)
    iface.addVectorLayer(output, '', 'ogr')
    #画面に出力します