基于MATLAB等差分纬线多圆锥投影变换.docVIP

基于MATLAB等差分纬线多圆锥投影变换 　　摘 要：在MATLAB平台下，综合利用等差分纬线多圆锥投影的正反解变换实现了投影后矢量和栅格数据的输出。输出后的成图不仅保证了数据的准确性，可靠性，并且可以供用户直接在地理信息系统软件ArcMap中读取和编辑，方便用户直接使用。 关键词：等差分纬线多圆锥投影 MATLAB 投影变换 TFW 中图分类号：P226 文献标识码：A 文章编号：1672-3791（2013）06（b）-0039-03 等差分纬线多圆锥投影是一种属于任意性质的多圆锥投影，由我国地图工作者于1963年根据我国同周围的国家和地区的位置和形状关系自主设计的一种投影。使用这种投影制作世界地图时，往往选取东经150°作为中央经线，这样不仅能完整显示太平洋和各洲大陆，还能使我国版图处于图上的相对居中位置[1，2]。因此该投影在中国大陆出版的各种世界地图专题图中被广泛应用，取得了较好的效果。 等差分纬线多圆锥投影的性质为：投影后的纬线（包括极点）为同轴圆弧，每个圆弧的圆心都位于中央经线上；中央经线的投影为一直线，其他经线的投影为曲线，并与中央经线对称，且离中央经线越远，其经线的间隔越成比例地递减；极点投影长度为赤道投影长度的一半。 由于等差分纬线多圆锥投影没有给出直接的投影公式，所以用户在使用该投影时往往需要自己来设置一些参考点才能够进行曲线拟合[3]，求解出该投影的正反解变换公式。 MATLAB是一种集数学计算功能、图形化显示功能等多种功能于一体的高级科学计算软件，能够高效快速地解决各种数学相关的科学问题，因此MATLAB被广泛应用于不同领域[4]。本文将在MATLAB平台下进行等差分纬线多圆锥投影的正反解变换，实现该投影下世界地图的矢量和栅格数据输出，以使这些数据可以直接在地理信息系统软件ArcMap中读取和编辑，方便用户直接使用。 1 等差分纬线多圆锥投影变换在MATLAB中的实现 1.1 投影变换使用的地图源数据 （1）世界地图的.shp格式文件（见图1）。 （2）像元行列比为1∶2的栅格世界地图（8192×4096像素）（见图2）。 1.2 投影变换使用的函数文件 以下均为在MATLAB软件平台下进行投影变换所用到的公式。这些公式可以在MATLAB软件中各自存储成一个.m类型的函数文件，以方便其他文件调用。 （1）正解变换公式：（xi，yi）=map_forwar d（long，lati） 式中，long为经度，范围在-180～180之间；lati为纬度，范围在-90-90之间；在经度值为long和纬度值为lati的条件下，xi，yi分别为经等差分纬线多圆锥投影正解变换[1，3，5，6]求得的x值和y值。若long和lati的值超过一定范围时，xi和yi均返回NaN值。 （2）反解变换公式：（long，lati）=map_inverse（xi，yi） 式中各个变量的定义和前面第2节（1）相同，在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的最大正整数。 该式通过指定的经纬度转换成栅格数据所对应的像元的行列号，把该行列号所在的像元值赋值给指定区域。因此该式主要用于栅格数据等差分纬线多圆锥投影的变换。 1.3 投影变换在MATLAB中的实现 （1）矢量数据的投影变换。 以前面所述的.shp矢量世界地图为例，设使用的矢量世界地图文件名为“world _map.shp”，读取、转换、显示和输出中央经线为东经150°的等差分纬线多圆锥投影世界地图代码如下[7]： S=shaperead（’world_map.shp’）； %读取shapefile文件。 Slength=length（S）； %读取属性长度。 for cou=1：Slength xlength=length（S（cou，1）.X）； for count=1：xlength long=S（cou，1）.X（count）；lati=S（cou，1）.Y（count）； %读取某一属性的经纬度坐标。