SRTMデータの可視化

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47

import sys
 
nrows = ncols = 1201
xdim = ydim = 0.000833
 
args = sys.argv
filename = args[1]
name, ext = os.path.splitext(filename)
outname = name + '.csv'
 
longitude = float(name[1:3])
latitude = float(name[4:7])
 
print('longitude: ' + name[0] + name[1:3])
print('latitude: ' + name[3] + name[4:7])
 
ulxmap = longitude + (xdim/2.0)
ulymap = (latitude+1.0) - (ydim/2.0)
 
a = array.array('h')
f = open(filename, 'rb')
 
a.fromfile(f, nrows*ncols)
f.close()
 
a.byteswap()
alist = a.tolist()
 
out = open(outname, 'w')
 
count = 0
for row in range(nrows):
  for col in range(ncols):
    x = ulxmap + xdim * col
    y = ulymap - ydim * row
    z = max(0.0, alist[count])
    out.write(str(x) + ',' + str(y) + ',' + str(z/10000.0) + '\n')
    count += 1
 
    if count%100000 == 0:
        print(str(int(100.0 * count/(nrows*ncols))) + '%')
 
print('100%')
print('node number: ' + str(count))
 
out.close()

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46

 
 
-
!
-
!
 
 
 
 
 
 
 
-
!
 
 
 
 
 
 
-
!
 
 
 
 
 
-
!
 
 
 
 
 
 
 
 
 
-
!
 
 
 
 
 

import os.path
import sys
#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
 
args = sys.argv
filename = args[1]
 
workingDir = './'
fileName = args[1]
 
# create a new 'CSV Reader'
csvReader = CSVReader(FileName=[workingDir + fileName + '.csv'])
csvReader.DetectNumericColumns = 1
csvReader.UseStringDelimiter = 1
csvReader.HaveHeaders = 0
csvReader.FieldDelimiterCharacters = ','
csvReader.MergeConsecutiveDelimiters = 0
 
# create a new 'Table To Structured Grid'
tableToStructuredGrid = TableToStructuredGrid(Input=csvReader)
tableToStructuredGrid.WholeExtent = [0, 1200, 0, 1200, 0, 0]
tableToStructuredGrid.XColumn = 'Field 0'
tableToStructuredGrid.YColumn = 'Field 1'
tableToStructuredGrid.ZColumn = 'Field 2'
 
# create a new 'Calculator'
calculator = Calculator(Input=tableToStructuredGrid)
calculator.AttributeMode = 'Point Data'
calculator.CoordinateResults = 0
calculator.ResultNormals = 0
calculator.ResultTCoords = 0
calculator.ResultArrayName = 'z'
calculator.Function = 'coordsZ'
calculator.ReplaceInvalidResults = 1
calculator.ReplacementValue = 0.0
 
# save data
SaveData(workingDir + fileName + '.vts', proxy=calculator, 
           Writealltimestepsasfileseries=0, 
           DataMode='Appended', 
           EncodeAppendedData=1, 
           CompressorType='ZLib')