十目十指: 電子国土ポータルに標高を載せる：標高検索

ExcelVBAでSRTM3を読込めるようにしたので、忘れないようにメモしておく。
（コメントは抜き）

２種類あるので注意。

（１）「ヤマレコ」さんのＰＨＰの焼直し版


'******************************************************************
'*　「ヤマレコ」さんのＰＨＰをＶＢＡに焼直したものです。
'*　　SRTMの生データではなく補正処理した結果が出力されます。
'******************************************************************

Function elevation(lat As Variant, lon As Variant) As Variant


  Dim H(1 To 4)   As Integer
  Dim I           As Integer
  Dim Cnt         As Variant

  Dim a1          As Byte
  Dim a2          As Byte

  Dim NV          As Variant
  Dim EV          As Variant
  Dim NI          As Integer
  Dim EI          As Integer

  Dim wDir        As String
  Dim fn          As String
  Dim ofs0        As Long
  Dim ofs         As Long

  Dim tmp         As Variant
  Dim La          As Variant
  Dim Lb          As Variant
  Dim Lc          As Variant
  Dim Ld          As Variant

  Dim ele1        As Variant
  Dim ele2        As Variant

  On Error GoTo errHandler

  NI = Int(lat)
  EI = Int(lon)
  wDir = ThisWorkbook.Path 'このBOOKの存在するディレクトリ
  fn = wDir & "\srtm\N" & NI & "E" & EI & ".hgt" 'そこに"SRTM"ディレクトリを作り、*hgt ファイルを入れる

  If Dir(fn) = "" Then
      MsgBox "file open error : 次のファイルが存在しません" & vbCr & fn
      Exit Function
  End If

  Open fn For Binary As #1
  NV = lat
  EV = lon

'計算式としては(1)が正しいように思えるが、「ヤマレコ」さんの計算結果を見ると
'(2)を使っている模様。
'ここでは、検証のため(2)を使うことにする。

  'ofs0 = (Fix((1 - (NV - NI)) * 1200) * 1201 + Fix((EV - EI) * 1200)) * 2 + 1 '-->(1)
  ofs0 = (Fix((1 - (NV - NI)) * 1200 - 1) * 1201 + Fix((EV - EI) * 1200)) * 2 + 1 '-->(2)

  ofs = ofs0
  Seek #1, ofs

  For I = 1 To 2 'I=1 => 左上の標高,I=2 => 右上の標高
      Get #1, , a1
      ofs = ofs + 1
      Get #1, , a2
      H(I) = a1 * 256 + a2
      ofs = ofs + 1
  Next I

  ofs = ofs0 + 2402
  Seek #1, ofs

  For I = 1 To 2 'I=1 => 左下の標高,I=2 => 右下の標高
      Get #1, , a1
      ofs = ofs + 1
      Get #1, , a2
      H(I + 2) = a1 * 256 + a2
      ofs = ofs + 1
  Next I

  tmp = (1 - (NV - NI)) * 1200
  Ld = tmp - Fix(tmp) '下との距離
  Lc = 1 - Ld '上との距離
  tmp = (EV - EI) * 1200
  La = tmp - Fix(tmp) '左との距離
  Lb = 1 - La '右との距離

  ele1 = La / (La + Lb) * (H(2) - H(1)) + H(1)
  ele2 = La / (La + Lb) * (H(4) - H(3)) + H(3)
  elevation = Round(Lc / (Lc + Ld) * (ele2 - ele1) + ele1, 0)

  Close #1
  Exit Function

errHandler:
  elevation = 0
  Resume Next

End Function

（２）SRTMの生データ出力版


'******************************************************************
'*　　SRTMの生データを検索した結果が出力されます。
'*　　（欠測値の補正が必要か？）
'******************************************************************

