基盤地図情報数値標高を読み込んでみる(修正版)

2012/3/5の記事の修正版です

Pythonで強引に基盤地図情報数値標高をQGISに取り込むぞというのを前に書いたのですが、中途半端な対応で開けないファイルがありました。開けないのは悲しいので、その辺だけ修正しておきます。

もう一回言っときますが、コマンドラインでやれよというような内容です。

処理手順

手順は同じ。
ここで紹介したScript Runnerプラグインを使用して、

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

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

ESRI ASCII Raster formatのDEM
等高線のshapefile

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

直したところ

拡張子は".asc"ですね。もしくは".txt"。
xx yy への対応がなかったので追加。

といったところです。

実行イメージ

前回と同じなので割愛

ソースです

dem20120330.py
Script Runnerプラグインのインストールもお忘れなく。
ソースの中身はこんな感じです。

# -*- coding: utf-8 -*-

import os, sys
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *

class demHandler():
  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.lstrip().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.lstrip().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.lstrip().startswith('</jps:sequencingRule>'):
        break
      if line.lstrip().startswith('<alti>'):
        self.alti.append(self._getValue(line))
        for i in range(5):
          file.next()
    
    for line in file:
      if line.lstrip().startswith('<jps:coordValues>'):
        val = self._getValue(line).split(' ')
        x = int(val[0])
        y = int(val[1])
        break
        
    file.close()
    
    for i in range(y*self.col+x):
      self.alti.insert(0, '-9999.')
    
    for i in range((self.col*self.row)-len(self.alti)):
      self.alti.append('-9999.')
    
  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]

class ascWriter():
  def __init__(self, dem):
    self.dem = dem
  
  def write(self, filename):
    #ESRI ASCII Raster
    file = open(filename, 'w')
    outbuff = []
    outbuff.append(' '.join(('ncols','%s\n' % self.dem.col)))
    outbuff.append(' '.join(('nrows','%s\n' % self.dem.row)))
    outbuff.append(' '.join(('xllcorner','%s\n' % self.dem.west)))
    outbuff.append(' '.join(('yllcorner','%s\n' % self.dem.south)))
    outbuff.append(' '.join(('cellsize','%.10f\n' % ((self.dem.east-self.dem.west)/self.dem.col))))
    outbuff.append(' '.join(('NODATA_value','-9999.\n')))
    for x in range(self.dem.row):
      outbuff.append(''.join((' '.join(self.dem.getAltiRow(x)), '\n')))
    file.writelines(outbuff)
    file.close()
    
    #prj
    fileInfo = QFileInfo(filename)
    prjfile = r'%s/%s.prj' % (fileInfo.path(), fileInfo.baseName())
    file = open(prjfile, '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()
    
def run_script(iface):
  demfile =  QFileDialog.getOpenFileName(iface.mainWindow(), ur'基盤地図情報 数値標高(JPGIS形式)', '', 'XML files (*.xml)')
  dem = demHandler(demfile)
  fileInfo = QFileInfo(demfile)
  ascfile = r'%s/%s.asc' % (fileInfo.path(), fileInfo.baseName())
  asc = ascWriter(dem)
  asc.write(ascfile)

  #contour shapefile
  shpfile = r'%s/%s.shp' % (fileInfo.path(), fileInfo.baseName())
  process = QProcess()
  process.start(r'gdal_contour -i 10.0 %s %s' % (ascfile, shpfile), QIODevice.ReadOnly)
  process.waitForFinished()
  
  demLayer = iface.addRasterLayer(ascfile,ascfile)
  demLayer.setDrawingStyle(QgsRasterLayer.SingleBandPseudoColor)
  demLayer.setColorShadingAlgorithm(QgsRasterLayer.PseudoColorShader)
  contLayer = iface.addVectorLayer(shpfile,shpfile,'ogr')
  iface.mapCanvas().refresh()