読者です 読者をやめる 読者になる 読者になる

waigani's diary

QGISを中心にFOSS4Gをいじくる

基盤地図情報 数値標高を読み込んでみる

応募用に考えていたんですが

何かネタ出したいなと思い考えていたのですが、いまいちパッとしません。
でも中途半端ですが折角作ったのでここにさらしておきます。
結論から言うと、"コマンドラインでやった方が早いよ"的なネタです。

処理手順

基盤地図情報を使うことが望ましいので、QGIS基盤地図情報 数値標高を読み込もうと思っていました。
手順としては、前回紹介したScript Runnerプラグインを使用して、

  • JPGIS形式そのままでは読めませんので、GDALで読めるフォーマットのうち簡単なものに変換(ESRI ASCII Raster formatにしました)
  • 変換したものをGDALを使用して開く
  • ついでなので、GDALで等高線のshapefileを作成して、それも開く

としました。
中間ファイルとして

が開いたファイルと同じディレクトリに作成されちゃいます。

実行イメージ

  • 読み込むJPGIS形式の数値標高データを選択

  • DEMが読み込まれます、おまけで等高線付き


ソースです

もし試してみたい方は下記からどうぞ。
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上でやる利点を考えると、やっぱり印刷系なのかなぁという気もしますし。
時間見つけてネタ思いついたら、また何か作ってみます。

追記 2012.3.29

見返してみるとやっつけ感が酷いですね。。

  • 拡張子は".asc"が正しいです。".xyz"で書こうとしていた時の名残が。
  • xx yy への対応が一切無いです。NODATAがデータの最初に続く場合に記述を省略するためのタグなのですが、対応が無いのでファイルによってはおかしくなります。

少なくとも基盤地図情報でダウンロード出来るデータには対応するように近々書き換えますm(__)m