基于MATLAB的等差分纬线多圆锥投影变换Word下载.docx

1、MATLAB是一种集数学计算功能、图形化显示功能等多种功能于一体的高级科学计算软件，能够高效快速地解决各种数学相关的科学问题，因此MATLAB被广泛应用于不同领域 4。本文将在MATLAB平台下进行等差分纬线多圆锥投影的正反解变换，实现该投影下世界地图的矢量和栅格数据输出，以使这些数据可以直接在地理信息系统软件ArcMap中读取和编辑，方便用户直接使用。2 等差分纬线多圆锥投影变换在MATLAB中的实现2.1 投影变换使用的地图源数据1） 世界地图的.shp格式文件（见图1）；2） 像元行列比为1：2的栅格世界地图（81924096像素）（见图2）。2.2 投影变换使用的函数文件以下均为在MA

2、TLAB软件平台下进行投影变换所用到的公式。这些公式可以在MATLAB软件中各自存储成一个.m类型的函数文件，以方便其他文件调用。（1）正解变换公式：（xi,yi）=map_forward（long,lati）式中，long为经度，范围在-180180之间；lati为纬度，范围在-9090之间；在经度值为long和纬度值为lati的条件下，xi，yi分别为经等差分纬线多圆锥投影正解变换1,3,5,6求得的x值和y值。若long和lati的值超过一定范围时，xi和yi均返回NaN值。（2）反解变换公式：（long,lati）=map_inverse（xi,yi）式中各个变量的定义和前面第2节（1

3、）相同，在x、y坐标值分别为xi和yi的条件下long和lati分别为经等差分纬线多圆锥投影反解变换1,6求得的经度值和纬度值。若xi和yi的值超过投影区域范围时，long和lati均返回NaN值。（3）经纬度转像元行列号的公式：（x,y）=image_position（long,lati,image）式中long和lati的定义和前面第2节（1）相同，image为读入的图像矩阵，在该函数的内部读取该图像的总行数imx和总列数imy。像元行列号x和y分别与long和lati的换算关系如下：（1）（2）式中，x代表不超过x的最大正整数。该式通过指定的经纬度转换成栅格数据所对应的像元的行列号，把该

4、行列号所在的像元值赋值给指定区域。因此该式主要用于栅格数据等差分纬线多圆锥投影的变换。2.3 投影变换在MATLAB中的实现（1）矢量数据的投影变换以前面所述的.shp矢量世界地图为例，设使用的矢量世界地图文件名为“world_map.shp”，读取、转换、显示和输出中央经线为东经150的等差分纬线多圆锥投影世界地图代码如下7：S=shaperead（´world_map.shp&）; %读取shapefile文件。Slength=length（S）;%读取属性长度。for cou=1:Slengthxlength=length（S（cou,1）.X）;for count=

5、1:xlengthlong=S（cou,1）.X（count）;lati=S（cou,1）.Y（count）; %读取某一属性的经纬度坐标。if long<=-30%此条件语句为输出中央经线为东经150的投影时long=long+210; %的经度变换。elselong=long-150;endxi,yi=map_forward（long,lati）; %遍历该Shapefile文件上的每一点，通过逐S（cou,1）.X（count）=yi; %点变换，最终生成该投影下的地图数据。S（cou,1）.Y（count）=xi;mapshow（S）;%显示投影变换后的地图数据。shapewri

6、te（S,&world_map_conversion.shp& %导出地图数据。在以上代码中，第二行执行后的S变量实际上是一个数组，且数组中的每一个元素都是结构体，上面显示了几何类型、属性名称、X坐标和Y坐标等信息。通过嵌套循环语句逐个访问了每一个点的经纬度坐标，并将变换后的X、Y坐标逐个赋值给S中原来的经纬度坐标。变换后还输出并保存为Shapefile格式文件。变换后的结果如图3所示。图3 等差分纬线多圆锥投影世界地图矢量数据（2）栅格数据的投影变换以前面所述的栅格世界地图数据为例，设使用的栅格世界地图文件名为“world_map.jpg”，输出分辨率为50004000的中央经线为东经150

