FAQ/ ASC形式からインポート可能なNetCDF形式への変換

テクニカルリファレンス 入出力機能/ 地形・粗度・盛土のインポート機能/ NetCDF形式 について、インポート可能なnetCDFファイルを作成する方法を教えてください。

回答

netCDFを作成する方法には、いくつかの方法があります。

• テキストファイルから ncgen.exe というツールを使って netCDF形式に変換する方法
• QGISなどのGISソフトを使って、メッシュデータをASC形式でエクスポートし、PowerShellなどのプログラミング言語を使って netCDF形式を作成する方法

ここでは、2番目の方法について説明します。 1番目の方法については、こちらをご参照ください。

1. GISソフトなどを使って地形・粗度・盛土などのデータを作成し、ESRI ASC Raster形式で作成します。
   サンプルデータを示します。これらのファイルを デスクトップに保存します。

   内容	サンプルファイル
   標高データ	glev_asc.txt
   粗度データ	roughness_asc.txt
   盛土天端データ	leveeHeight_asc.txt
   盛土フラグデータ	leveeFlag_asc.txt

2. PowerShellをインストールします。
   Microsoft Store から、PowerShell をインストールします。
   Microsoft Store を使わずにインストールする場合は Windows への PowerShell のインストール を参照してください。
   なお、ここで指している「PowerShell」とは、Windowsで標準インストールされている「Windows PowerShell」のことではございません。スクリプトの実行にあたってはPowerShell からインストールしたものを使用してください。

3. netCDFライブラリをインストールします。
   NetCDF の公式サイトの Installing and Using netCDF-C Libraries in Windows から netCDF4 64-bit版を取得します。
   2021年10月現在の最新バージョンは、netCDF4.8.0-NC4-64.exe です。
   以降の説明では、このバージョン4.8.0をインストールしたものとして説明します。
   既定では C:\Program Files\netCDF 4.8.0\ に netCDFライブラリがインストールされます。

4. PowerShellのためのnetCDFライブラリをインストールします。
   NetcdfCs.dll をダウンロードし、C:\Program Files\netCDF 4.8.0\bin にコピーします。
   必要に応じて、.NET 7.0 をインストールします。Windows に .NET をインストールする
   なお、NetcdfCs.dllについて、Windowsセキュリティにより実行が拒否される場合はWindowsの設定からセキュリティの許可をしてください。

