こんな図を作ります
ググっても出てこなかったので作りました．↓
ラスターデータの外枠をcartopy上で描写
全球-180, 180, -90, 90の範囲で
解像度5分(赤道で9km)の0と1で構成されたラスターデータを考えます
x軸方向に4320
y軸方向に2160
グリッド存在する空間だと思ってください
1の部分は東京の都市域を表す都市マスクです
全球のうち138-141E, 34-37Nの範囲をcartopyで描写します
描写された空間内に1で示された都市マスクの外枠をハイライト
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
one = 138
two = 141
three = 34
four = 37
upperindex = (90-four)*12
lowerindex = (90-three)*12
leftindex = (180+one)*12
rightindex = (180+two)*12
img_extent = (one, two, three, four)
delta = 360/4320
x_left = img_extent[0]
y_top = img_extent[3]
raster_data = mask[upperindex:lowerindex, leftindex:rightindex]
# outer
upper_indices = []
lower_indices = []
left_indices = []
right_indices = []
rows, cols = raster_data.shape
for i in range(rows):
for j in range(cols):
if raster_data[i, j] == 1:
if raster_data[i-1, j] == 0:
upper_indices.append((i, j))
if raster_data[i+1, j] == 0:
lower_indices.append((i, j))
if raster_data[i, j-1] == 0:
left_indices.append((i, j))
if raster_data[i, j+1] == 0:
right_indices.append((i, j))
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(projection=ccrs.PlateCaree())
ax.coastlines()
ax.set_extent(img_extent)
ax.gridlines(draw_labels=True)
for up in upper_indices:
x_index = x_left + delta*(up[1])
y_index = y_top - delta*(up[0])
ax.plot([x_index, x_index+delta], [y_index, y_index], color='red', transform=ccrs.PlateCarree())
for lo in lower_indices:
x_index = x_left + delta*(lo[1])
y_index = y_top - delta*(lo[0]+1)
ax.plot([x_index, x_index+delta], [y_index, y_index], color='red', transform=ccrs.PlateCarree())
for le in left_indices:
x_index = x_left + delta*(le[1])
y_index = y_top - delta*(le[0]+1)
ax.plot([x_index, x_index], [y_index, y_index+delta], color='red', transform=ccrs.PlateCarree())
for ri in right_indices:
x_index = x_left + delta*(ri[1]+1)
y_index = y_top - delta*(ri[0]+1)
ax.plot([x_index, x_index], [y_index, y_index+delta], color='red', transform=ccrs.PlateCarree())
plt.show()
追記 (20260323)
LLM全盛期になったので上記のscriptをreviewしてもらった。
結果 > np.stackを使えば全球でも一瞬で境界線を取得できた(すごい)
- 横線をすべて取得
- 縦線をすべて取得
- 最後にそれらをconcatinateする
from __future__ import annotations
from pathlib import Path
import time
import numpy as np
# ============================================================
# Settings
# ============================================================
H08DIR = "/your/directory/name"
SUF = ".bin"
NX = 4320
NY = 2160
DX = 360.0 / NX # 5/60 degrees
DY = 180.0 / NY
L = NX * NY
OUT_DIR = Path(H08DIR) / "URB/out"
OUT_PATH = OUT_DIR / "boundary_segments.npy"
cty_msk_path = h08 / "map/dat/cty_msk_" / f"GPW_____00000000{SUF}"
P0MIS = np.float32(1.0e20)
# ============================================================
# Helper functions
# ============================================================
def to_2d_grid(data_1d, l2x, l2y):
grid = np.full((NY, NX), np.nan)
grid[l2y, l2x] = data_1d
return grid
# ============================================================
# Main
# ============================================================
def main():
h08 = Path(H08DIR)
# --- Load city mask ---
print("Loading city mask ...")
cty_raw = np.fromfile(str(cty_msk_path), dtype="<f4", count=L)
cty_id = np.where(
(cty_raw >= np.float32(P0MIS) * 0.9) | ~np.isfinite(cty_raw),
0,
cty_raw.astype(np.int32),
)
# --- Map to 2D grid ---
print("Mapping to 2D grid ...")
cty_float = np.where(cty_id > 0, cty_id.astype(np.float64), 0.0)
cty_2d = cty_float.reshape(NY, NX)
cty_2d = np.nan_to_num(cty_2d, nan=0.0)
delta = DX # = DY = 5/60 degrees
# Cell edge coordinates
lon_edge = -180.0 + np.arange(NX + 1) * delta
lat_edge = 90.0 - np.arange(NY + 1) * delta
# Vertical edges (between column x and x+1)
diff_v = (cty_2d[:, :-1] != cty_2d[:, 1:])
has_city_v = (cty_2d[:, :-1] > 0) | (cty_2d[:, 1:] > 0)
ys_v, xs_v = np.where(diff_v & has_city_v)
x_pos = lon_edge[xs_v + 1]
y_top = lat_edge[ys_v]
y_bot = lat_edge[ys_v + 1]
v_segs = np.stack([
np.column_stack([x_pos, y_top]),
np.column_stack([x_pos, y_bot]),
], axis=1)
# Horizontal edges (between row y and y+1)
diff_h = (cty_2d[:-1, :] != cty_2d[1:, :])
has_city_h = (cty_2d[:-1, :] > 0) | (cty_2d[1:, :] > 0)
ys_h, xs_h = np.where(diff_h & has_city_h)
x_left = lon_edge[xs_h]
x_right = lon_edge[xs_h + 1]
y_pos = lat_edge[ys_h + 1]
h_segs = np.stack([
np.column_stack([x_left, y_pos]),
np.column_stack([x_right, y_pos]),
], axis=1)
segments = np.concatenate([v_segs, h_segs], axis=0)
# --- Save ---
OUT_DIR.mkdir(parents=True, exist_ok=True)
np.save(str(OUT_PATH), segments)
print(f"\nSaved: {OUT_PATH}")
if __name__ == "__main__":
main()
下記画像はDelhiとその周辺の都市をzoom-upして表示した例です。
ちゃんと都市境界が赤線で表示されています。
描写はLineCollectionで線の集合として一気に描写できてしまいます。
cartopy上では下記のようにplotします。
def _overlay_boundary(ax, segments, transform,
linewidth=0.8, color="black"):
if segments is None or len(segments) == 0:
return
# Clip to axes extent for performance
extent = ax.get_extent(crs=transform) # (x0, x1, y0, y1)
x0, x1, y0, y1 = extent
seg_x = segments[:, :, 0] # (N, 2)
seg_y = segments[:, :, 1] # (N, 2)
in_extent = (
(seg_x.max(axis=1) >= x0) & (seg_x.min(axis=1) <= x1) &
(seg_y.max(axis=1) >= y0) & (seg_y.min(axis=1) <= y1)
)
visible = segments[in_extent]
if len(visible) == 0:
return
lc = LineCollection(visible, colors=color, linewidths=linewidth,
transform=transform, zorder=3)
ax.add_collection(lc)