python非均勻網格等高線 [Kurt]

你可以使用Python中的NumPy、SciPy和Matplotlib库来实现这个需求。具体步骤如下：

1. 計算梯度（高低變化速度）: 使用numpy.gradient計算出每個點在x和y方向上的梯度，並結合這兩個梯度計算總的高低變化速度。
2. 判斷變化速度是否超過閾值: 設定一個閾值，將變化速度超過該閾值的點標記出來。
3. 繪製等高線圖和標記位置: 使用matplotlib來繪製等高線圖，並將標記出的點疊加在圖上。

以下是實現這一需求的Python代碼示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# 假設已有非均勻網格上的數據
# x, y 為非均勻的網格點坐標
# z 為對應的高度值

# 生成範例數據 (假設已知x, y, z)
x = np.random.rand(100) * 10
y = np.random.rand(100) * 10
z = np.sin(x) + np.cos(y)

# 將非均勻數據插值到均勻網格
xi = np.linspace(min(x), max(x), 200)
yi = np.linspace(min(y), max(y), 200)
xi, yi = np.meshgrid(xi, yi)
zi = griddata((x, y), z, (xi, yi), method='cubic')

# 計算梯度（高低變化速度）
dzdx, dzdy = np.gradient(zi, xi[0,:], yi[:,0])

# 計算總的變化速度
speed = np.sqrt(dzdx**2 + dzdy**2)

# 設置閾值
threshold = 1.0  # 根據實際情況調整閾值

# 標記超過閾值的點
markers = np.where(speed > threshold)

# 繪製等高線圖
plt.figure()
cp = plt.contourf(xi, yi, zi, cmap='viridis')
plt.colorbar(cp)

# 在圖上標記變化速度超過閾值的位置
plt.scatter(xi[markers], yi[markers], color='red', marker='x')

plt.title('Contour plot with marked rapid elevation changes')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

1. 插值到均勻網格: 由於你有非均勻的x、y資料點，為了方便後續計算，使用scipy.interpolate.griddata將數據插值到均勻網格上。
2. 計算梯度: 使用numpy.gradient來計算z值在x和y方向上的梯度，即高低變化的速度。
3. 標記超過閾值的點: 根據指定的閾值，使用numpy.where找出變化速度超過該閾值的點。
4. 繪圖與標記: 使用matplotlib來繪製等高線圖，並在圖中標出變化速度超過閾值的位置。

這樣，你就能以最少的運算標記出高低變化速度超過某閾值的位置，並將其顯示在等高線圖上。

=============================================

要在三维地形数据中标记出高度变化速度超过某个阈值的位置，可以遵循以下步骤：

1. 计算邻近点之间的高度变化速度：通过计算每个点与其相邻点之间的高度差异，然后根据距离计算出速度。
2. 标记超过阈值的位置：对每个点的高度变化速度与给定的阈值进行比较，标记出超过阈值的位置。

以下是一个Python代码示例，假设x, y, z为一维数组表示地形的坐标和高度：

import numpy as np

def compute_slope(x, y, z):
    # 计算x, y方向上的距离
    dx = np.diff(x)
    dy = np.diff(y)
    # 计算斜率 (高度差除以距离)
    dz = np.diff(z)
    dist = np.sqrt(dx**2 + dy**2)
    slope = dz / dist
    return slope

def mark_high_slope_regions(x, y, z, threshold):
    # 计算所有点之间的斜率
    slope = compute_slope(x, y, z)
    
    # 标记斜率超过阈值的点
    high_slope_indices = np.where(np.abs(slope) > threshold)[0]
    
    # 返回标记的点 (考虑到diff减少了一个点的索引)
    return high_slope_indices

# 示例数据
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 1, 1, 0])
z = np.array([0, 1, 3, 2, 1])

# 设置斜率阈值
threshold = 1.0 # 45度, 60度 tan(60deg) = 1.732

