Lucy’s blog

プログラムで気になったことの覚書

Vimで改行コード^Mを\nに置換

Windowsから送られてきたASCII形式のファイルをVimで見ると、
改行コードが全て^Mとなっていたので、
Macの改行コード\nに置換します。

Vimでファイルを開いて以下を普通に入力
:%s/^M/\n/g

すると、エラーが返ってきました。
E486: Pattern not found: ^M

^Mが見つからないらしいです。

調べたところ、この^Mを入力するには以下のようにしないといけないみたいです。
理由はよくわかりません。ご存じの方どなたか教えてください。

^Mの部分を、Ctrl+V、Ctrl+M で入力

:%s/^M/\n/g
これで置換するとうまく改行できました。

dbfファイルをcsvに変換

国土数値情報からDLしたシェープファイル形式のデータの中にある、

dbfファイルのデータを使いたくて、

これを一気に(全国分だと170ファイルくらい)csvに変換したいときのシェルスクリプトbash)。

少数ファイルだけの処理の場合はExcelでファイル形式をcsvに変更して保存をすれば事足ります。

 

あらかじめ全てのファイル名(今回は1次メッシュ番号)を記載したテキストファイルを作成しておきます。

$head meshNum.txt

3035
3036
3622
3623
3624
3631
3641
3653
3724
3725

 

環境:OSX Sierra

 #!/bin/bash
meshNum="meshNum.txt"
mesh=(`cat ${meshNum} | awk '{print $1}'`)


for meshCode in ${mesh[@]}
  do
    echo ${meshCode};

    #入力ファイル名を取得

    bfile="../DL/L03b-91_${meshCode}_LandUseSubdivisionMesh.dbf";
 
    if [ -e ${bfile} ] ; then  #ファイルが存在するとき実行

  #文字列があって邪魔なときは sed '/mojiretsu/d'で削除

        strings ${bfile} | awk '{for(i=1;i<=NF;i++)print substr($i,1,10)","substr($i,11,14)}' | sort -k1 | sed '/L03/d' > ${meshCode}_tky.csv

    fi
done
echo "done"

 

100mメッシュを緯度経度に変換する

今回利用したデータはこちら

土地利用細分メッシュデータ

 

参考にしたのはこちらの記事たちです。

標準メッシュの体系とコード

white-bear.info

1/10(100m)メッシュコードを緯度経度に変換する関数をつくります。

環境:OSXSierra、Anaconda、python3.6

 ============

import pandas as pd
import numpy as np

def CoordinateLatLon( MeshCode ):
    #変数を初期化

    LatArr =
    LonArr =

    for i in MeshCode:

  #文字列型に変換
        Mesh = str(i)
        AB  = Mesh[0:2]
        CD  = Mesh[2:4]
        E   = Mesh[4]
        F   = Mesh[5]
        G   = Mesh[6]
        H   = Mesh[7]
        I   = Mesh[8]
        J   = Mesh[9]
        AB  = float(AB)
        CD  = float(CD)
        E   = float(E)
        F   = float(F)
        G   = float(G)
        H   = float(H)
        I   = float(I)
        J   = float(J)
        Lat = (AB/1.5*3600+E*5*60+G*30+I*3)/3600
        Lon = ((CD+100)*3600+F*7.5*60+H*45+J*4.5)/3600
        LatArr.append(Lat)
        LonArr.append(Lon)
    return (LatArr,LonArr)

 

#関数を利用する。

#データの読み込み(pandasを使っています)

df  = pd.read_csv('5339.csv',header=None,names=['MeshCode'])

MeshCode = np.array(df.MeshCode)

print(df.MeshCode)

>[5339000000 5339000001 ...]

Lat,Lon = CoordinateLatLon(MeshCode)

緯度経度を用いてmatplotlibで地図をplotする

ざっくりですが覚書です。

環境:OSX,Anaconda,Python3.5

#必要なものをインポート

#matplotlibとBasemap (OSXではどちらもpip installで入手可能だった気がします。違ったらスミマセン。。)

 

 

import pandas as pd

import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
from mpl_toolkits.basemap import Basemap,cm

 

#データを用意

#csvファイルをpandas dataframeで読み込む

#(今回は1/10メッシュコードの土地利用データを緯度経度変換したもの)

df = pd.read_csv('test.csv',header=0)

----------

Lat,Lon,MeshCode_100m,Landuse

35.33333333333333,139.0,5339000000,100

35.33333333333333,139.00125,5339000001,1400

...

----------

 

#図の大きさ等指定

fig = plt.figure(figsize=(10,5)) #図の大きさ

fig.patch.set_alpha(0.0) #図全体の背景透明度

 

