IDL编程分块读取影像_张卓宇.pdf

2023.4电脑编程技巧与维护1概述遥感是通过主动或被动方式探测地物电磁辐射特性信息来表征的。近年来，随着天地一体化对地观测系统和智能计算的快速发展，遥感技术正在进入大数据时代1。遥感信息提取是遥感研究的重要分支之一，遥感影像提取方式有两大类：一类是目视解译，即人工通过观察影像的色调、大小、形状、纹理等完成对地物范围进行圈定；另一类是计算机辅助影像提取，它的传统做法是采用自动或半自动的提取方式获取地物的外部轮廓。IDL是一种数据分析和图像化的应用程序及编程语言，由美国ExelisVis公司所有。IDL为用户提供了可视化的数据分析解决方案，在地球科学、医学影像、图像处理、软件开发等多个领域有着广泛的应用。图像处理中经常用到的ENVI是使用IDL编写的一个经典软件2。开发“环南极冰山遥感影像提取算法”时，一般用IDL封装的函数read_tiff读取整个数组data，然后通过size（data/dimensions）封装函数获取影像的行数与列数，然而当读取的影像过大时，往往需要花大量的时间读取整幅影像，常会出现内存不足的情况。查阅了很多相关的开发文献及资料，都没有找到解决的办法。在反复查阅IDL官方开发帮助文档和不断测试后，终于找到分块读取遥感影像的函数，并实验解决了“快速获取大型遥感影像行列号”的难点，用IDL实现了对行列号的计算。2数据合成孔径雷达（SAR）图像已成为冰山探测的主要数据来源，即使在极夜也能穿透云层进行高度的空间和时间监控3。Sentinel-1是一种合成孔径雷达任务，支持在不受云层遮掩或缺少光照的波长下工作，并且支持在白天或晚上通过站点获取数据资料。提供C波段连续的全天候昼夜图像（中心频率：5.405 GHz），4路运行独特的成像模式，具有不同的空间分辨率和覆盖范围，由共享同一轨道平面的两颗卫星Sentinel-1A和Sentinel-1B组成，单星的重访期限为12天，两颗卫星组成星座协同工作4。环南极冰山监测基于Google Earth Engine平台5合成的20152020年8月距离1号时间差最小的HH极化方式的Sentinel-1数据集。Google Earth Engine是一种地理空间处理服务，可以在谷歌云平台的支持下进行大规模、大范围的地理空间处理，可为大规格地理空间算法开发提供交互式平台，达到高主导、高数据带动的科学，在牵涉庞大地理空间数据集的全球考验方面取得实际性进展6。Google Earth Engine平台展示如图1所示。Google Earth Engine平台中Sentinel-1图像采集系统中的数据由一级地面距离探测（GRD）影像组成，经过SNAP Toolbox预处理，具体步骤包括：轨道文件应用、GRD影像边缘去除、热噪声去除、辐射定标和地形正射校正，最终处理为以分贝（dB）为单位的后向散射系数()。将SNAP预处理后的HH极化影像分块导出为Geotiff格式影像，最终在ArcMap10.7中拼接和重投影为极投影坐标影像。20152020年8月间合成影作者简介：张卓宇（1997），女，助理工程师，硕士，研究方向为遥感数据处理。IDL 编程分块读取影像张卓宇（西安煤航遥感信息有限公司，西安710100）摘要：在开发“环南极冰山遥感影像提取算法”时，发现遥感影像需要分块进行图像增强，而 ENVI+IDL 软件仅能通过读取整幅影像获取影像的行列号。当遥感影像像幅过大时，往往会出现读取失败或耗时较长的问题，并且目前可查阅的相关文献及资料中都没有找到解决方案。介绍了基于交互式数据语言 IDL 影像的行列号快速获取方式，并解决了影像分块读取问题，给出相关程序代码。关键词：遥感；冰山；IDL 技术；行列号；分块读取图1Google Engine Earth平台展示147DOI:10.16184/prg.2023.04.0112023.4电脑编程技巧与维护像，如图2所示。3解决思路及代码实现3.1求解行号计算行数的思路：根据设定的尺长参数n，读取尺长长度的tif影像，结合catch算法，捕获错误状态并直到尺长小于等于行数时结束循环，获取适合该影像的初始尺长参数。参数l为累计尺长-1，count为当前尺长的累计尺数，当尺长大于等于1时开始循环。分块读取0,l,1,1位置的tif影像失败即下标超限时，尺长n变为原来的一半n/2，累计尺长l变为原始累计尺长l-尺长n；当下标未超限时，读取0,l+n,1,1位置的下标数据，累计尺数+1，累计尺长更新为原始累计尺长l+现有尺长n。;求tif文件的行数，fn为tif路径，ruler为尺长function tifnl,fn,rulern=rulerl=n-1CATCH,Error_status;捕捉错误状态repeat beginIf(Error_status)ne 0 then begin;当前下标读取超限时n=n/2l=n-1temp_0=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=0,l,1,1)CATCH,Error_statusendif else beginl=ntemp_0=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=0,l,1,1)endelseendrep until Error_status eq 0;以上代码是寻找匹配图大小的尺长;以下统计行数n=n;尺长l=long(n-1);累计总长度-1count=long(0);累计不同尺长下的尺数while(n ge 1)do begin;当尺长大于等于1时开始循环Temp 0=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=0,l,1,1)CATCH,Error_statusif(Error_status)ne 0 then begin;当前下标读取超限时count=0n=n/2l=l-nendif else if(Error_status)eq 0 then begin;当前下;标读取未超限时temp_1=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=0,l+n,1,1)count+=1l=l+nendifendwhile;结束循环return,l+1end3.2求解列号计算列数的思路：根据设定的尺长参数n，读取尺长长度的tif影像，结合catch算法，捕获错误状态并直到尺长小于等于列数时结束循环，获取适合该影像的初始尺长参数。参数l为累计尺长-1，count为当前尺长的累计尺数，当尺长大于等于1时开始循环。当分块读取l,0,1,1位置的tif影像失败即下标超限时，尺长n变为原来的一半n/2，累计尺长l变为原始累计尺长l-尺长n；当下标未超限时，读取l+n,0,1,1位置的下标数据，累计尺数+1，累计尺长更新为原始累计尺长l+现有尺长n。function tifns,fn,rulern=rulerl=n-1CATCH,Error_statusrepeat beginif(Error_status)ne 0 then begin;当前下标读取超限时n=n/2l=n-1temp_0=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=l,0,1,1)CATCH,Error_statusendif else beginl=ntemp_0=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=l,0,1,1)endelseendrep until Error_status eq 0;以上代码是寻找匹配图大小的尺度;以下为统计列数图220152020年8月间合成影像2015201620172018201920201482023.4电脑编程技巧与维护n=n;尺长l=long(n-1);累计总长度count=long(0);累计不同尺度下的尺数while(n ge 1)do begin temp_0=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=l,0,1,1)CATCH,Error_status;设置错误捕获主要是为了捕获下标超限的情况if(Error_status)ne 0 then begin;当前下标读取超限时count=0n=n/2l=l-nendif else if(Error_status)eq 0 then begin;当前下;标读取未超限时temp_1=READ_TIFF(fn,geotiff=getinfo,SUB_RECT=l+n,0,1,1)count+=1l=l+nendifendwhilereturn,l+1End3.3分块滑动读取数据分块滑动读取遥感影像的步骤如下：窗口行数或窗口列数小于0或者大于影像行数或列数时，不读取影像；分块读取影像时，计算行块数和列块数；滑块未到达最右侧、最下侧和右下角时，或到达最右侧、最下侧和右下角时记录每个滑块4个角点坐标，其中，若滑块在边缘超出影像时采用反向取块的方法。;ns为影像的列数;nl为影像的行数;wds为窗口列数;wdl为窗口行数;steps为列方向滑动步长;stepl为行方向滑动步长;dim为扫描方式，1为沿着列滑动(上下)，2为沿着行(左;右)滑动;为使程序更加简洁，最右边和最下边采用反向取块function move_block,ns,nl,w ds,wdl,steps,stepl,dim=dimif wds gt ns or wdl gt nl then beginprint,WindowSize ge imageSizereturn,0endifif wds le 0 or wdl le 0 then beginprint,WindowSize should ge 0;若不是行或者列扫描，则进行非滑动endifif(0 ge(steps0;左上列坐标slide_data1,num=wdl*j0slide_data2,num=wds*(i+1)-1;右下列坐标slide_data3,num=wdl*(j+1)-1max_ns=wds*(i+1)-1max_nl=wdl*(j+1)-1;当滑块到右端，还未到下端if max_ns ge ns and max_nl lt nl then beginslide_data0,num=ns-wds;左上列坐标slide_data1,num=wdl*j0slide_data2,num=ns-1;右下列坐标slide_data3,num=wdl*(j+1)-1endif;当滑块到下端，还未到右端if max_ns lt ns and max_nl ge nl then beginslide_data0,num=wds*i0;左上列坐标slide_data1,num=nl-wdlslide_data2,num=wds*(i+1)-1;右下列坐标slide_data3,num=nl-1endif;当滑块到右下角if max_ns ge ns and max_nl ge nl then beginslide_data0,num=ns-wdsslide_data1,num=nl-wdsslide_data2,num=ns-1slide_data3,num=nl-1endifnum=num+1endfor1492023.4电脑编程技巧与维护endforendif else begin;取整用ceil是防止滑动时错过一些数据case 1 of;上下滑动dim eq 1:beginsChunk=ceil(float(ns-wds)/steps)+1lChunk=ceil(float(nl-wdl)/stepl)+1;列行方向窗;口个数k=sChunk*lChunk;总块数slide_data=lonarr(4,k)num=0;计数器，逐行写入数组for i=0,(sChunk*steps-1),steps do begin;沿列;方向进行滑动,列之间的窗口有重叠-1是为了防止i;溢出列数for j