Function elevation2(lat As Variant, lon As Variant) As Variant

  Dim I           As Integer

  Dim a1          As Byte 'SRTM上位 1 BYTE
  Dim a2          As Byte 'SRTM下位 1 BYTE

  Dim NV          As Variant
  Dim EV          As Variant
  Dim NI          As Integer
  Dim EI          As Integer

  Dim wDir        As String
  Dim fn          As String
  Dim ofs0        As Long
  Dim ofs         As Long


  On Error GoTo errHandler

  NI = Int(lat) 'SRTM 該当データ存在チェック用（緯度）
  EI = Int(lon) 'SRTM 該当データ存在チェック用（経度）
  wDir = ThisWorkbook.Path 'このBOOKの存在するディレクトリ
  fn = wDir & "\srtm\N" & NI & "E" & EI & ".hgt" 'そこに"SRTM"ディレクトリを作り、*hgt ファイルを入れる

  If Dir(fn) = "" Then
      MsgBox "file open error : 次のファイルが存在しません" & vbCr & fn
      Exit Function
  End If

  Open fn For Binary As #1
  NV = lat
  EV = lon

  '入力データ位置を計算
  ofs = (Fix((1 - (NV - NI)) * 1200) * 1201 + Fix((EV - EI) * 1200)) * 2 + 1

  '入力データにポジショニング
  Seek #1, ofs

  'データ入力（２BYTE）
  Get #1, , a1
  Get #1, , a2
  elevation2 = a1 * 256 + a2 '標高算出

  Close #1
  Exit Function

'欠測データ(-32768)に対するエラー処理
errHandler:
  elevation2 = -9999
  Resume Next

End Function

【使用例】Excelに富士山付近の標高を１分幅で表示する



'******************************************************************
'* 標高計算使用事例
'*　　 富士山のある「N35E138」の標高を緯経度１分間隔で表示する
'******************************************************************

Sub Fuji_mesh()

  Dim dd1         As Integer
  Dim mm1         As Integer
  Dim ss1         As Integer
  Dim dc1         As Integer '開始緯度　北西端の緯度(度)。ファイル名称N35+1度
  Dim mc1         As Integer '北西端の緯度(分)
  Dim sc1         As Variant '北西端の緯度(秒)
  Dim dx1         As Variant
  Dim sx1         As Variant

  Dim dd2         As Integer
  Dim mm2         As Integer
  Dim ss2         As Integer
  Dim dc2         As Integer '開始経度　北西端の経度(度)。ファイル名称E138度
  Dim mc2         As Integer '　　　　　北西端の経度(分)
  Dim sc2         As Variant '　　　　　北西端の経度(秒)
  Dim dx2         As Variant
  Dim sx2         As Variant

  Dim I, J, K     As Integer

  Sheet2.Activate
  Range("A1:B65535").NumberFormatLocal = "0.0000000"

  dc1 = 36: mc1 = 0: sc1 = 0
  dc2 = 138: mc2 = 0: ss2 = 0
  sc1 = dc1 * CLng(3600) + mc1 * 60 + sc1
  sc2 = dc2 * CLng(3600) + mc2 * 60 + sc2

  For I = 0 To 60
      ss1 = I * 60
      sx1 = sc1 - ss1
      dd1 = sx1 \ 3600
      mm1 = (sx1 - dd1 * CLng(3600)) \ 60
      ss1 = sx1 - dd1 * CLng(3600) - mm1 * 60
      dx1 = dd1 + mm1 / 60 + ss1 / 3600

      For J = 0 To 59
          ss2 = J * 60
          sx2 = sc2 + ss2
          dd2 = sx2 \ 3600
          mm2 = (sx2 - dd2 * CLng(3600)) \ 60
          ss2 = sx2 - dd2 * CLng(3600) - mm2 * 60
          dx2 = dd2 + mm2 / 60 + ss2 / 3600

          Cells(I * 60 + J + 1, 1) = dx1
          Cells(I * 60 + J + 1, 2) = dx2
          Cells(I * 60 + J + 1, 3) = elevation(dx1, dx2)
          Cells(I * 60 + J + 1, 4) = elevation2(dx1, dx2)
      Next J
  Next I
End Sub

十目十指

2008年4月7日月曜日

電子国土ポータルに標高を載せる：標高検索

0 件のコメント:

ただいま作業中

ラベル