5. PowerShellスクリプトを作成します。 サンプルを示します。

   # 入力ファイル
   $GlevAscFile = Get-Item "glev_asc.txt"
   $roughnessAscFile = Get-Item "roughness_asc.txt"
   $leveeHeightAscFile = Get-Item "leveeheight_asc.txt"
   $leveeFlagAscFile = Get-Item "leveeflag_asc.txt"
   
   # 出力ファイル
   $ncFileName = Join-Path (Get-Location) "DiovistaDomainData.nc"
   
   # dllファイル
   $dllFileName = Get-Item "C:\Program Files\netCDF 4.8.0\bin\netcdfCs.dll"
   [void]([System.Reflection.Assembly]::LoadFile($dllFileName))
   
   # エラー時停止
   $ErrorActionPreference = "Stop"
   
   enum NcType {
       NC_BYTE = 1
       NC_SHORT = 3
       NC_INT = 4
   }
   class AscData {
       [System.IO.FileInfo]$file
       [int]$ncols
       [int]$nrows
       [double]$xllcorner
       [double]$yllcorner
       [double]$dx
       [double]$dy
       [double]$nodataValue
       [double[]] $data
       [int[]] $intData
       [int16[]] $shortData
       [byte[]] $byteData
       [int]$validMin
       [double]$scaleFactor
   
       # ctor
       AscData([System.IO.FileInfo]$file, [NcType]$ncType, [int]$validMin, [double]$scaleFactor) {
           $this.file = $file
           $this.validMin = $validMin
           $this.scaleFactor = $scaleFactor
           $scaleFactorInverse = 1.0 / $scaleFactor
           $headerLines = 7
           $lines = Get-Content -Encoding utf8 $this.file.FullName
           foreach ($headerLine in $lines[0..($headerLines - 1)]) {
               $header = $headerLine -split " +"
               if ($header[0] -like "NCOLS") {
                   $this.ncols = $header[1]
               }
               elseif ($header[0] -like "NROWS") {
                   $this.nrows = $header[1]
               }
               elseif ($header[0] -like "XLLCORNER") {
                   $this.xllcorner = $header[1]
               }
               elseif ($header[0] -like "YLLCORNER") {
                   $this.yllcorner = $header[1]
               }
               elseif ($header[0] -like "DX") {
                   $this.dx = $header[1]
               }
               elseif ($header[0] -like "DY") {
                   $this.dy = $header[1]
               }
               elseif ($header[0] -like "NODATA_VALUE") {
                   $this.nodataValue = $header[1]
               }
               else {
                   throw [System.Exception] (("{0}: found not supported header: {1}" -f $this.file, $header[0]))
               }
           }
           $this.DumpHeader()
           $this.data = New-Object double[] ($this.ncols * $this.nrows)
           for ($j = 0; $j -lt $this.nrows; $j++) {
               $ws = $lines[$j + $headerLines] -split " +"
               $i0 = $j * $this.ncols
               for ($i = 0; $i -lt $this.ncols; $i++) {
                   $this.data[$i0 + $i] = $ws[$i]
               }
           }
           if ($ncType -eq [NcType]::NC_BYTE) {
               $this.byteData = New-Object byte[] ($this.ncols * $this.nrows)
               for ($k = 0; $k -lt $this.data.Length; $k++) {
                   $value = $this.data[$k]
                   if ($value -eq $this.nodataValue) {
                       $this.byteData[$k] = 0;
                   } else {
                       $this.byteData[$k] = $value
                   }
               }
           }
           elseif($ncType -eq [NcType]::NC_SHORT) {
               $this.shortData = New-Object int16[] ($this.ncols * $this.nrows)
               for ($k = 0; $k -lt $this.data.Length; $k++) {
                   $value = $this.data[$k]
                   if ($value -eq $this.nodataValue) {
                       $this.shortData[$k] = $validMin - 1;
                   } else {
                       $this.shortData[$k] = [int16]($value * $scaleFactorInverse)
                   }
               }
           }
           elseif($ncType -eq [NcType]::NC_INT) {
               $this.intData = New-Object int[] ($this.ncols * $this.nrows)
               for ($k = 0; $k -lt $this.data.Length; $k++) {
                   $value = $this.data[$k]
                   if ($value -eq $this.nodataValue) {
                       $this.intData[$k] = $validMin - 1;
                   } else {
                       $this.intData[$k] = [int]($value * $scaleFactorInverse)
                   }
               }
           }
   
       }
       [void] DumpHeader() {
           $message = "file: {0}`n" -f $this.file.FullName
           $message += "ncols: {0}`n" -f $this.ncols
           $message += "nrows: {0}`n" -f $this.nrows
           $message += "xllcorner: {0}`n" -f $this.xllcorner
           $message += "yllcorner: {0}`n" -f $this.yllcorner
           $message += "dx: {0}`n" -f $this.dx
           $message += "dy: {0}`n" -f $this.dy
           $message += "nodataValue: {0}`n" -f $this.nodataValue
           Write-Host $message
       }
   }
   
   # load asc files
   $GlevAsc = New-Object AscData($GlevAscFile, [NcType]::NC_INT, -100000, 0.01)
   $roughnessAsc = New-Object AscData($roughnessAscFile, [NcType]::NC_SHORT, 0, 0.001)
   $leveeHeightAsc = New-Object AscData($leveeHeightAscFile, [NcType]::NC_INT, -100000, 0.01)
   $leveeFlagAsc = New-Object AscData($leveeFlagAscFile, [NcType]::NC_BYTE, 0, 1.0)
   $asc = $GlevAsc
   
   # create netcdf file
   $ncFile = New-Object netcdfCs.NcFile
   if ((Test-Path $ncFileName) -eq $true) {
       Remove-Item $ncFileName
   }
   # Options(bool isNetCDF4, bool forceOverwrite, int compressionLevel
   $ncOption = New-Object NetcdfCs.NcFile+Options($true, $true, 6)
   $ncFile.Create($ncFileName, $ncOption)
   Write-Host ("Create {0}" -f $ncFileName)
   
   # global attributes
   # 表 9(1)ファイル識別子
   [void]$ncFile.AddAttribute("title", [netcdfCs.NcType]::NC_CHAR, "DioVISTA calculation domain data")
   # 使用するnetcdf規約のバージョン（CF-1.6）
   [void]$ncFile.AddAttribute("Conventions", [netcdfCs.NcType]::NC_CHAR, "CF-1.6")
   # その他コメント
   [void]$ncFile.AddAttribute("comment", [netcdfCs.NcType]::NC_CHAR, ("generated on {0} by netcdfPs1_Convert-From-Asc-To-DiovistaDomainNc.ps1" -f (Get-Date)))
   
   # dimensions
   $dimLon = $ncFile.AddDimension("lon", $asc.ncols)
   $dimLat = $ncFile.AddDimension("lat", $asc.nrows)
   
   # variables
   $varLon = $ncFile.AddVariable("lon", [netcdfCs.NcType]::NC_SHORT, $dimLon.DimId)
   [void]$varLon.AddAttribute("add_offset", [netcdfCs.NcType]::NC_DOUBLE, $asc.xllcorner + $asc.dx * 0.5)
   [void]$varLon.AddAttribute("scale_factor", [netcdfCs.NcType]::NC_DOUBLE, $asc.dx)
   [void]$varLon.AddAttribute("units", [netcdfCs.NcType]::NC_CHAR, "degrees_east")
   
   $varLat = $ncFile.AddVariable("lat", [netcdfCs.NcType]::NC_SHORT, $dimLat.DimId)
   [void]$varLat.AddAttribute("add_offset", [netcdfCs.NcType]::NC_DOUBLE, $asc.yllcorner + $asc.dy * 0.5)
   [void]$varLat.AddAttribute("scale_factor", [netcdfCs.NcType]::NC_DOUBLE, $asc.dy)
   [void]$varLat.AddAttribute("units", [netcdfCs.NcType]::NC_CHAR, "degrees_north")
   
   $varGlev = $ncFile.AddVariable("glev", [netcdfCs.NcType]::NC_INT, ($dimLat.DimId, $dimLon.DimId))
   [void]$varGlev.AddAttribute("valid_min", [netcdfCs.NcType]::NC_INT, [int]$GlevAsc.validMin)
   [void]$varGlev.AddAttribute("scale_factor", [netcdfCs.NcType]::NC_DOUBLE, $GlevAsc.scaleFactor)
   [void]$varGlev.AddAttribute("units", [netcdfCs.NcType]::NC_CHAR, "m")
   
   $varRoughness = $ncFile.AddVariable("roughness", [netcdfCs.NcType]::NC_SHORT, ($dimLat.DimId, $dimLon.DimId))
   [void]$varRoughness.AddAttribute("valid_min", [netcdfCs.NcType]::NC_SHORT, [int16]$roughnessAsc.validMin)
   [void]$varRoughness.AddAttribute("scale_factor", [netcdfCs.NcType]::NC_DOUBLE, $roughnessAsc.scaleFactor)
   [void]$varRoughness.AddAttribute("units", [netcdfCs.NcType]::NC_CHAR, "m^(-1/3).s")
   
   $varLeveeHeight = $ncFile.AddVariable("leveeHeight", [netcdfCs.NcType]::NC_INT, ($dimLat.DimId, $dimLon.DimId))
   [void]$varLeveeHeight.AddAttribute("valid_min", [netcdfCs.NcType]::NC_INT, [int]$leveeHeightAsc.validMin)
   [void]$varLeveeHeight.AddAttribute("scale_factor", [netcdfCs.NcType]::NC_DOUBLE, $leveeHeightAsc.scaleFactor)
   [void]$varLeveeHeight.AddAttribute("units", [netcdfCs.NcType]::NC_CHAR, "m")
   
   $varLeveeFlag = $ncFile.AddVariable("leveeFlag", [netcdfCs.NcType]::NC_BYTE, ($dimLat.DimId, $dimLon.DimId))
   [void]$varLeveeFlag.AddAttribute("flag_meanings", [netcdfCs.NcType]::NC_CHAR, "east north")
   $flagMasks = New-Object byte[] (2)
   $flagMasks[0] = 1
   $flagMasks[1] = 2
   [void]$varLeveeFlag.AddAttribute("flag_masks", [netcdfCs.NcType]::NC_BYTE, $flagMasks)
   
   $ncFile.SaveDefinitions()
   # Write-Host "Save Definitions"
   
   # set variable data
   # x
   $shortData = New-Object int16[] $dimLon.Length
   for ($i = 0; $i -lt $asc.ncols; $i++) {
       $shortData[$i] = $i
   }
   $varLon.Values = $shortData
   $varLon.SaveData($ncFile.NcId)
   
   # y
   $shortData = New-Object int16[] $dimLat.Length
   for ($j = 0; $j -lt $asc.nrows; $j++) {
       $shortData[$j] = $asc.nrows - 1 - $j
   }
   $varLat.Values = $shortData
   $varLat.SaveData($ncFile.NcId)
   
   # Glev
   $varGlev.Values = $GlevAsc.intData
   $varGlev.SaveData($ncFile.NcId)
   
   # roughness
   $varRoughness.Values = $roughnessAsc.shortData
   $varRoughness.SaveData($ncFile.NcId)
   
   # leveeheight
   $varLeveeHeight.Values = $leveeheightAsc.intData
   $varLeveeHeight.SaveData($ncFile.NcId)
   
   # leveeheight
   $varLeveeFlag.Values = $leveeflagAsc.byteData
   $varLeveeFlag.SaveData($ncFile.NcId)
   
   $ncFile.Close()
   # Write-Host ("Close {0}" -f $ncFileName)
   

   • 2~5行目に、入力するファイルを指定します。
   • 8行目に、出力するファイルをフルパスで指定します。ここでは、デスクトップに DiovistaDomainData.nc を出力する設定をしています。
   • 11行目に、ステップ2でインストールした、PowerShell用のnetCDFライブラリのパスを指定します。ここでは、"C:\Program Files\netCDF 4.8.0\bin\netcdfCs.dll" を指定しています。
   • スクリプトをテキストエディタに貼り付け、デスクトップにファイル名 netcdfPs1_Convert-From-Asc-To-DiovistaDomainNc.ps1 として保存します。

6. PowerShell でスクリプトを実行します。
   + R を押し、出現したダイアログ [ファイル名を指定して実行] に pwsh と入力します。
   PowerShell が立ち上がるので、次の文字列をコピーして貼り付けて実行します。

   cd $env:USERPROFILE\Desktop
   Set-ExecutionPolicy RemoteSigned -Scope Process
   .\netcdfPs1_Convert-From-Asc-To-DiovistaDomainNc.ps1
   

   デスクトップに DiovistaDomainData.nc が作成されます。

7. 作成した netCDF をDioVISTA に取り込みます。

   1. DioVISTAを立ち上げます。

   2. メニュー > [ファイル] > [テンプレートから新規作成] > [緯度経度座標] を選択します。
      

   3. プロジェクト > [計算領域] > [NetCDFからインポート] を選択します。
      

   4. 出現したダイアログのファイル名に、 DiovistaDomainData.nc を指定します。
      メッシュサイズに 25 m を指定します。
      

   5. DioVISTAにデータが取り込まれました。
      

---

最終更新日: 2023-09-19