todashuta · August 20, 2024 16:48 · Aug 20, 2024 · Jul 28, 2020 · Jul 28, 2020
diff --git a/lem2gtif.py b/lem2gtif.py
@@ -2,7 +2,7 @@
 import re
 import glob
 import time
-import gdal, osr
+from osgeo import (gdal, osr)
 import numpy as np
 
 

diff --git a/README.md b/README.md
@@ -5,10 +5,11 @@ lem ファイルを GeoTIFF に変換するスクリプトです。
 
 * lem ファイル仕様
   * https://www.gsi.go.jp/MAP/CD-ROM/dem5m/doc/info5m1.htm
+  * ただし、データ間隔によりレコード数等は説明と合致しないことがあります。
 
 ## 使い方
 
-現時点では） lem2gtif.py が存在するフォルダ内の lem ファイルを GeoTIFF に変換するようになっているので、このファイルを対象のフォルダにコピーし、実行してください。
+（現時点では） lem2gtif.py が存在するフォルダ内の lem ファイルを GeoTIFF に変換するようになっているので、このファイルを対象のフォルダにコピーし、実行してください。
 
 ### 使い方（ QGIS をインストールした Windows ）
 
@@ -29,5 +30,5 @@ C:\>cd C:\Users\username\Downloads\05mdtm02epsg2454
 最後に `python3 lem2gtif.py` とスクリプトを実行します。
 
 ```
-C:\>python3 lem2gtif.py
+C:\Users\username\Downloads\05mdtm02epsg2454>python3 lem2gtif.py
 ```
