桥接桥梁（R 到 Python 到 IDL）
原文链接：https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/bridge-to-bridge-r-to-python-to-idl
19045 对此文章评分：
暂无评分
桥接桥梁（R 到 Python 到 IDL）
匿名作者 2016年2月11日 星期四
在 ENVI 和 IDL 中，有许多不同的分类方案和函数可用于辅助分析数据（ENVI 分类文档）。然而，有时您可能会注意到我们尚未有机会将您最喜欢的方法编码进去，或者您希望使用的方法非常前沿，略微领先于我们的开发周期。如果是这种情况，那么您想要的方法可能在不断扩展的 R 统计软件包中找到。正如您上周在 IDL Data Point 博客上可能读到的，我们现在可以通过使用 IDL 到 Python 桥接来调用 R！本周，我想向您展示另一个如何实现此功能的示例，以及一些帮助您入门的技巧。
如果您尚未按照步骤设置您的 IDL、Python 和 R 环境，那么现在最好完成这些设置。一切就绪后，我们就可以深入探讨本文的核心内容。本文的目的不是详细探讨这些桥接的所有细微差别，而是通过一个有趣的示例展示如何可视化分类树，并在此过程中指出一些潜在问题。
首先，我只是加载了一些数据，对其进行了子集化，并对数据进行了手动分类。我选择了4个不同的类别：非光合植被（NPV）、植被（VEG）、城市（Urban）和水体（Water）。
;启动 ENVI
e = envi(/current)
if ~OBJ_VALID(e) then e = envi()
; 打开测试图像
file1 = FILEPATH('qb_boulder_msi', ROOT_DIR=e.ROOT_DIR, $
SUBDIRECTORY = ['data'])
oRaster = e.OpenRaster(file1)
; 对图像进行子集化
oSubSet = oRaster.Subset(Sub_Rect = [207,705,420,838])
; 获取图像的偏移量
xo = oSubSet.METADATA['X START']
yo = oSubSet.METADATA['Y START']
; 提取 ROI 数据
Water = oSubSet.GetData(Sub_Rect = [300-xo,800-yo,312-xo,808-yo], interleave='bip')
NPV = oSubSet.GetData(Sub_Rect = [242-xo,758-yo,249-xo,765-yo], interleave='bip')
Urban1 = oSubSet.GetData(Sub_Rect = [321-xo,736-yo,326-xo,740-yo], interleave='bip')
Urban2 = oSubSet.GetData(Sub_Rect = [363-xo,713-yo,366-xo,718-yo], interleave='bip')
Veg = oSubSet.GetData(Sub_Rect = [369-xo,738-yo,372-xo,742-yo], interleave='bip')
roi_data = list(Water,NPV,Urban1,Urban2,Veg)
然后，我以更易于导入 R 环境的方式构建了数据。
; 获取输出数组的大小
npts = 0
for i = 0 , n_elements(roi_data)-1 do npts = npts + product((size(roi_data[i], /DIMENSIONS))[1:2])
training_data = intarr(oSubSet.nb,npts)
class = strarr(npts)
classes = ['Water', 'NPV', 'Urban', 'Urban', 'Veg']
; 构建训练数据
count = 0
for i = 0 , n_elements(roi_data)-1 do begin
training_data[*, count : count + product((size(roi_data[i], /DIMENSIONS))[1:2])-1] = reform(roi_data[i], oSubSet.nb, product((size(roi_data[i], /DIMENSIONS))[1:2]))
class[count : count + product((size(roi_data[i], /DIMENSIONS))[1:2])-1] = classes[i]
count = count + product((size(roi_data[i], /DIMENSIONS))[1:2])
endfor
; 收集输入信息
classes = strjoin(class,",",/SINGLE)
band1 = training_data[0,*]
band2 = training_data[1,*]
band3 = training_data[2,*]
band4 = training_data[3,*]
接下来是棘手的部分。我从 Python 内部调用 rpy2 对象，并将 robject 返回到 IDL。
; 为新的 R 数据框准备好 Python
!null = Python.run('import rpy2.robjects as robjects')
robjects = Python.robjects
一旦我在 Python 中准备好使用 R，我就在 IDL 中动态编写了一些简单的 R 代码，以便可视化构成我手动图像分类的不同波段阈值。
; 定义 R 函数
train_rf = "robjects.r('''" + $
"train_rf <- function(band1, band2, band3, band4, classes) { \n" + $
"require(rpart) \n" + $
"classes = c(unlist(strsplit(classes,','))) \n" + $
"train_df = data.frame(band1, band2, band3, band4, classes) \n" + $
"fit = rpart(classes ~ band1 + band2 + band3 + band4, method='class', data=train_df) \n" + $
"plot(fit, uniform = TRUE, main = 'Classification Tree for Training Data') \n" + $
"text(fit, use.n=TRUE, all=TRUE, cex=.8) \n" + $
"}" + $
"''')"
这里一个可能的障碍是，R 必须安装有 rpart 包，此代码才能运行。我将其安装在主 R 库目录中，以确保我机器上的任何用户都能找到这个包（例如，C:\Program Files\R\R-3.1.3\library）。只要安装了这个包，代码中的 "require(rpart)" 部分应该不会给您带来任何麻烦。
最后，您需要将新定义的 train_rf 函数发布到 Python。
!null = Python.run(train_rf)
一旦它被 Python 识别，您就可以将其引入 IDL 并像原生函数一样运行它。
train_rf = robjects.globalenv['train_rf']
; 从 IDL 内部调用 R 函数
result = train_rf(band1, band2,Band3, band4, classes)
运行代码时，它应该生成类似这样的内容！
