2地点間の緯度経度から直線距離でもっとも近いデータを取得する話 geodist関数

親戚の子供の宿題を手伝っていたんですが、以下の問題に頭抱えました。

「県庁所在地で、直線距離がもっとも近いのは何県と何県?
地図帳とものさしを使って考えてみよう」

うわぁ~、凄いめんどい。

「名古屋と岐阜であってるよね」って言われても、そんなのわかんない。
大阪-神戸だと思ったけど自信なし

あたりをつけて、さし当ててみるけど、結構難しいぞこれ。
わかります??

困った時のSAS On Demand
ネットさえつながればどこでも自由にSAS使えるので、ちょっとPC借りて
以下のコードをパチパチ。

data q1;
length no 8. shi $10. ido keido 8.;
input no shi ido keido;
cards;
1 札幌市 43.063968 141.347899
2 青森市 40.824623 140.740593
3 盛岡市 39.703531 141.152667
4 仙台市 38.268839 140.872103
5 秋田市 39.7186 140.102334
6 山形市 38.240437 140.363634
7 福島市 37.750299 140.467521
8 水戸市 36.341813 140.446793
9 宇都宮市 36.565725 139.883565
10 前橋市 36.391208 139.060156
11 さいたま市 35.857428 139.648933
12 千葉市 35.605058 140.123308
13 新宿区 35.689521 139.691704
14 横浜市 35.447753 139.642514
15 新潟市 37.902418 139.023221
16 富山市 36.69529 137.211338
17 金沢市 36.594682 136.625573
18 福井市 36.065219 136.221642
19 甲府市 35.664158 138.568449
20 長野市 36.651289 138.181224
21 岐阜市 35.391227 136.722291
22 静岡市 34.976978 138.383054
23 名古屋市 35.180188 136.906565
24 津市 34.730283 136.508591
25 大津市 35.004531 135.86859
26 京都市 35.021004 135.755608
27 大阪市 34.686316 135.519711
28 神戸市 34.691279 135.183025
29 奈良市 34.685333 135.832744
30 和歌山市 34.226034 135.167506
31 鳥取市 35.503869 134.237672
32 松江市 35.472297 133.050499
33 岡山市 34.661772 133.934675
34 広島市 34.39656 132.459622
35 山口市 34.186121 131.4705
36 徳島市 34.06577 134.559303
37 高松市 34.340149 134.043444
38 松山市 33.84166 132.765362
39 高知市 33.559705 133.53108
40 福岡市 33.606785 130.418314
41 佐賀市 33.249367 130.298822
42 長崎市 32.744839 129.873756
43 熊本市 32.789828 130.741667
44 大分市 33.238194 131.612591
45 宮崎市 31.91109 131.423855
46 鹿児島市 31.560148 130.557981
47 那覇市 26.212401 127.680932
;
run;

県庁所在地と緯度経度をネットから適当に拾ってきて
geodist関数にあてます。
geodist関数は(地点1の緯度,地点1の経度,地点2の緯度,地点2の経度)で、距離を返します(デフォルトなら単位はキロメートル、オプションでマイルにもできます)

ちなみに緯度経度には日本測地系と世界測地系があって、
微妙に違うので、できれば世界に変換して同じ測地系で比べましょう。

変換法はPHPだけど式は同じなので以下のページで紹介されているコードなんかがいいと
思います
http://d.hatena.ne.jp/nakamura001/20080501/1209660263


で、書いたコードは以下のとおり、これぐらいのオブザベーション数なら
直積作ってもよかったですけどハッシュでもわりとあっさり書けました。

data a1;
length _shi $10. _ido _keido 8.;
if _N_=1 then do;
call missing(_shi,_ido,_keido);
declare hash h1(dataset:'q1(rename=(shi=_shi ido=_ido keido=_keido))');
h1.definekey('no');
h1.definedata('_shi','_ido','_keido');
h1.definedone();
end;
set q1;
do no= 1 to 47;
rc=h1.find();
dist=geodist(ido,keido,_ido,_keido);
if dist^= 0 and (dist<mindist or mindist=.) then do;
mindist=dist;
moyori=_shi;
end;
end;
keep shi moyori mindist;
run;

proc sort data=A1;
by mindist;
run;
proc print;
run;

結果は



















































































滋賀の大津市と京都の京都市か~、確かに近いな!
10キロちょっとか

神戸にとっての最寄は大阪だけど大阪の最寄は奈良か。
へ~、そういう感じだったんですね。
勉強になりました。

2 件のコメント:

  1. SASYAMAさん、geodist関数面白いですね。
    そして和歌山の件も意外でした笑
    SQL直積でブログ書かせて頂きました。
    http://d.hatena.ne.jp/O_Kohsuke/20160204/1454592643

    返信削除
    返信
    1. O_Kohsukeさん

      コメント有難うございます。
      またSQLでの記事も拝見させていただきました。
      これぐらいであればSQLの方が楽でいいですね。

      和歌山市と直線で一番違いのが神戸というのは、本当に意外でした

      削除