ラスターデータの外枠

こんな図を作ります

ググっても出てこなかったので作りました．↓

ラスターデータの外枠を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()