diff --git a/README.md b/README.md
@@ -0,0 +1,33 @@
+# lem2gtif.py
+
+## 概要
+lem ファイルを GeoTIFF に変換するスクリプトです。
+
+* lem ファイル仕様
+  * https://www.gsi.go.jp/MAP/CD-ROM/dem5m/doc/info5m1.htm
+
+## 使い方
+
+現時点では） lem2gtif.py が存在するフォルダ内の lem ファイルを GeoTIFF に変換するようになっているので、このファイルを対象のフォルダにコピーし、実行してください。
+
+### 使い方（ QGIS をインストールした Windows ）
+
+lem2gtif.py を、 lem ファイルが入っているフォルダに入れ、スタートメニューから「OSGeo4W Shell」を起動します。
+
+まず `py3_env` を実行し、 Python3 を使えるように環境設定を行います。
+
+```
+C:\>py3_env
+```
+
+次に、 lem2gtif.py があるフォルダに `cd` コマンドで移動します。エクスプローラのアドレスバーからアドレスをコピーすると楽です。
+
+```
+C:\>cd C:\Users\username\Downloads\05mdtm02epsg2454
+```
+
+最後に `python3 lem2gtif.py` とスクリプトを実行します。
+
+```
+C:\>python3 lem2gtif.py
+```
diff --git a/lem2gtif.py b/lem2gtif.py
@@ -0,0 +1,136 @@
+import os
+import re
+import glob
+import time
+import gdal, osr
+import numpy as np
+
+
+"""
+lem ファイルを GeoTIFF に変換するスクリプトです。
+
+lem ファイルの仕様
+https://www.gsi.go.jp/MAP/CD-ROM/dem5m/doc/info5m1.htm
+
+
+（現時点では） lem2gtif.py が存在するフォルダ内の lem ファイルを
+GeoTIFF に変換するようになっているので、このファイルを対象の
+フォルダにコピーし、実行してください。
+"""
+
+
+
+def main():
+    # lem2gtif.py と同じフォルダ内に存在する lem ファイルを走査
+    basedir = os.path.dirname(os.path.abspath(__file__))
+    for lem_filepath in glob.glob(os.path.join(basedir, '*.lem')):
+        try:
+            # lem ファイルを GeoTIFF に変換
+            lem_to_gtif(lem_filepath)
+        except (ValueError, FileNotFoundError) as err:
+            print(f"Error: {err}")
+            print()
+
+
+
+def lem_to_gtif(lem_filepath):
+    if not re.search(r'\.lem$', lem_filepath):
+        raise ValueError("lem ファイル以外が指定されました。")
+
+    # lem ファイルと同じ場所に csv 形式のヘッダファイルがあるはず
+    csv_filepath = re.sub(r'\.lem$', '.csv', lem_filepath)
+    # lem ファイルと同じ場所に変換した TIFF ファイルを出力する
+    tif_filepath = re.sub(r'\.lem$', '.tif', lem_filepath)
+
+    lem_filename = os.path.basename(lem_filepath)
+    csv_filename = os.path.basename(csv_filepath)
+    tif_filename = os.path.basename(tif_filepath)
+
+    now = time.strftime('%H:%M:%S', time.localtime())
+    print(f"[{now}] convert to {tif_filename}")
+
+    # ファイルが存在しているかチェック
+    if not os.path.exists(lem_filepath):
+        raise FileNotFoundError(f"{lem_filename} がみつかりません。")
+    if not os.path.exists(csv_filepath):
+        raise FileNotFoundError(f"{csv_filename} がみつかりません。")
+
+    # ヘッダファイルを読み取る
+    header = read_header_file(csv_filepath)
+    if header is None:
+        raise ValueError(f"{csv_filename} のパースに失敗しました。")
+
+    # GeoTIFF ファイル作成
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(header['crs'])
+    driver = gdal.GetDriverByName('GTiff')
+    raster = driver.Create(tif_filepath, header['cols'], header['rows'], 1, gdal.GDT_Float32)
+    raster.SetGeoTransform((header['left'], header['xres'], 0, header['top'], 0, -header['yres']))
+    raster.SetProjection(srs.ExportToWkt())
+    band = raster.GetRasterBand(1)
+    band.SetNoDataValue(-9999.0)
+
+    with open(lem_filepath) as lem_f:
+        row = 0
+        # lem ファイルから1行読み取り、 GeoTIFF ファイルに1行書き込む
+        for line in lem_f:
+            line = line.replace('-1111', '-9999')
+            array = np.array([[-9999.0 if line[5*i+10:5*i+15] == '-9999' else int(line[5*i+10:5*i+15])/10 for i in range(header['cols'])]])
+            band.WriteArray(array, yoff=row)
+            row += 1
+        band.FlushCache()
+
+
+
+def read_header_file(csv_filepath):
+    csv_filename = os.path.basename(csv_filepath)
+    header = {}
+
+    # ヘッダファイルを読み込む
+    with open(csv_filepath, 'r', encoding='cp932') as f:
+        for line in f:
+            k, v = line.strip().split(',')
+            if k == '東西方向の点数':
+                header['cols'] = int(v)
+            elif k == '南北方向の点数':
+                header['rows'] = int(v)
+            elif k == '東西方向のデータ間隔':
+                header['xres'] = float(v)
+            elif k == '南北方向のデータ間隔':
+                header['yres'] = float(v)
+            elif k == '平面直角座標系番号' or k == '座標系番号':
+                header['crs'] = int(v) + 6668
+            elif k == '区画左下X座標' or k == '区画左下Ｘ座標':
+                header['bottom'] = int(v) / 100
+            elif k == '区画左下Y座標' or k == '区画左下Ｙ座標':
+                header['left'] = int(v) / 100
+            elif k == '区画右上X座標' or k == '区画右上Ｘ座標':
+                header['top'] = int(v) / 100
+            elif k == '区画右上Y座標' or k == '区画右上Ｙ座標':
+                header['right'] = int(v) / 100
+
+    # 必要な情報が読み取れなかったら例外を投げる
+    if 'cols' not in header:
+        raise ValueError(f"{csv_filename} に「東西方向の点数」がありません。")
+    if 'rows' not in header:
+        raise ValueError(f"{csv_filename} に「南北方向の点数」がありません。")
+    if 'xres' not in header:
+        raise ValueError(f"{csv_filename} に「東西方向のデータ間隔」がありません。")
+    if 'yres' not in header:
+        raise ValueError(f"{csv_filename} に「南北方向のデータ間隔」がありません。")
+    if 'crs' not in header:
+        raise ValueError(f"{csv_filename} に「座標系番号」がありません。")
+    if 'bottom' not in header:
+        raise ValueError(f"{csv_filename} に「区画左下Ｘ座標」がありません。")
+    if 'left' not in header:
+        raise ValueError(f"{csv_filename} に「区画左下Ｙ座標」がありません。")
+    if 'top' not in header:
+        raise ValueError(f"{csv_filename} に「区画右上Ｘ座標」がありません。")
+    if 'right' not in header:
+        raise ValueError(f"{csv_filename} に「区画右上Ｙ座標」がありません。")
+
+    return header
+
+
+
+main()
No results found