はじめに

皆さん使ってますか？Natural Earth

（個人）気候変化情報研究室では、主に海岸線データを使わせていただいています。
以下の記事では、様々な投影法で日本を中心とした世界地図を描く方法を紹介しました。

おっと、Natural Earthから入手した海岸線データを描こうとすると、ところどころに横線が・・・

これは、地球を円筒に投影して筒を開くとき、鋏を入れる線上にある地物がきちんと切れていないことが原因らしい。

なら、切れ目に位置する地物をきちんと切ってあげればよい。
ということで、GeoPandasを使ってやってみましょう。

世界地図ベクタデータの入手

たくさんの地図データをpublic domainで提供している、Natural Earthから入手します。

今回は、1:50m解像度のCoastlineとLandをダウンロードしました。

折り返し部分に切れ目を入れる

東経135度を中心にした世界地図を作る場合、折り返しは西経45度になりますので、GeoPandasを使って切れ目を入れます。

まずは海岸線データから。

Python

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString, Polygon

# 処理するシェープファイルを読み込む
os.chdir("/mnt/d/Users/kazuh/natural_earth/")
src_gdf = gpd.read_file("./physical/ne_50m_coastline.shp")

# ジオメトリが西経-45度に引っかかる地物を抽出
intersects_gdf = src_gdf.loc[
    src_gdf.geometry.intersects(
        LineString([(-45, -90), (-45, 90)])
    )
]

# 西経45度で切った西側部分の地物
W_split_gdf = gpd.GeoDataFrame(
    intersects_gdf.drop("geometry",axis=1),
    geometry=intersects_gdf.geometry.intersection(
        Polygon([(-180,-90),(-45, -90),(-45, 90),
                (-180,90),(-180,-90)])
    )
)

# 西経45度で切った東側部分の地物
E_split_gdf = gpd.GeoDataFrame(
    intersects_gdf.drop("geometry",axis=1),
    geometry=intersects_gdf.geometry.intersection(
        Polygon([(180,-90),(-45, -90),(-45, 90),
                (180,90),(180,-90)])
    )
)

#　西経45度に引っかからない地物と統合
dst_gdf = pd.concat([
        src_gdf.drop(intersects_gdf.index),
        W_split_gdf,
        E_split_gdf
        ],
        ignore_index=True).reset_index(drop=True)

# 西経45度より西側の地物を抽出
W_part_gdf = dst_gdf.cx[:-45.01,:]

# 地球一回転分（360度）移動
dst_gdf.loc[W_part_gdf.index,"geometry"] = W_part_gdf.geometry.translate(360,0)

# 投影法を更新（180度以上を処理できるようにする）
dst_gdf.crs = "+proj=longlat +datum=WGS84 +over +no_defs"

# geopackageに保存
dst_gdf.to_file("JP_natural_earth.gpkg",
            layer="ne_50m_coastline",
            driver="GPKG",
            overwrite=True)

陸地データも加工します。
折り返しの西経45度に切れ目を入れ、１８０度のラインで西と東に分かれている地物を融合します。

Python

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely import segmentize
from shapely.ops import polygonize
from shapely.geometry import LineString, Polygon

# 処理するシェープファイルを読み込む
os.chdir("/mnt/d/Users/kazuh/natural_earth/")
src_gdf = gpd.read_file("./physical/ne_50m_land.shp")

# ジオメトリが西経45度に引っかかる地物を抽出
intersects_gdf = src_gdf.loc[
    src_gdf.geometry.intersects(
        LineString([(-45, -90), (-45, 90)])
    )]

# 西経45度を境に西側と東側のエリアの境界線ポリゴンを作成する
C_line = segmentize(
            LineString([(-45,-90),(-45,90)]),
            max_segment_length=1
        )

W_bbox = list(polygonize([
    C_line,
    LineString([(-45,90),(-180,90),(-180,-90),(-45,-90)])
    ]))[0]

E_bbox = list(polygonize([
    C_line,
    LineString([(-45,90),(180,90),(180,-90),(-45,-90)])
    ]))[0]

# 西経45度で切った西側部分の地物
W_split_gdf = gpd.GeoDataFrame(
    intersects_gdf.drop("geometry",axis=1),
    geometry=intersects_gdf.geometry.intersection(W_bbox)
    )

# 西経45度で切った東側部分の地物
E_split_gdf = gpd.GeoDataFrame(
    intersects_gdf.drop("geometry",axis=1),
    geometry=intersects_gdf.geometry.intersection(E_bbox)
    )

# 西経45度に引っかからない地物と統合
src2_gdf = pd.concat([
    src_gdf.drop(intersects_gdf.index),
    W_split_gdf,
    E_split_gdf
    ],
    ignore_index=True).reset_index(drop=True)
 
# ジオメトリが西経180度に引っかかる地物を抽出
W_intersects_gdf = src2_gdf.loc[
    src2_gdf.geometry.intersects(
        LineString([(-180, -90), (-180, 90)])
    )]

# 西側の地物を地球一回転分（360度）移動させる
W_intersects_gdf = gpd.GeoDataFrame(
    W_intersects_gdf.drop("geometry",axis=1),
    geometry=W_intersects_gdf.geometry.translate(360,0)
    )

