下準備
地図が3つのエリアに分かれているので、それぞれのエリアについて座標と緯度経度の関係を調査します
<本州・四国・九州>
x y long lat
青森大間崎 (548, 5) = 140.91, 41.55
鹿児島佐多岬 ( 62, 621) = 130.66, 31.00
差 (486, 616) = 10.25, 10.55
scalex = 10.25/486 = 0.021
scaley = 10.55/616 = 0.017
long = 130.66 + (x - 62) * 0.021
lat = 31.00 + (621 - y) * 0.017
-------------------
<北海道>
x y long, lat
白神岬 ( 28, 272) = 140.18, 41.40
long = 140.18 + ( x - 28 ) * 0.021
lat = 41.40 + (272 - y) * 0.017
-------------------
<沖縄>
x y long, lat
喜屋武岬 (533, 578) = 127.67, 26.08
long = 127.67 + ( x - 533 ) * 0.021
lat = 26.08 + (578 - y) * 0.017
プログラム
内外判定はsako氏のものを使わせていただいた。
%config InlineBackend.figure_formats = {'png', 'retina'}
%matplotlib
import matplotlib.pyplot as plt
from japanmap import picture, get_data, pref_points
import numpy as np#http://sak12.blogspot.com/2013/09/blog-post.html
def inpolygon(sx, sy, x, y):
'''
x[:], y[:]: polygon
sx, sy: point
'''
np = len(x)
inside = False
for i1 in range(np):
i2 = (i1+1)%np
if min(x[i1], x[i2]) < sx < max(x[i1], x[i2]):
if (y[i1] + (y[i2]-y[i1])/(x[i2]-x[i1])*(sx-x[i1]) - sy) > 0:
inside = not inside
return insidedef x2long(x):
return 130.66 + (x - 62) * 0.021
def y2lat(y):
return 31.00 + (621 - y) * 0.017
def x2longhokkaido(x):
return 140.18 + (x - 28) * 0.021def y2lathokkaido(y):
return 41.40 + (272 - y) * 0.017def x2longokinawa(x):
return 127.67 + (x -533) * 0.021def y2latokinawa(y):
return 26.08 + (578 - y) * 0.017def in_hokkaido(x, y):
pnt = pnts[0] # hokkaido
longs = [row[0] for row in pnt]
lats = [row[1] for row in pnt]
long = x2longhokkaido(x)
lat = y2lathokkaido(y)
if inpolygon(long, lat, longs, lats):
return 1 # hokkaido
else:
return 0def in_okinawa(x, y):
pnt = pnts[46] # okinawa
longs = [row[0] for row in pnt]
lats = [row[1] for row in pnt]
long = x2longokinawa(x)
lat = y2latokinawa(y)
if inpolygon(long, lat, longs, lats):
return 47 # okinawa
else:
return 0def in_ken(x, y):
long = x2long(x)
lat = y2lat(y)
for i in range(1, 47): #1, ... 46
pnt = pnts[i]
longs = [row[0] for row in pnt]
lats = [row[1] for row in pnt]
if inpolygon(long, lat, longs, lats):
return i + 1
return 0
plt.rcParams['figure.figsize'] = 6, 6
plt.imshow(picture());
qpqo = get_data()
pnts = pref_points(qpqo) #県の緯度経度データ
ken = 0try:
while True:
a=plt.ginput(n=1)[0] #クリックして座標を取得
x = a[0]
y = a[1]
if x < 314 and y < 284: # hokkaido_area
ken = in_hokkaido(x, y)
elif x > 508 and y > 520: # okinawa_area
ken = in_okinawa(x, y)
else: # honshu and shikoku area
ken = in_ken(x, y)
if ken > 0:
plt.imshow(picture({ken: 'blue'}))
except Exception as e:
print(e)
コメント
最新を表示する
NG表示方式
NGID一覧