基盤地図情報 数値標高を読み込んでみる
応募用に考えていたんですが
何かネタ出したいなと思い考えていたのですが、いまいちパッとしません。
でも中途半端ですが折角作ったのでここにさらしておきます。
結論から言うと、"コマンドラインでやった方が早いよ"的なネタです。
処理手順
基盤地図情報を使うことが望ましいので、QGISで基盤地図情報 数値標高を読み込もうと思っていました。
手順としては、前回紹介したScript Runnerプラグインを使用して、
- JPGIS形式そのままでは読めませんので、GDALで読めるフォーマットのうち簡単なものに変換(ESRI ASCII Raster formatにしました)
- 変換したものをGDALを使用して開く
- ついでなので、GDALで等高線のshapefileを作成して、それも開く
としました。
中間ファイルとして
が開いたファイルと同じディレクトリに作成されちゃいます。
ソースです
もし試してみたい方は下記からどうぞ。
dem.py 
Script Runnerプラグインのインストールもお忘れなく。
ソースの中身はこんな感じです。
特化した処理になっていますので、あしからず。
すいません、表示が一部壊れていますが直し方がわからないのでそのまま置いておきます
直せました。スーパーpre記法にすればいいのね。
# -*- coding: utf-8 -*-
import os, sys
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
class demHandle():
def __init__(self, filename):
self.row = 0
self.col = 0
self.west = 0
self.east = 0
self.north = 0
self.south = 0
self.alti = []
self._readFile(filename)
def _readFile(self, filename):
file = open(filename, 'r')
for line in file:
if line.startswith('<jps:westBoundLongitude>'):
self.west = float(self._getValue(line))
self.east = float(self._getValue(file.next()))
self.south = float(self._getValue(file.next()))
self.north = float(self._getValue(file.next()))
break
for line in file:
if line.startswith('<jps:high>'):
val = (self._getValue(file.next())).split(' ')
self.col = int(val[0])+1
self.row = int(val[1])+1
break
for line in file:
if line.startswith('<alti>'):
self.alti.append(self._getValue(line))
for i in range(5):
file.next()
file.close()
def _getValue(self, line):
return ((line.split('>'))[1].split('<'))[0]
def getAlti(self):
return self.alti.pop(0)
def getAltiRow(self, row):
return self.alti[row*self.col:(row+1)*self.col]
def run_script(iface):
demfile = QFileDialog.getOpenFileName(iface.mainWindow(), ur'基盤地図情報 数値標高(JPGIS形式)', '', 'XML files (*.xml)')
dem = demHandle(demfile)
#ESRI ASCII Raster
xyz = demfile.left(demfile.length()-4)+r'.xyz'
file = open(xyz, 'w')
outbuff = []
outbuff.append(' '.join(('ncols','%s\n' % dem.col)))
outbuff.append(' '.join(('nrows','%s\n' % dem.row)))
outbuff.append(' '.join(('xllcorner','%s\n' % dem.west)))
outbuff.append(' '.join(('yllcorner','%s\n' % dem.south)))
outbuff.append(' '.join(('cellsize','%.10f\n' % ((dem.east-dem.west)/dem.col))))
outbuff.append(' '.join(('NODATA_value','-9999\n')))
for x in range(dem.row):
outbuff.append(''.join((' '.join(dem.getAltiRow(x)), '\n')))
file.writelines(outbuff)
file.close()
#prj
prj = demfile.left(demfile.length()-4)+r'.prj'
file = open(prj, 'w')
outbuff = [r'GEOGCS["GCS_JGD_2000",DATUM["D_JGD_2000",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]']
file.writelines(outbuff)
file.close()
#contour shapefile
shp = demfile.left(demfile.length()-4)+r'.shp'
process = QProcess()
process.start(r'gdal_contour -i 10.0 %s %s' % (xyz, shp), QIODevice.ReadOnly)
demLayer = iface.addRasterLayer(xyz,xyz)
demLayer.setDrawingStyle(QgsRasterLayer.SingleBandPseudoColor)
demLayer.setColorShadingAlgorithm(QgsRasterLayer.PseudoColorShader)
contLayer = iface.addVectorLayer(shp,shp,'ogr')
iface.mapCanvas().refresh()
役に立つのかな？
どんな処理をまとめておけば役に立つのかさっぱりわかってません。
QGIS上でやる利点を考えると、やっぱり印刷系なのかなぁという気もしますし。
時間見つけてネタ思いついたら、また何か作ってみます。