#描画領域の指定

north = max(df.Lat)
south = min(df.Lat)
east = max(df.Lon)
west = min(df.Lon)

m=Basemap(projection='merc',llcrnrlat=south,urcrnrlat=north,llcrnrlon=west,urcrnrlon=east,resolution='l')
m.drawparallels(np.arange(24.0, 48.1, 2.0), labels = [1,0,0,0], fontsize=12) #緯線を引く
m.drawmeridians(np.arange(125.0, 145.1, 4.0), labels = [0,0,0,1], fontsize=12) #経線を引く

 

#散布図を描く準備。流れは以下のような感じ。

##########################

#lonlist=list(df.Lon)
#latlist=list(df.Lat)
#data = list(df.data)
#x, y = m(lonlist,latlist)
#nx = len(lonlist)
#ny = len(latlist)
#lons,lats = m.makegrid(400,400)

###########################

 

#データフレームから必要なデータを取得

df0 = df[(df.Landuse == 100)]
df1 = df[(df.Landuse == 1400)]

 

#緯度経度をリスト化

lonall=list(df.Lon)
latall=list(df.Lat)

lon0=list(df0.Lon)
lat0=list(df0.Lat)
lon1=list(df1.Lon)
lat1=list(df1.Lat)

 

#プロットしたいデータの緯度経度を指定

xall,yall = m(lonall,latall)

x0, y0 = m(lon0,lat0)

x1, y1 = m(lon1,lat1)

 

#色を設定

ccolor=['#008000','#00FF00'] #16進数カラーコードで指定してもよし、色名で指定してもよし。

#ラベルを設定
cname= ['田','海浜']

 

#マーカーのサイズを指定

msize=8.0

#マーカーの透明度を指定

alpha=0.5

 

#散布図をプロット

plt.scatter(xall,yall,c='0.85',alpha=1.0,s=2.0,label='九十九里町',marker='o')

plt.scatter(x0,y0,c=ccolor[0],alpha=alpha,s=msize,label=cname[0],marker='o')#100
plt.scatter(x1,y1,c='olive',alpha=alpha,s=msize,label=cname[1],marker='o')#1400

 

#枠線をグレーにする

ax = plt.gca()                        # get current axis
ax.spines["right"].set_color("lightgray")  # 右枠
ax.spines["top"].set_color("lightgray")    # 上枠
ax.spines["left"].set_color("lightgray")      # 左枠
ax.spines["bottom"].set_color("lightgray")    # 下枠
ax.patch.set_alpha(0.0) #描画領域の中の背景透明度

 

#フォントを指定

fp = FontProperties(fname='/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc')

 

#タイトルを指定
plt.title('九十九里町',fontdict={'fontproperties':fp})

 

#凡例の設定

leg=plt.legend(loc='upper left',prop=fp,markerscale=1.0,bbox_to_anchor(-0.5,0.8),fancybox=1)
leg.get_frame().set_alpha(0.0) #凡例を透明にする←これやったら透明にはなったけどlegendの枠線が消えてしまいました。。

 

#ファイルを保存する
plt.savefig('example.png')

 

#描画を表示する

plt.show()

 

だいたいこんな感じの図が描けます。

f:id:harutaromaru:20170712103949p:plain

 

 

測地系のはなし

GoogleMapからある施設の緯度経度を求めて、Pythonでplotするという事案がありました。

 

そこでGoogleMapの測地系を調べてみたところ、WGS84世界測地系)となっていました。

WGS84とは何か、JGD2000と違うのか、と思い調べるとこんな記事がありました。

4.測地系 - QGIS入門

2~3メートルの誤差が許容される程度のものであれば、
いつのWGS84でもJGD2000同じと考えてかまいません。

 

ちなみにJGD2011というのは東日本大震災のあとに生じたズレを加えたものだそうです。こちらも1〜2mとかのズレを気にしない場合はJGD2000と同じとしてよさそうです。

 

GoogleMapで緯度経度情報を出す

1.GoogleMapに目当ての施設を入力(例:東京駅)

f:id:harutaromaru:20170711091038p:plain

2.赤いマークを右クリックして、「この場所について」をクリック

f:id:harutaromaru:20170711091313p:plain

3.下の方に緯度経度情報が出てきます。

 

日本測地系へ変換

 今回の地図は日本測地系(Tokyo)のデータを使っているので、

いいか悪いかは謎ですが、まあざっくり位置が合えばいいので、以下のサイトで 世界測地系日本測地系 に変換します。

Web版TKY2JGD - 国土地理院

 

東京駅

緯度:35.681259 → 35.678034633

軽度:139.767062 → 139.770285697

 

次回はこれをPythonでplotしていきます。