利用高程数据计算坡度、坡向及山体阴影

MeteoInfo

2024年05月28日 14:14

收录于文集

共51篇

高程数据体现了地形高低变化，可以通过读取高程数据二维数组，用imshow将高程用不同颜色表示出来绘图显示地形。

 代码块
Python
自动换行
复制代码
fn = &#39;D:/Temp/nc/etopo2_new.nc&#39;
f = addfile(fn)
data = f[&#39;z&#39;][&#39;25:50&#39;,&#39;72:115&#39;]

#Plot
axesm()
geoshow(&#39;country&#39;, edgecolor=&#39;k&#39;)
imshow(data, 40, cmap=&#39;MPL_terrain&#39;)
colorbar()复制成功

上图中地形高低变化很清楚，但还缺乏立体感。在二维图中体现地形立体感需要用到山体阴影，要先计算出每个格点的坡度、坡向。

首先读取高程数据，etopo2_new.nc文件中包含了全球0.33度左右分辨率的高程数据，下面的语句通过设定经纬度范围读取部分数据，并将数据的经纬度也读取出来。

 代码块
Python
自动换行
复制代码
fn = &#39;D:/Temp/nc/etopo2_new.nc&#39;
f = addfile(fn)
data = f[&#39;z&#39;][&#39;25:50&#39;,&#39;72:115&#39;]
data = data.astype(&#39;float&#39;)
lon = data.dimvalue(1)
lat = data.dimvalue(0)复制成功

坡度计算需要计算格点高程变化和水平距离的比值，再取反正切。x, y方向高程的变化可以用梯度函数gradient计算，格点水平距离0.33度大约为3000米（单位和高程单位保持一致），坡度用下面的代码计算，zf 是高程的放大系数，更高的zf值会带来更大的立体感。

 代码块
Python
自动换行
复制代码
x, y = np.gradient(data, 3000, 3000)
zf = 1.    #slope zoom factor
slope = np.arctan(zf*np.sqrt(x*x + y*y))复制成功

绘制坡度图：

 代码块
Python
自动换行
复制代码
axesm()
geoshow(&#39;country&#39;, edgecolor=&#39;k&#39;)
levs = arange(0.01, 0.61, 0.01)
imshow(lon, lat, slope, levs, cmap=&#39;MPL_Greys&#39;)复制成功

高程放大系数 zf = 5 的坡度图：

坡向用 arctan2 函数计算：

 代码块
Python
自动换行
复制代码
aspect = np.arctan2(-x, y)复制成功

坡向图：

计算山体阴影需要设置太阳高度角和方位角：

 代码块
Python
自动换行
复制代码
altitude = np.deg2rad(30)
azimuth = np.deg2rad(30)
zenith = np.pi/4 - altitude

shaded = 255.0 * ((cos(zenith)* cos(slope)) +
         (sin(zenith) * sin(slope) * cos(azimuth- aspect)))
shaded[shaded&lt;0] = 0复制成功

山体阴影绘图：

 代码块
Python
自动换行
复制代码
levs = arange(200., 255, 1)
imshow(lon, lat, shaded, levs, cmap=&#39;MPL_Greys_r&#39;)复制成功

叠加上带有一定透明度的高程颜色，绘制出有立体感的高程变化图：

 代码块
Python
自动换行
复制代码
imshow(lon, lat, data, 40, cmap=&#39;MPL_terrain&#39;, alpha=0.6, zorder=1)复制成功

完整的代码如下：

 代码块
Python
自动换行
复制代码
fn = &#39;D:/Temp/nc/etopo2_new.nc&#39;
f = addfile(fn)
data = f[&#39;z&#39;][&#39;25:50&#39;,&#39;72:115&#39;]
data = data.astype(&#39;float&#39;)
lon = data.dimvalue(1)
lat = data.dimvalue(0)
x, y = np.gradient(data, 3000, 3000)
zf = 5.    #slope zoom factor
slope = np.arctan(zf*np.sqrt(x*x + y*y))

aspect = np.arctan2(-x, y)

altitude = np.deg2rad(30)
azimuth = np.deg2rad(30)
zenith = np.pi/4 - altitude

shaded = 255.0 * ((cos(zenith)* cos(slope)) +
         (sin(zenith) * sin(slope) * cos(azimuth- aspect)))
shaded[shaded&lt;0] = 0

#Plot
axesm()
geoshow(&#39;country&#39;, edgecolor=&#39;k&#39;)
#imshow(lon, lat, aspect, 40, cmap=&#39;MPL_Greys&#39;)
#levs = arange(0.01, 0.61, 0.01)
#imshow(lon, lat, slope, levs, cmap=&#39;MPL_Greys&#39;)
levs = arange(200., 255, 1)
imshow(lon, lat, shaded, levs, cmap=&#39;MPL_Greys_r&#39;)
imshow(lon, lat, data, 40, cmap=&#39;MPL_terrain&#39;, alpha=0.6, zorder=1)
colorbar()复制成功

还有一种更简单的方式估算山体阴影，不用计算坡度、坡向等，效果也不错：

 代码块
Python
自动换行
复制代码
fn = &#39;D:/Temp/nc/etopo2_new.nc&#39;
f = addfile(fn)
data = f[&#39;z&#39;][&#39;25:50&#39;,&#39;72:115&#39;]
lon = data.dimvalue(1)
lat = data.dimvalue(0)
x, y = np.gradient(data)
shaded = x * 0.7 + y * 0.5

#Plot
axesm()
geoshow(&#39;country&#39;, edgecolor=&#39;k&#39;)
geoshow(&#39;cn_province&#39;)
levs = linspace(-100, 100, 40)
imshow(lon, lat, shaded, levs, cmap=&#39;MPL_gist_yarg_r&#39;)
imshow(lon, lat, data, 40, cmap=&#39;MPL_terrain&#39;, alpha=0.6, zorder=2)
colorbar()复制成功

本文为我原创，未经授权禁止转载

cv34857440

分享至

投诉或建议