# 标记高斜率区域
high_slope_indices = mark_high_slope_regions(x, y, z, threshold)
print("High slope indices:", high_slope_indices)

解释：

1. 计算斜率：
2. • 使用np.diff()计算每个点和其下一个点之间的差异（在x, y, z三个维度上）。
   • 使用欧几里得距离公式计算两个点之间的距离。
   • 斜率即为高度变化除以距离。
2. 标记高斜率区域：
3. • 使用np.where()函数找出斜率超过阈值的位置。

注意事项：

• 这个方法计算的是相邻点之间的斜率。如果你的数据不是均匀网格，你可能需要考虑更多邻居（如2D或3D网格上的多个方向）。
• 你可以根据需要调整threshold以标记不同的变化速度。

这个方法只需要一次遍历数据，计算复杂度为O(n)，非常高效。

============================================

找出水平及垂直方向（誤差1度）高低變化速度超過某閾值的位置，你可以採取以下步驟。這裡的「水平及垂直方向」應指與X軸、Y軸平行或近似平行的方向，並且允許一度的誤差。

步驟：

1. 計算梯度（高低變化速度）：使用numpy.gradient計算每個點在X方向和Y方向上的梯度。
2. 計算角度：根據梯度來計算梯度方向與X軸和Y軸的夾角，確保這些方向與水平或垂直方向之間的誤差不超過1度。
3. 過濾符合條件的點：過濾出梯度超過某閾值並且方向接近水平或垂直的點。
4. 標記並繪圖：使用matplotlib繪製等高線圖，並將這些點標記在圖上。

以下是具體的Python代碼實現：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# 假設已有非均勻網格上的數據
x = np.random.rand(100) * 10
y = np.random.rand(100) * 10
z = np.sin(x) + np.cos(y)

# 將非均勻數據插值到均勻網格
xi = np.linspace(min(x), max(x), 200)
yi = np.linspace(min(y), max(y), 200)
xi, yi = np.meshgrid(xi, yi)
zi = griddata((x, y), z, (xi, yi), method='cubic')

# 計算梯度（高低變化速度）
dzdx, dzdy = np.gradient(zi, xi[0,:], yi[:,0])

# 計算變化速度
speed = np.sqrt(dzdx**2 + dzdy**2)

# 計算梯度方向的角度（與X軸的夾角）
angles = np.arctan2(dzdy, dzdx) * 180 / np.pi

# 設定角度和閾值
angle_threshold = 1  # 誤差角度
speed_threshold = 1.0  # 變化速度閾值

# 過濾出符合條件的點（方向接近水平或垂直，並且速度超過閾值）
markers = np.where(
    (speed > speed_threshold) &
    ((np.abs(angles) <= angle_threshold) | (np.abs(angles - 90) <= angle_threshold) | 
     (np.abs(angles + 90) <= angle_threshold) | (np.abs(angles - 180) <= angle_threshold))
)

# 繪製等高線圖
plt.figure()
cp = plt.contourf(xi, yi, zi, cmap='viridis')
plt.colorbar(cp)

# 在圖上標記變化速度超過閾值的位置
plt.scatter(xi[markers], yi[markers], color='red', marker='x')

plt.title('Contour plot with marked rapid elevation changes (near horizontal or vertical)')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

代碼說明：

1. 計算梯度：dzdx和dzdy分別是z在X方向和Y方向的梯度。
2. 計算角度：使用np.arctan2計算梯度方向相對於X軸的角度，並轉換為度數。
3. 過濾符合條件的點：判斷角度是否在允許的誤差範圍內，並且高低變化速度是否超過閾值。
4. 繪製圖像：使用matplotlib繪製等高線圖，並將符合條件的點標記在圖上。

這個方法通過計算梯度和角度，能夠高效地找出接近水平或垂直方向的高低變化點，並將它們可視化。

=============================================