# ジオメトリが東経180度に引っかかる地物を抽出
E_intersects_gdf = src2_gdf.loc[
    src2_gdf.geometry.intersects(
        LineString([(180, -90), (180, 90)])
    )]

# 接続関係を整理
sjoin_gdf = E_intersects_gdf.sjoin(W_intersects_gdf,
    how="inner",predicate="touches")
sjoin_gdf.index.name = "index_left"
sjoin_pair = sjoin_gdf["index_right"].reset_index()

# ジオメトリを結合する
E_geom = E_intersects_gdf.geometry.loc[sjoin_pair["index_left"]]
W_geom = W_intersects_gdf.geometry.loc[sjoin_pair["index_right"]]
union_geom = E_geom.union(W_geom,align=False)

# 結合した地物のジオメトリを更新
dst_gdf = src2_gdf.drop(sjoin_pair["index_right"])
dst_gdf.loc[union_geom.index,"geometry"] = union_geom

# 西経45度より西側の地物を抽出
W_part_gdf = dst_gdf.cx[:-45.01,:]

# 地球一回転分（360度）移動
dst_gdf.loc[W_part_gdf.index,"geometry"] = W_part_gdf.geometry.translate(360,0)

# 投影法を更新（180度以上を処理できるようにする）
dst_gdf.crs = "+proj=longlat +datum=WGS84 +over +no_defs"

# geopackageに保存
dst_gdf.to_file("JP_natural_earth.gpkg",
    layer="ne_50m_land",
    driver="GPKG",
    overwrite=True)

背景をぬりつぶすための境界線も作成しておきます。

Python

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely import segmentize
from shapely.ops import polygonize
from shapely.geometry import LineString, Polygon

# 東西南北のラインを作成
E_line = segmentize(
            LineString([(315,-90),(315,90)]),
            max_segment_length=1
        )
N_line = segmentize(
            LineString([(315,90),(-45,90)]),
            max_segment_length=30
        )
W_line = segmentize(
            LineString([(-45,90),(-45,-90)]),
            max_segment_length=1
        )
S_line = segmentize(
            LineString([(-45,-90), (315,-90)]),
            max_segment_length=30
        )

# ポリゴンに変換
JP_bbox = polygonize([
    E_line,N_line,W_line,S_line
    ])

# GeoDataFrameを作成
JP_bbox_gdf = gpd.GeoDataFrame(
                index=[0],
                geometry=list(JP_bbox),
                crs="+proj=longlat +datum=WGS84 +over +no_defs"
                )

# geopackageに保存
JP_bbox_gdf.to_file(
        "JP_natural_earth.gpkg",
        layer="bbox",
        driver="GPKG",
        overwrite=True
        )

陰影起伏図を東経135度を中心としたデータに作り直す。

入手したラスタはグリニッジを中心としたデータなので、東経135度を中心としたデータに作り直します。

Python

import os
import numpy as np
from osgeo import gdal, gdalconst, gdal_array 
from osgeo import osr

# ラスタを読み込む
os.chdir("/mnt/c/Users/hoge/natural_earth")
relief_ds = gdal.Open("./MSR_50M/MSR_50M.tif",gdalconst.GA_ReadOnly)

# データ配列を取り出す
relief_arr = relief_ds.GetRasterBand(1).ReadAsArray()

# 折り返す位置を特定
x_break = int(10800*(180-45)/360)

# データの折り返し
JP_relief_arr =np.concatenate([relief_arr[:,x_break:],relief_arr[:,:x_break]],axis=1)

# メモリ上にGdalDataSetを作成
x_size = relief_ds.RasterXSize
y_size = relief_ds.RasterYSize

JP_relief = gdal.GetDriverByName('MEM').Create("",x_size,y_size,1,gdal.GDT_Byte)

# 投影パラメータを設定
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromProj4("+proj=longlat +datum=WGS84 +over +no_defs")

JP_relief.SetProjection(outRasterSRS.ExportToWkt())
JP_relief.SetGeoTransform((-45.,0.03333333333333, 0,90.0, 0, -0.03333333333333))

#　海岸線データを収録したgpkgファイルに追加
output = gdal.GetDriverByName('GPKG').CreateCopy("JP_natural_earth.gpkg",JP_relief,options=["RASTER_TABLE=MSR_50M","APPEND_SUBDATASET=YES"])

outband = output.GetRasterBand(1)
outband.WriteArray(JP_relief_arr)

outband.FlushCache()
output =None

QGISに読み込んでみる

WGS84では普通に読み込めます。

日本を通る経度を中心にした地図については以下を参考にさせていただいています。

ロビンソン図法

proj.4 string

+proj=robin +lon_0=135 +x_0=0 +y_0=0 +datum=WGS84 +units=m +over +no_defs

モルワイデ図法

proj.4 string

+proj=moll +lon_0=135 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +over

サンソン図法

proj.4 string

+proj=sinu +lon_0=135 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +over

様々な投影法で、東経１３５度を中心にした世界地図を描くことができました。

　
（個人）気候変化情報研究室では、できるだけ簡単に、将来、地域の気候がどのように変化するかを示すグラフ・地図を作成する方法を開発し、紹介しています。
興味のある方はぜひ読んでみてくださいね～。