李国春

随着高分辨率对地观测卫星发射的日益增多，对数据处理软件的要求也越来越高。通常每个系列卫星都有自己的数据特点并需要专门的处理软件，但卫星数量的增加为每种卫星单独设计软件的压力越来越大。本文介绍的一种处理方案旨在能够正射校正处理大多数的光学遥感卫星的观测数据。

将卫星运营部门发布的卫星数据处理成用户数据最重要的处理成本是正射校正部分。这里提供的处理方案可以快速无人工干预处理由有理式（RPC）参数描述的一级数据。

本文演示处理了几种卫星数据，更多类别的卫星（尤其是近年大量发射的新卫星）数据用户可以用本文程序试验处理，遇有问题可以及时反馈以便修正处理。

一、99脚本语言处理代码

1. 准备正射校正用DEM数据

	Print("开始生成中国区 SRTM 90m DEM 数据...");
	DWORD tc = GetTickCount();
	//创建 DEM 数据，正射校正用///
	double		demGeoBox[4];				//排列依次为：最小纬度、最小经度、最大纬度和最大经度
	longlong	demHeight = 36000;
	longlong	demWidth = 42000;
	short		demBuffer[demHeight][demWidth];
	double		dbOrgLat = 46.0;			//设置左上角位置的经纬度
	double		dbOrgLon = 95.0;
	short		defaltElv = 0.0;	
	double		demPixSize = 5.0/6000.0;	//SRTM90m = 0.0008333。ASTGTM2、3 DEM 为 0.00027778度/像元.pixsize = 1.0/3600.0;
	DemSrtm90Buffer(demBuffer,dbOrgLat,dbOrgLon,demHeight,demWidth,defaltElv,demGeoBox);
	double dmtm = (GetTickCount() - tc) / 1000.0;
	Print("生成中国区 SRTM 90m DEM 数据用时 %f 秒。",dmtm);

这段代码从SRTM90读来高程数据，存放在DemSrtm90Buffer指向的缓存。数据左上角纬度46.0°N，经度95°E。东西向42000个像元，南北向36000个像元。每个像元5/6000°。这样根据经纬度位置能够定位到对应的高程。

正射校正用DEM数据可以使用SRTM90、ASTGDEM2、3等通用DEM，也可以使用用户自行定义的DEM数据，只需指定开始经纬度位置、像元大小和数据范围，并将数据预存到一个数据缓存中即可。

2. 选择待校正的数据集

	//选择待进行正射校正的卫星数据文件，自动查找对应的RPC数据文件 ///
	longlong	width,height;
	int			bands,datatype;
	STRING tifName = OpenFileDialog(TRUE,"*.tif;*.tiff");//打开一个tif文件
	STRING	rpcName = tifName - '.';
	rpcName = rpcName + ".rpb";
	if(!PathFileExists(rpcName))
	{
		rpcName = tifName - '.';
		rpcName = rpcName + "_rpc.txt";
		if(!PathFileExists(rpcName))
		{
			Print("没有找到多光谱数据 %s 对应的 RPC 文件。",tifName);
			return FALSE;
		}
	}

这段代码用OpenFileDialog函数通过对话框选择一个待处理的文件tifName，并通过文件名查找对应的有理式系数（RPC）文件rpcName。先查找.rpb文件，如果该文件不存在再查找.txt文件。

3. 选择结果输出目录并生成一个输出文件名

	//选择一个结果文件的输出路径，并指定导出文件的类型（.img或者.tif）/
	STRING outPath;
	BrowseFolder(outPath);
	STRING mainName = FileNameMain(tifName);
	STRING outName = outPath + mainName + ".tif";//正射校正后数据的导出文件名

输出文件名指定.tif扩展名自动保存为TIFF文件，指定为.img时自动保存为带.hdr的IMG格式。

4. 开始计时

	Print("开始对%s进行正射校正...",tifName);//开始进行正射校正/
	DWORD tc = GetTickCount();

