前言:栅格数据可视化 RVT 软件的小小使用说明。
1.什么是 RVT
RVT 的全称是 Relief Visualization Toolbox,中文意思大概就是地形可视化工具箱(Relief这个词在地图制图方面没有一个非常准确的翻译)。
首先这个工具箱不是 ArcGIS 工具箱,是使用 IDL 语言编写的独立应用程序,非常轻量化、专一化,这是我非常喜欢的一点,出自于 Research Centre of the Slovenian Academy of Sciences and Arts(斯洛文尼亚科学与艺术学院研究中心)。
RVT 有4种主要功能:地表显示、格式转换、镶嵌、图像混合,其中地表显示涉及到 hilshade、SVF、ASVF、sky illumination、SLRM、Openness-Positive等;格式转换和镶嵌是很常见的功能;图像混合最多可以添加五张栅格数据,进行图像混合处理,用于凸显目标物体或用于其他解译。
2.为什么要小小的说明这个工具?
读者或许读到过关于
RED RELIEF IMAGE MAP (图1)的文章,其中就是使用了 RVT 来实现的;也可能了解过或者制作过
天空光照模型 (图2),天空光照模型可以使用 ArcGIS 的插件制作,但是你懂的,本来计算量就特别大,时间按小时算,最后算完了还经常出错,而 RVT 也可以制作天空光照模型,并且稳定(虽然也很慢)。
图1 RED RELIEF IMAGE MAP
图2 Sky ModelRVT 工具官网以及下载:
官网、下载:https://www.zrc-sazu.si/en/rvt
如果进不去可公众号后台回复:rvt 可获得工具箱以及说明手册。
RVT 不仅可以制作 RED RELIEF IMAGE MAP 相关的图层、天空光照模型,还能制作多方向加权叠加 Hillshade、SVF 及其相关的变种等数据。这是我使用它的另一个理由。
RVT 界面橙色选中框从上到下依次为:添加数据、可视化、天空光照模型、开始。
可视化的右边三个功能选项卡分别是:转换、镶嵌、混合。
3.使用 RVT
以制作天空光照模型为例。
你可以看到在 RVT 主界面的下半部分有个叫 Sky Illumination,这个就是天空光照模型,勾选上,然后选择天空模型(Sky model),目前内置了 overcast 和 uniform 模型。
选好后,添加 DEM 源数据,然后点击最下面的 Start 即可。
4.RVT库
软件虽好,却不能尽兴,也耗费时间,使用 RVT 自带的 python 包可以帮助到我们,自动实现流程。
4.1使用coda 安装 RVT 库
创建名为 RVT 的新环境:conda create -n RVT
激活(进入)新环境:conda activate RVT
安装 rvt 库:conda install -c zmigyyy rvt_py
查看库当前的 Python 版本:python -V Python 3.8.12
4.2使用 pypi 安装
如果没有使用 conda,可以执行 pip 命令来安装。
pip install rvt-py
RVT 在 pypi 的介绍网站:https://pypi.org/project/rvt-py/
5.使用 RVT 库
为了让处理流程更加灵活强大,比如在调用 RVT 之前对栅格数据做一些预处理,比如镶嵌、投影、裁剪或者修改像素类型,建议引入 GDAL 一起使用。
下面是个人写的流程代码,使用了 GDAL 库和 RVT 库,前者用于实现栅格的投影,后者用来快速生成 SVF 等数据产品。
代码比较简单,GDAL 使用 Warp 方法,只有一行;而 RVT 的 Python 库也基本上是开箱即用,封装程度很高,所以只是看起来长,其方法形式都是一样的。
# -*- coding:utf-8 -*-# -------------------------------------------# Name: example1# Author: Hygnic# Created on: 2021/12/30 12:35# Version: # Reference: """Description: 使用 RVT 库生成 hillshade\slope\SVF等产品Usage: """# -------------------------------------------"""--------IMPORT_SETTINGS"""import sysimport osimport rvt.vis # fro calculating visualizationsimport rvt.default # for loading/saving rasterstry: from osgeo import gdalexcept: sys.exit('ERROR: cannot find GDAL/OGR modules')print(">>"*50)print('Python: {}'.format(sys.version))print('GDAL: {}'.format(gdal.__version__))print(">>"*50)"""--------DATA"""filepath = r"E:\rasterdate"filename = os.path.join(filepath, "SRTMv3_1.tif")output_raster = os.path.join(filepath, "projected.tif")os.chdir(filepath)"""--------REPROJ"""input_raster = gdal.Open(filename)warp = gdal.Warp(output_raster,input_raster,dstSRS='EPSG:4542',resampleAlg="cubic")input_raster = None # Closes the files# CGCS2000_3_Degree_GK_CM_105E# WKID: 4544 权限: EPSG# CGCS2000 / 3-degree Gauss-Kruger CM 99E# 4542"""--------RASTER_INFO"""# dem_dataset = filenamedem_dataset = output_rasterdem_dict = rvt.default.get_raster_arr(dem_dataset)dem_arr = dem_dict["array"] # numpy array of DEMdem_resolution = dem_dict["resolution"]dem_res_x = dem_resolution[0] # resolution in X directiondem_res_y = dem_resolution[1] # resolution in Y directionsun_azimuth = 315 # Solar azimuth angle (clockwise from North) in degreessun_elevation = 45 # Solar vertical angle (above the horizon) in degrees"""--------Hillshade"""hillshade_arr = rvt.vis.hillshade( dem=dem_arr, resolution_x=dem_res_x, resolution_y=dem_res_y,sun_azimuth=sun_azimuth, sun_elevation=sun_elevation, ve_factor=1)hillshade_path = "Hillshade_315_45_f1.tif"rvt.default.save_raster( src_raster_path=dem_dataset, out_raster_path=hillshade_path, out_raster_arr=hillshade_arr,e_type=6)print("Hillshade Exported")"""--------Multiple directions hillshade"""nr_directions = 16 # Number of solar azimuth angles (clockwise from North)# (number of directions, number of bands)sun_elevation = 45 # Solar vertical angle (above the horizon) in degreesmulti_hillshade_arr = rvt.vis.multi_hillshade( dem=dem_arr, resolution_x=dem_res_x, resolution_y=dem_res_y, nr_directions=nr_directions, sun_elevation=sun_elevation, ve_factor=1)multi_hillshade_path = "MultiHS.tif"rvt.default.save_raster( src_raster_path=dem_dataset, out_raster_path=multi_hillshade_path, out_raster_arr=multi_hillshade_arr,e_type=6)print("Multi_Hillshade Exported")"""--------------------Slope"""dict_slope_aspect = rvt.vis.slope_aspect( dem=dem_arr, resolution_x=dem_res_x, resolution_y=dem_res_y, output_units="degree", ve_factor=1)slope_arr = dict_slope_aspect["slope"]slope_name = "slope.tif"rvt.default.save_raster(src_raster_path=dem_dataset, out_raster_path=slope_name, out_raster_arr=slope_arr, e_type=6)print("Slope Exported")"""--------Positive_openness & SVF"""svf_n_dir = 16 # number of directionssvf_r_max = 10 # max search radius in pixelssvf_noise = 3 # level of noise remove (0-don't remove, 1-low, 2-med, 3-high)dict_svf = rvt.vis.sky_view_factor( dem=dem_arr, resolution=dem_res_x, compute_svf=True, compute_asvf=False, compute_opns=False,svf_n_dir=svf_n_dir, svf_r_max=svf_r_max, svf_noise=svf_noise)SVF_arr = dict_svf["svf"]svf_arr_name = "SVF_denoise3.tif"rvt.default.save_raster(src_raster_path=dem_dataset, out_raster_path=svf_arr_name, out_raster_arr=SVF_arr, e_type=6)print("SVF_3 Exported")svf_noise = 1 # level of noise remove (0-don't remove, 1-low, 2-med, 3-high)dict_svf = rvt.vis.sky_view_factor( dem=dem_arr, resolution=dem_res_x, compute_svf=True, compute_asvf=False, compute_opns=False,svf_n_dir=svf_n_dir, svf_r_max=svf_r_max, svf_noise=svf_noise)SVF_arr = dict_svf["svf"]svf_arr_name = "SVF_denoise1.tif"rvt.default.save_raster(src_raster_path=dem_dataset, out_raster_path=svf_arr_name, out_raster_arr=SVF_arr, e_type=6)print("SVF_1 Exported")
最后
抛砖引玉吧,毕竟这个东西也是非常不错的,官网有更详尽的英文资料。
https://www.zrc-sazu.si/en/rvt
荟GIS精粹,关注公众号:GIS荟纯粹分享,只因热爱,你的转发是对我最大的鼓励!不然点个赞和在看也好