要將一組非均勻網格的地形資料（x, y, z）在excel內用VBA轉換成網格格式並輸出到新工作表，您可以使用VBA來自動化這個過程。以下是步驟和範例VBA程式碼。

步驟：

1. 準備數據：確保您的原始數據存放在Excel的工作表中。例如，x座標在A列，y座標在B列，z高度值在C列。
2. 編寫VBA腳本：此腳本將會自動讀取x, y, z資料，然後將它們轉換成均勻的網格格式。
3. 創建新工作表：將轉換後的網格數據輸出到一個新的工作表。

VBA 範例程式碼：

Sub ConvertToGrid()
    Dim ws As Worksheet
    Dim newWs As Worksheet
    Dim xRange As Range, yRange As Range, zRange As Range
    Dim xValues As Object, yValues As Object, zDict As Object
    Dim x As Double, y As Double, z As Double
    Dim row As Long, col As Long
    Dim xMin As Double, xMax As Double, yMin As Double, yMax As Double
    Dim xStep As Double, yStep As Double
    
    ' 設定來源工作表
    Set ws = ThisWorkbook.Sheets("Sheet1") ' 替換為您的工作表名稱
    
    ' 獲取範圍
    Set xRange = ws.Range("A2:A" & ws.Cells(ws.Rows.Count, "A").End(xlUp).Row)
    Set yRange = ws.Range("B2:B" & ws.Cells(ws.Rows.Count, "B").End(xlUp).Row)
    Set zRange = ws.Range("C2:C" & ws.Cells(ws.Rows.Count, "C").End(xlUp).Row)
    
    ' 創建新工作表
    Set newWs = ThisWorkbook.Sheets.Add
    newWs.Name = "GridData"
    
    ' 收集x, y, z資料到字典
    Set xValues = CreateObject("Scripting.Dictionary")
    Set yValues = CreateObject("Scripting.Dictionary")
    Set zDict = CreateObject("Scripting.Dictionary")
    
    For row = 1 To xRange.Rows.Count
        x = xRange.Cells(row, 1).Value
        y = yRange.Cells(row, 1).Value
        z = zRange.Cells(row, 1).Value
        
        If Not xValues.exists(x) Then xValues.Add x, x
        If Not yValues.exists(y) Then yValues.Add y, y
        
        zDict(x & "," & y) = z
    Next row
    
    ' 轉換為陣列
    xMin = Application.Min(xValues.keys)
    xMax = Application.Max(xValues.keys)
    yMin = Application.Min(yValues.keys)
    yMax = Application.Max(yValues.keys)
    
    ' 設定網格步長
    xStep = (xMax - xMin) / (xValues.Count - 1)
    yStep = (yMax - yMin) / (yValues.Count - 1)
    
    ' 填充新工作表的標題
    newWs.Cells(1, 1).Value = "X\Y"
    For col = 1 To yValues.Count
        newWs.Cells(1, col + 1).Value = yMin + (col - 1) * yStep
    Next col
    
    For row = 1 To xValues.Count
        newWs.Cells(row + 1, 1).Value = xMin + (row - 1) * xStep
        For col = 1 To yValues.Count
            x = xMin + (row - 1) * xStep
            y = yMin + (col - 1) * yStep
            If zDict.exists(x & "," & y) Then
                newWs.Cells(row + 1, col + 1).Value = zDict(x & "," & y)
            Else
                newWs.Cells(row + 1, col + 1).Value = "" ' 未知值
            End If
        Next col
    Next row

End Sub

說明：

• xRange、yRange、zRange：定義x, y, z資料的範圍。
• xValues、yValues、zDict：使用字典來存儲x和y的唯一值，以及x,y對應的z值。
• xStep、yStep：計算網格中的步長。
• newWs：將結果輸出到一個新的工作表。

此程式碼會將非均勻網格的資料轉換成均勻的網格，並在新工作表中列出所有的x和y值對應的z值。未找到的z值會顯示為空白。您可以根據需求調整這個程式碼。