7、的等差分纬线多圆锥投影世界地图。实现这一操作的代码如下：Mapc=uint8（ones（4000,5000,3）; %建立三维矩阵，图像矩阵类型为uint8。imagec=imread（&world_map.jpg&%读取栅格世界地图文件。Mapc=Mapc*255; %使Mapc输出的图片背景为白色。p,q,r=size（Mapc）;%读取Mapc矩阵的尺寸。for i=1:pxi=（i-2000）*7000; %第6行和第8行：通过一系列变换，使xi和yifor j=1:q%的范围能包含投影范围。yi=（j-2500）*7000;long,lati=map_inverse（xi,yi）;

8、%投影反解变换。if isnan（long）&isnan（lati） %若long和lati的值为NaN，不进行颜色赋值，continue; %并且直接跳入下一轮循环。=30&long>=-180 %此条件语句为输出中央经线为东经150的投影。long=long+150;=180long=long-210;row,column=image_position（long,lati,imagec）; %读取经纬度对应栅格世界for k=1:r %地图的行列号。Mapc（i,j,k）=imagec（row,column,k）; %把世界地图的任意一点的end %像元值赋值给投影地图的end %对

9、应位置上。image（Mapc）; %显示地图图像。该代码首先是建立50004000的空白图像，然后通过嵌套循环语句逐点扫描该空白图像，把每一个像元的行列号（i,j）经一系列的变换转化成该投影下的x、y坐标，利用反解变换求解对应的经纬度坐标，之后读取经纬度坐标所对应的栅格世界地图图像所对应的行列号（row,column），把世界地图中行列号为（row,column）的像元值赋值给50004000图像的行列号为（i,j）的像元。当循环语句结束后，显示投影变换后的图像并保存。由于保存后的图像已缩小7000倍，这样如果在ArcMap中直接打开该图像，会导致数据显示不正确。但“GeoTiff”解决了这

10、一瓶颈。tfw文件是TIFF文件坐标信息的文本文件，ArcInfo、MicroStation、AutoCAD等软件均支持该格式的坐标信息文件。此文件定义了栅格图像素坐标与实际大地坐标的仿射关系8，9。一个.tfw格式的文件构成如下10：第1行：A：地理x坐标中的像元分辨率第2行：D：y轴旋转系数第3行：B：x轴旋转系数第4行：E：地理y坐标中的像元分辨率第5行：C：左上角第1个像元中心的x坐标值第6行：F：左上角第1个像元中心的y坐标值各个参数的图解由此保存分辨率为50004000的投影后的世界地图为TIF格式，用记事本编写对应的tfw格式文件。最后在ArcMap中加载这个世界地图。3 结语本

11、文依托MATLAB软件在图形处理方面的优势，综合利用等差分纬线多圆锥投影的正反解变换公式，还利用了像素长宽比为2:1的世界地图原始数据中经纬度与行列号的关系，实现了矢量数据和栅格数据的等差分纬线多圆锥投影变换。此外本文以桌面GIS软件ArcMap为例，使用了.tfw格式文件让输出后栅格数据等差分纬线多圆锥投影世界地图准确显示。本文中的投影变换呈现出如下特点：1） 输出后的矢量和栅格地图数据可以供用户在ArcMap、MapInfo等GIS软件中被直接加载，方便用户进行编辑；2） 输出后的矢量和栅格地图数据具有准确性、可靠性。与文献3相比，整个变换过程无需人工干预，节省了人工计算配准点的繁重劳动，

12、减少了栅格地图数据地理配准中可能出现的误差。参考文献1 董曼，李胜乐世界地图等差分纬线多圆锥投影的正反解变换J大地测量与地球动力学，2008，28（2）：95-992 祝国瑞地图学M武汉：武汉大学出版社，20043 叶远智栅格数据的等差分纬线多圆锥投影转换J测绘通报，2012（9）：68-704 刘正君MATLAB科学计算宝典M北京：电子工业出版社，20125 胡毓钜地图投影M北京：测绘出版社，19816 杨启和地图投影变换原理和方法M北京：解放军出版社，19897 苏金明，王永利MATLAB图形图像M北京：电子工业出版社，20058 余建军，蔡以雷，洪景峰一种利用数字线划图快速生产数字栅格地图的方法J测绘通报，2009（7）：51-549 王宪民ArcGIS在栅格数据处理中的应用初探J测绘技术装备，2005，1（7）： 32-34+2110 吴铭杰基于AutoCAD的TIFF world（TFW）图像的应用研究A史照良第十四届华东六省一市测绘学会学术交流会论文集C南京：现代测绘编辑部，2012：157-158