等值線圖DXF檔的產生

Table of contents

• 背景
• 程式說明
  • 匯入必要的庫
  • 讀取數據文件
  • 篩選DTM
  • 生成 DXF 文件
  • 程式碼
• 程式下載
• 圖檔結果比較

---

背景

• 等高線dxf圖檔事實上是matplotlib.pyplot.contour過程的中間產物，由每條折線上的逐點座標所組成，一般程式設計者是不會需要知道每條等高線上的詳細座標的。
• 除了最終成果的控制目的之外，視覺化驗證也是本項工作一個很重要的理由。
• 處理好的龐大矩陣，也需要藉由本程式來測試看看讀取、篩選過程的工作效率到底怎樣、是否還需要進一步優化。
• 這支程式只有函式版本，可以在API伺服器者、或其他程式場合呼叫

程式說明

這段程式碼讀取存儲在 .dat 文件中的經緯度和數據，並根據給定的地理界限生成等高線圖，最終保存為 PNG 文件。它使用了 matplotlib 庫來繪製圖形，並使用 tempfile 來生成臨時文件名。以下是代碼的詳細解釋：

匯入必要的庫

• numpy 用於數據處理。
• matplotlib 用於繪圖。
• tempfile 和 BytesIO 用於處理臨時文件和內存中的文件。
• ezdxf程式庫，用於創建和操作 DXF 文件。

讀取數據文件

• rd_mem 函數讀取經緯度和DTM數據的 .dat 文件，返回一個包含這3項數據的序列。
• 使用記憶體映射的方式讀取
  • 定義矩陣形狀與讀取檔案同一個指令完成
  • 副程式無法執行exec()指令，GPT建議使用dict或list方式依序讀取、回應呼叫的程式。

篩選DTM

• 雖然等高線圖不是本次專案的目標，但是因為DXF檔案檢視不易，DTM數據量又很大，需要有效率的工具來進行過程驗證。
• dxf 函數根據圖框給定的地理邊界 swLL（西南角）和 neLL（東北角），在數據中以np.where()篩選符合條件的經緯度索引(idx)。
• 如果沒有找到符合條件的數據，返回錯誤信息。

生成 DXF 文件

• dxf 函數根據給定的地理界限 swLL（西南角）和 neLL（東北角），在數據中查找符合條件的經緯度索引。
• 如果沒有找到符合條件的數據，返回錯誤信息。
• 根據找到的經緯度索引計算繪圖的範圍和等高線圖的級別。
• 使用 matplotlib 繪製等高線圖，並將這些線條和文本標籤添加到 DXF 文件中。
• 使用 ezdxf 創建一個新的 DXF 文檔，並保存等高線和標籤。
• 將 DXF 文件保存到內存中，並返回文件名和內存中的文件對象。

程式碼

以下是完整的代碼：

# 匯入必要的庫
import numpy as np
import matplotlib
matplotlib.use('Agg')  # 使用 Agg 后端
import matplotlib.pyplot as plt
import tempfile as tf
from io import BytesIO

# ezdxf程式庫
import ezdxf
from ezdxf import colors
from ezdxf.enums import TextEntityAlignment

# 讀取數據文件
def rd_mem(shape):
    fnames = ['lat', 'lon', 'data']
    d = []
    for f in fnames:
        filename = f + '.dat'
        d.append(np.memmap(filename, dtype='float32', mode='r', shape=shape))
    return d

