japanmapでクリックした県に色を塗る

ページ名:japanmapでクリックした県に色を塗る

下準備

地図が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 inside

def 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.021

def y2lathokkaido(y):
    return 41.40 + (272 - y) * 0.017

def x2longokinawa(x):
    return 127.67 + (x -533) * 0.021

def y2latokinawa(y):
    return 26.08 + (578 - y) * 0.017

def 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 0

def 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 0   

def 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 = 0

try:
  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)

 

シェアボタン: このページをSNSに投稿するのに便利です。

コメント

返信元返信をやめる

※ 悪質なユーザーの書き込みは制限します。

最新を表示する

NG表示方式

NGID一覧