5. 获取数据行数、列数、波段数、数据类型和四角经纬度信息

	//检测数据集基本信息///
	longlong	width,height;
	int			bands,datatype;
	double quad[8]; //数据四角经纬度，每角一对，书写顺序
	int b = TifInqDataInfo(tifName,height,width,bands,datatype);	//读取 TIFF 文件的数据范围、数据类型和波段数
	int bo2 = SwRpcGeoQuad(rpcName,quad,height,width,demBuffer,demHeight,demWidth,SINT16BL,demGeoBox,demPixSize);//四角坐标

6. 指定投影参数、四角经纬度转换成Northing、Easting坐标

	STRING projstr = "";
	double box[4];
	ProjGeoCoordTransForEx(projstr,2,2,quad,box);//四角经纬度转换成投影坐标 N、E

投影参数是RSD投影串，格式为：

投影串定义

/ ID = <投影ID>, //投影ID，例如UTM、ALBERS_CONIC_EQAREA等

/ a = <椭球体半长轴>, //椭球体半长轴，如WGS84 为6378137

/ f = <椭球体扁率>, //椭球体扁率(实为扁率的倒数1/f),如WGS84 为298.2572,

/ SP1 = <第1标准纬线>, //圆锥投影的第1标准纬线，麦卡托的割线，UTM、高斯投影的带号。

/ SP2 = <第2标准纬线>, //圆锥投影的第2标准纬线。

/ org_lat = <原点纬度>, //投影的纬度原点。

/ org_lon = <中央经线>, //投影的中央经线。

/ FN = <伪北>, //False Northing，伪北坐标。

/ FE = <伪东>, //False Easting，伪东坐标。

/ K = <系数>, //变换系数。

/ DATUM = <大地基准>, //缺省为WGS84。

这里STRING projstr = "";定义了一个空串，为缺省的UTM投影，投影参数由数据范围自动探测得到。

7. 估算像元尺度和正射后图像大小

	//估算像元大小，南北方向和东西方向都估算，取平均值，然后取整。（或者根据卫星已知信息直接指定像元尺度）
	double nsize,esize;
	double pixsize = SwEvalPixelSize(quad,height,width,nsize,esize);	//估算像元尺度的函数

	//正射校正后的图像大小/
	height = (box[2]-box[0]) / pixsize;
	width = (box[3]-box[1]) / pixsize;

8. 正射校正

	//正射校正，///
	int precis = TRUE;	//精度选项，TRUE时全程使用 8 字节浮点精度。	
	int resamp = CC;//重采样像元重定位方式 FW=前向重采样，NN=最近邻点法，BL=二次线性插值法，CC=三次褶积法。
	WORD retBuffer[height][width][bands]; //哑数组，接收处理结果
	int boo = Sw7WristbandsBuffer3(tifName,rpcName,retBuffer,box,projstr,pixsize,demBuffer,demGeoBox,demPixSize,precis,resamp);

正射校正结果保存在retBuffer缓存中。

9. 统计处理时间

	double tm = (GetTickCount() - tc) / 1000.0;
	Print("正射校正结束，并保存在数组retBuffer中。总计用时 %f 秒。",tm);

10. 输出结果、结束程序

	//正射校正结果导出为TIFF数据格式文件。
	Print("开始导出正射校正结果...");
	tc = GetTickCount();
	int rank,dims,datatype;
	int dims[16];
	InqVariable(retBuffer,rank,dims,datatype);//读正射校正后PAN数据的数据范围、维数和数据类型信息
	height = dims[0];
	width = dims[1];
	bands = dims[2];
	GdSaveTiffFile(outName,retBuffer,height,width,bands,datatype,BIP,projstr,box,pixsize);
	double outpuTime = (GetTickCount() - tc) / 1000.0;
	Print("数据保存为 %s 用时 %f 秒。",outName,outpuTime);

	Print("主程序结束");
	return 1;