# 生成(swLL, neLL)範圍的等高線圖
def dxf(swLL, neLL):

    # 讀取GeoTiff檔案的基本參數、用在產生網格點的TWD97座標
    with open('params.txt', 'r') as f:
        line = [i.strip('\n') for i in f][0]
    x0, y0, nx, ny, dx, dy = (float(i) for i in line.split())
    nx, ny = int(nx), int(ny)
    shape = (ny, nx)

    # 讀取全島緯度、精度、以及高程等3項矩陣
    lat, lon, data = rd_mem(shape)

    # 負值歸0，避免拉大最低、最高值的區間
    data = np.where(data < 0, 0, data)

    # 利用經緯度進行範圍切割：共有4個條件都必須符合。
    idx = np.where((lat >= swLL[0]) & (lat <= neLL[0]) & (lon >= swLL[1]) & (lon <= neLL[1]))

    # 'LL not right!'是關鍵字，app.py中會檢核。
    if len(idx[0]) == 0:
        return 'LL not right!', list(swLL) + list(neLL)
    
    # TWD97座標值：
    ## 1維
    x = [x0 + dx * i for i in range(nx)]
    y = [y0 + dy * i for i in range(ny)]
    ## 2維（y,x）
    xg, yg = np.meshgrid(x, y)

    # 取4個方向的TWD97座標極值，用於標定位置
    bounds = [np.min(xg[idx[0], idx[1]]), np.max(xg[idx[0], idx[1]]), np.min(yg[idx[0], idx[1]]), np.max(yg[idx[0], idx[1]])]
    
    # 起訖點的索引
    ib = [x.index(bounds[0]), x.index(bounds[1]), y.index(bounds[2]), y.index(bounds[3])]

    # 範圍內地形的極值、用於計算間距(固定間距5M或1M)
    cmin = data[ib[2]:ib[3]+1, ib[0]:ib[1]+1].min()
    cmax = data[ib[2]:ib[3]+1, ib[0]:ib[1]+1].max()

    # 開始新的DXF檔案
    doc = ezdxf.new(dxfversion="R2010")
    clrs = {5: colors.RED, 1: colors.BLUE}
    
    for intv in [5, 1]:
        N = int((cmax - cmin) // intv)
        if N == 0:
            continue
        levels = np.linspace(cmin, cmax, N)

        # 新的圖面(使用Agg避免顯示器干擾)
        fig, ax = plt.subplots()

        ## 繪圖
        cs = ax.contour(x[ib[0]:ib[1]+1], y[ib[2]:ib[3]+1], data[ib[2]:ib[3]+1, ib[0]:ib[1]+1], levels=levels)
        N = len(cs.collections[:])
        if N == 0:
            continue

        # 新的圖層
        doc.layers.add("TEXTLAYER" + str(intv), color=clrs[intv])
        msp = doc.modelspace()
        for l in range(N):
            lines = cs.collections[l].get_paths()
            for line in lines:
                verts = line.vertices
                xv = verts[:, 0]
                yv = verts[:, 1]
                if len(xv) <= 1:
                    continue
                points = [p for p in zip(xv, yv)]

                # 繪製等高(折)線
                for i in range(len(points) - 1):
                    msp.add_line(points[i], points[i + 1], dxfattribs={"color": clrs[intv]})

            # 每條等高線加上數字標籤
            msp.add_text(
                str(levels[l]),
                dxfattribs={
                    "layer": "TEXTLAYER" + str(intv)
                }).set_placement(points[-1], align=TextEntityAlignment.CENTER)

    ## 存檔、輸出到BytesIO供前後端API傳輸
    ran = tf.NamedTemporaryFile().name.replace('/', '').replace('tmp', '')
    fname = 'terr_' + ran + '.dxf'
    output = BytesIO()
    doc.write(output, fmt='bin')
    output.seek(0)  # 重置指針位置
    doc.saveas('./dxfs/' + fname)
    
    return fname, output

這段程式碼在給定的地理界限內從 .dat 文件中讀取數據，並生成一個包含等高線圖的 DXF 文件。具體步驟如下：

這段代碼的關鍵步驟包括：

• 讀取地理數據。
• 查找給定界限內的數據。
• 繪製等高線圖
• 生成DXF 文件，其中包含了數據的等高線和對應的標籤，便於在 CAD 軟件中進行查看和分析。

程式下載

• 副程式版本<div markdown="span" class="alert alert-success" role="alert"> Download: mem2dxf.py</div>

，供API呼叫、自動回復。

圖檔結果比較

matplotlib等高線圖	dxf檔案