精华内容
下载资源
问答
  • 参数坐标转换 投影换带 大地坐标系、投影坐标系、空间直角坐标系的任意互转 程序特点: 支持任意椭球参数 程序内定义了北京54、西安80、WGS84、CGCS2000等常用椭球基准,并且可以采用自定义椭球,只需要...

    由于工作需要,编写了大地测量方面坐标转换的程序,主要实现以下功能:

    1. 投影正反算
    2. 布尔沙七参数计算
    3. 七参数坐标转换
    4. 投影换带
    5. 大地坐标系、投影坐标系、空间直角坐标系的任意互转

    程序特点:

    1. 支持任意椭球参数
      程序内定义了北京54、西安80、WGS84、CGCS2000等常用椭球基准,并且可以采用自定义椭球,只需要输入椭球长轴半径和扁率的倒数两个最常用的参数,即可完成椭球的定义。方便快捷。
    2. 支持单点转换和文件转换两种方式
      两种方式共用一个按键,两种方式同时执行。单点转换可以可视化的观察坐标转换的结果,文件转换可以大批量的进行数据转换,无数量限制。
    3. 支持多种角度单位
      无论是“度”或者dd.mmss的经纬度格式,均可支持。也可以用来进行“度”和dd.mmss格式的转换。

    经过大量的实验数据验证,程序稳定可靠。精度可达到亚毫米级别。

    程序的可执行文件可以免费分享和使用。在CSDN博客我的下载中有。
    源代码有偿提供。
    联系GeoPhoto@126.com

    软件中的基本坐标转换都是通过以下几个核心函数组合实现的。

    //度数和弧度的转换
    double Deg2Rad(double deg);
    double Rad2Deg(double rad);
    
    //度转度分秒
    double Deg2DMS(double deg);
    //度分秒转度
    double DMS2Deg(double DMS);
    
    //高斯投影正反算
    void GaussForward(double a,double f,double B,double L,double L0,double &x,double &y,double east_plus);
    void GaussBackword(double a,double f,double x,double y,double L0,double &B,double &L,double east_plus);
    
    //UTM投影正反算
    void UTMForward(double a,double f,double B,double L,double L0,double &x,double &y,double east_plus);
    void UTMBackword(double a,double f,double x,double y,double L0,double &B,double &L,double east_plus);
    
    //大地坐标系和空间直角坐标系的相互转换
    void BLH2XYZ(double a,double f,double B,double L,double H, double &X,double &Y,double &Z);
    void XYZ2BLH(double a,double f,double X,double Y,double Z, double &B,double &L,double &H);```
    
    

    以上代码每次只能转换一个点的坐标,为了实现点的批量转换,构建下面的类:
    单点类

    /******************坐标系坐标****************
    成员变量有点号、北、东、高
    可以用来表示空间直角坐标系坐标,对应为XYZ
    可以用来表示投影坐标系坐标,对应为北东高
    可以用来表示大地坐标,对应为纬度、经度、椭球高
    ********************************************/
    class Coor
    {
    public:
    	string ID;
    	double east;
    	double north;
    	double h;
    public:
    	Coor();
    	Coor(string ID, double east,double north, double h);
    	void PrintData();
    }
    

    点集类
    通过该类进行批量坐标的转换

    /******************点位坐标的集合****************
    
    ************************************************/
    class CoorVector
    {
    public:
    	vector <Coor> v_coor;
    public:
    	CoorVector();
    	//将上面的单点转换写成批量转换
    	void GaussForward_cv(double a,double f,double L0,double east_plus = 500000);
    	void GaussBackword_cv(double a,double f,double L0,double east_plus = 500000);
    	void UTMBackword_cv(double a,double f,double L0,double east_plus = 500000);
    	void UTMForward_cv(double a,double f,double L0,double east_plus = 500000);
    	void Deg2Rad_cv();
    	void Rad2Deg_cv();
    	void BLH2XYZ_cv(double a,double f);
    	void XYZ2BLH_cv(double a,double f);
    	void DMS2Deg_cv();
    	void Deg2DMS_cv();
    	//读文件
    	bool ReadData(string fileName);
    	//打印数据到屏幕
    	void PrintData();
    	//数据写入文件
    	void WriteData(string fileName);
    };
    

    到此为止,软件可以进行单点和批量的坐标转换,包括大地坐标、高斯投影、UTM投影坐标,空间直角坐标系之间的任意转换,涉及北京54、西安80、WGS84、CGCS2000和自定义坐标系。
    那么给定两组坐标,如何进行布尔沙七参数的反算呢?以及如何进行不同椭球之家附带七参数转换关系的坐标转换呢?下面一个新的类可以解决问题:
    七参数模型类:

    /*****************布尔沙七参数模型*******************
    平移参数DX,DY,DZ 单位米
    旋转参数WX,WY,WZ 单位秒
    尺度参数K,单位ppm
    初值为0,这样即使进行了七参数转换,依然不改变坐标值
    ****************************************************/
    class Brusa
    {
    public:
    	double DX;
    	double DY;
    	double DZ;
    	double WX;
    	double WY;
    	double WZ;
    	double K;
    public:
    	Brusa();
    	
    	/**************根据七参数进行坐标计算****************
    	输入为两组坐标变量,一组转换前,一组转换后
    	坐标系为空间直角坐标系
    	****************************************************/
    	void TransCoor(double X_befor,double Y_befor,double Z_befor,double &X_after,double &Y_after,double &Z_after);
    
    	/*****************布尔沙七参数反算*******************
    	平移参数DX,DY,DZ 单位米
    	旋转参数WX,WY,WZ 单位秒
    	尺度参数K,单位ppm
    	运用最小二乘计算七参数,公式为Ax=B,无需迭代
    	输入为两组点坐标,坐标系为空间直角坐标系,计算完更新
    	单位为米、弧度、单位1
    	两组点数量不一致或者数量少于3对,返回false
    	不进行点名的匹配,同名点按照顺序匹配
    	****************************************************/
    	bool CalBrusa(CoorVector vc_befor, CoorVector vc_after);
    	
    	void PrintData();
    };
    

    该类中一个函数进行七参数反算,一个函数进行使用七参数的坐标转换。下面只贴出最核心的七参数反算代码,该代码进行最小二乘计算,无需迭代。涉及矩阵运算,矩阵运算有自己写的类库,但是太水了,后面改成Engin类库。使用方法自行百度,很简单,类库源代码代码包中有。

    //参数为转换前点集、转换后点集,该函数其实就是解一个Ax=B的矩阵。
    //在代码中A=MA,x=Mx,B=MB.
    bool CalBrusa(CoorVector vc_befor, CoorVector vc_after)
    {
    		//两组坐标点数量不同,不进行计算
    		if (vc_befor.v_coor.size() != vc_after.v_coor.size())
    			return false;
    
    		//点集数量少于3组,无法计算
    		if (vc_befor.v_coor.size() < 3)
    			return false;
    
    		int num = vc_befor.v_coor.size();
    
    		MatrixXd MB(num*3,1);
    		MatrixXd MA(num*3,7);
    
    		for (int i=0;i<num;i++)
    		{
    
    			MB(i*3 + 0, 0) = vc_after.v_coor[i].north - vc_befor.v_coor[i].north;
    			MB(i*3 + 1, 0) = vc_after.v_coor[i].east - vc_befor.v_coor[i].east;
    			MB(i*3 + 2, 0) = vc_after.v_coor[i].h - vc_befor.v_coor[i].h;
    			
    			MA(i*3 + 0, 0) = 1;
    			MA(i*3 + 0, 1) = 0;
    			MA(i*3 + 0, 2) = 0;
    			MA(i*3 + 0, 3) = 0;
    			MA(i*3 + 0, 4) = 0 - vc_befor.v_coor[i].h;
    			MA(i*3 + 0, 5) = vc_befor.v_coor[i].east;
    			MA(i*3 + 0, 6) = vc_befor.v_coor[i].north;
    
    			MA(i*3 + 1, 0) = 0;
    			MA(i*3 + 1, 1) = 1;
    			MA(i*3 + 1, 2) = 0;
    			MA(i*3 + 1, 3) = vc_befor.v_coor[i].h;
    			MA(i*3 + 1, 4) = 0;
    			MA(i*3 + 1, 5) = 0 - vc_befor.v_coor[i].north;
    			MA(i*3 + 1, 6) = vc_befor.v_coor[i].east;
    
    			MA(i*3 + 2, 0) = 0;
    			MA(i*3 + 2, 1) = 0;
    			MA(i*3 + 2, 2) = 1;
    			MA(i*3 + 2, 3) = 0 - vc_befor.v_coor[i].east;
    			MA(i*3 + 2, 4) = vc_befor.v_coor[i].north;
    			MA(i*3 + 2, 5) = 0;
    			MA(i*3 + 2, 6) = vc_befor.v_coor[i].h;
    
    		}
    
    		//double **x = LeastSquare(A,num*3,7,B,num*3,1);
    		MatrixXd Mx(7,1);
    		Mx = MA.colPivHouseholderQr().solve(MB);
    		//Mx = MA.jacobiSvd(ComputeThinU|ComputeThinV).solve(MB);
    
    		DX = Mx(0,0);
    		DY = Mx(1,0);
    		DZ = Mx(2,0);
    		WX = Mx(3,0);
    		WY = Mx(4,0);
    		WZ = Mx(5,0);
    		K  = Mx(6,0);
    	}
    

    软件截图如下:
    在这里插入图片描述
    测试数据:
    下面简单介绍软件使用:
    下面三个点是WGS84坐标系下大地坐标
    坐标排列方式在软件面板上写着:
    p1 30.0 120.0 1000
    p1 30.5 120.5 1000
    p3 30.8 120.8 1000


    选取源坐标,选取源数据的坐标系,选取目标坐标系,选取写入文件,点击坐标转换按钮,即可完成计算。如果批量计算,可以采用文件方式,如果单点计算,面板上的数据框也可以进行单点计算,而不用采用文件方式。点击转换按钮,可以同时计算单点和文件。如果没有文件或者没有单点,也没有关系,忽略即可。

    转换后坐标,UTM投影,中央经线120度。
    p1 3318785.352607667900 500000.000000000000 999.999999991618
    p1 3374297.773515739000 547980.561370084410 999.999999989755
    p3 3407710.849772922700 576532.969669877670 999.999999989755

    下面有刚才的两种数据进行七参数计算。
    在这里插入图片描述
    1选取源数据,2选取转换后数据,34选取数据对应的坐标系,5进行计算。计算完成后,如果需要使用该七参数进行坐标转,需要按6启用七参,再通过坐标转换按钮进行计算。
    七参计算后,在程序文件夹下回生成一个report报告:
    在这里插入图片描述
    包含七参数的值,以及三个方向的误差,以及按照七参数将源数据进行转换,和转换后的数据进行比较,很直观看出来点位的精度和偏差。

    注意:
    七参计算和坐标转换共享一套坐标系选项。七参数只是表示连个椭球之间的空间关系,和坐标系,投影方式以及分度带无关。七参计算完毕后,只要椭球没变,可以任意改变投影的方式。
    另外:
    如果你有七参数,可以填到七参数文本框中,点击启用,计算的时候会将七参数加入计算过程。依然可以得到正确的结果。
    由于数据保密原因,随手写了几个坐标进行转换,读者可以使用真实数据进行验证。该软件进行过大量数据验证,请放心使用。

    展开全文
  • 通过建立虚拟相机坐标系以及分析条纹信息在投射器坐标系相机坐标系之间的转换关系, 在相机坐标系中建立了从相位到高度的一一映射模型, 并在映射模型中校正了相机镜头的畸变。该模型结构简单, 便于结合查表(LUT)法...
  • 研究感知亮度光亮度、色度之间的关系,并针对激光投影电视设计了两次异色感知亮度匹配实验。实验一中基于两台具有相同参数的单色激光电视,研究了色饱和度、色调对感知亮度的影响。结果表明,当在同一色调,光亮度相同...
  • 输入和输出的地理坐标系相同:直接转换投影坐标系。 栅格数据:ArcToolbox—Data management tools—Projections andtransfomations— Raster—Project raster 矢量数据:ArcToolbox—Data management tools—...

    前言

    在处理地理数据时,对坐标系的概念一直很模糊,特别是操作arcgis时经常出现坐标错误。搜了好多资料都觉得不清不楚。最后看了一些视频才大概明白怎么操作,特此记录关于坐标系的知识


    提示:推荐两个视频:传送门1传送门2
    再推荐一个官方文档吧:传送门3,除了理解起来困难点,官方文档还是权威的。

    一、基本概念

    ArcGIS中预定义了两套坐标系统,地理坐标系(Geographic coordinate system)和投影坐标系(Projectedcoordinate system)。
    在这里插入图片描述

    1.地理坐标系

    地理坐标系定义:地理坐标系是以椭球体面为参考面,以法线为依据,用经纬度表示地面点在椭球表面的位置的坐标系统。

    • 简单点来说,地理坐标系是用经纬度来表示地球表面物体的位置,如下图所示。
    • 不同的地理坐标系的区别就在于用于拟合地球大地水准面的椭球大小和位置有关。
    • 我国常用的地理坐标系有GCS_Beijing_1954, GCS_Xian_1980,CS_WGS_1984,GCS_CN_2000)

    2.投影坐标系

    投影:将球面坐标转为平面坐标的过程。
    投影坐标系的实质是平面坐标系统,地图单位通常为(米)。投影坐标系=地理坐标系+投影函数算法
    投影坐标系通俗来讲就是把地球椭球展平在地图上画出来,既然已经有地理坐标系能够准确表达物体位置了,为啥还需要投影坐标系呢?主要是制图需要,还可以进行长度和面积的量测等。

    投影坐标系有一些知识需要了解

    1. 我国基本比例尺地形图主要采用高斯-克吕格投影(Gauss-Kruger)(1:100万除外,1:100万采用正轴等角割圆锥投影,又叫兰勃特投影(Lambert Conformal Conic)),根据不同的比例尺采用不同的投影分带。
      根据不同的比例尺采用不同的投影分带
    2. 投影分带:高斯克吕格投影采用分带来减小投影的形变,所以选择投影的时候根据中央经线来选取分带号。这个代号是可以用公式计算的,查表也行。6度分带中中国处于13带到23带中,共12个带之间。在3度分带中,中国处于24带到45带共22带之间。
    3. ArcMap中投影坐标系的表示方法:如下图所示。需要注意的是,有带号的坐标系都是8位,前两位表示的是带号。即有带号的投影坐标系坐标是在中央经度投影坐标系的坐标前面加上带号。

    二、相关操作

    讲了那么多,原理再清楚,要知道怎么用才是关键,而且如果不是专业绘图制图的人员,原理不太清楚也没关系。

    1.查看数据框坐标信息

    在arcmap中加载数据时,如果在加载数据之前没有对数据框的坐标进行设置,第一个加载的数据的坐标信息会作为数据框的坐标。如果后来加进来的数据坐标与数据框坐标不符,则arcmap软件会采用动态投影方法基于默认参数将数据加在一起,会存在一些误差。

    常见问题1:加载数据进来,数据不显示

    检查数据框的坐标系与加进来数据坐标系是否相同,如果不同则进行投影变换(ArcToolbox—Data management tools—Projections andtransfomations— Raster/Feature—Project)更改数据坐标系或者重新定义数据框的坐标系。

    常见问题2:想更改数据坐标显示的格式

    在arcgis软件右下角会显示鼠标位置的坐标信息,有时候我们需要显示经纬度信息,有时候可能需要显示平面坐标信息,该怎么设置呢?如图所示找到数据框属性设置display。
    在这里插入图片描述

    2.查看图层坐标信息

    找到要查看的图层或者在catalog里找到要查看的要素,右键properties,在source字段里查看图层的坐标系信息。
    在这里插入图片描述
    在框A(Extent)中,标明了该图层的坐标范围,这个框里的信息可以检查是不是投影坐标系(如果框内的范围是经纬度信息,而单位确是米m,则说明投影错了,转到第三步–删除坐标信息再进行操作)。
    在框B(Data Source)中,显示数据的投影坐标系和地理坐标系,可以查看是不是自己需要的坐标系,如果不是则需要进行处理。

    3.删除坐标信息

    很多时候会遇到一个问题就是:命名两个图层的数据是同一个区域的,加载进来或者进行投影变换之后,两个图层数据不是在一个区域内了。出现上述问题说明数据坐标定义或者投影定义错了,需要将坐标的信息进行删除。
    删除的步骤如下:
    1)打开catalog,找到对应的图层,右键properties打开属性对话框。会显示选中的图层是什么坐标系。

    2)点击上图红色框框里的图标,选择clear,确定即可。地理坐标和投影坐标就被删除了。

    常见问题3:在catalog里右键属性直接设置坐标可以吗?

    不行,catalog只能查看和删除坐标系,虽然在里面选一个其他的坐标系点击确定之后会显示相关的属性,但是文件会出现缩放问题,这不是正确的定义坐标系方法。

    4.地理坐标系定义与转换

    坐标系定义

    如果某一个数据集的坐标系未知、被删除了或是不正确,可以定义地理坐标系。注意使用这个工具的前提是知道该数据的正确坐标系
    工具位置:ArcToolbox—Data management tools—Projections andtransfomations—Define Projections。
    在这里插入图片描述

    地理坐标系转换

    1. 地理坐标系转换是涉密的,不同椭球之间坐标系转换要先自己定义(ArcToolbox—Data management tools—Projections and transfomations-Create Custom Geographic Transformation)
      在这里插入图片描述
    2. 在方法选择栏里:三参数转换选GEOCENTRIC_TRANSLATION,七参数选POSITION_VECTOR,然后输入相关参数。
    3. 在投影(ArcToolbox—Data management tools—Projections andtransfomations— 栅格或者矢量—Project)里选择自己定义的参数转换方法,便可以完成地理坐标系的转换。

    5.投影坐标系定义与转换

    投影变换可以理解为不同坐标系之间的转换,按照坐标系的不同可以分为:
    输入和输出的地理坐标系不同:先转换地理坐标系(如4的方法),再转换为投影坐标系。
    输入和输出的地理坐标系相同:直接转换投影坐标系。

    • 栅格数据:ArcToolbox—Data management tools—Projections andtransfomations— Raster—Project raster
    • 矢量数据:ArcToolbox—Data management tools—Projections andtransfomations—Feature—Project

    总结

    很多时候,需要自己掌握原理和一些知识,一直问别人是不靠谱的。
    展开全文
  • 栅格数据投影转换

    万次阅读 2018-06-01 23:24:15
    栅格数据投影转换 作者:阿振 邮箱:tanzhenyugis@163.com 博客:https://blog.csdn.net/theonegis/article/details/80089375 修改时间:2018-06-01 声明:本文为博主原创文章,转载请注明原文出处 使用...

    栅格数据投影转换

    作者:阿振

    邮箱:tanzhenyugis@163.com

    博客:https://blog.csdn.net/theonegis/article/details/80089375

    修改时间:2018-06-01

    声明:本文为博主原创文章,转载请注明原文出处


    使用GDAL提供的命令行工具进行转换

    GDAL提供了gdalwarp命令可以方便地让我们进行影像拼接,重投影,裁剪,格式转换等功能

    比如,我们需要将MODIS数据的Sinusoidal投影转为UTM投影,我们可以这样操作。

    我需要转换的地区位于UTM的49度带内,我查看了一下其EPSG的编码为:EPSG:32649(WGS 84 / UTM zone 49N)

    注:推荐大家一个网站,可以查阅各种投影的定义:http://spatialreference.org

    然后,终端中执行如下命令:

    gdalinfo MOD09A1.A2017361.h28v06.006.2018005034659.hdf (用于查看MODIS数据中的波段名称与地址,这里我们只转换第一波段)

    gdalwarp -t_srs EPSG:32649 HDF4_EOS:EOS_GRID:"MOD09A1.A2017361.h28v06.006.2018005034659.hdf":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01 MODSI_WARP_32649.tif-t_srs参数用于指定输出投影信息,可以是EPSG,或者OGC WKT,或者PROJ4格式,后面分别是输入数据和输出数据文件名)

    使用代码进行转换

    使用命令行转换,当然有两种方法啦:

    第一,直接在代码中调用命令行工具的接口(比较懒的人,像我,当然直接用第一种啦,有现成的工具为什么不用);

    第二,自己做投影转换之后的坐标计算,主要是计算重投影之后的GeoTransform参数,有了GeoTransform参数以及投影的定义,我们就可以通过SetGeoTransform()SetProjection()投影转换了.

    下面我给出具体的实现代码:

    第一种方法直接调用gdal.Warp()方法,该方法其实就是对gdalwarp命令的封装,第一个参数是输出文件,第二个参数是输入文件或者输入的Dataset,后面的都是可选参数(dstSRS参数指定输出投影)

    from osgeo import gdal
    
    root_ds = gdal.Open('MOD09A1.A2017361.h28v06.006.2018005034659.hdf')
    # 返回结果是一个list,list中的每个元素是一个tuple,每个tuple中包含了对数据集的路径,元数据等的描述信息
    # tuple中的第一个元素描述的是数据子集的全路径
    ds_list = root_ds.GetSubDatasets()
    
    # 取出第1个数据子集(MODIS反射率产品的第一个波段)进行转换
    # 第一个参数是输出数据,第二个参数是输入数据,后面可以跟多个可选项
    gdal.Warp('reprojection.tif', ds_list[0][0], dstSRS='EPSG:32649')
    
    # 关闭数据集
    root_ds = None

    在介绍第二种方法之前,我们有必要回忆一下之前说过的GDAL反射变换的六参数模型:

    放射变换使用如下的公式表示栅格图上坐标和地理坐标的关系:

    Xgeo=GT(0)+XpixelGT(1)+YlineGT(2)Ygeo=GT(3)+XpixelGT(4)+YlineGT(5) X g e o = G T ( 0 ) + X p i x e l ∗ G T ( 1 ) + Y l i n e ∗ G T ( 2 ) Y g e o = G T ( 3 ) + X p i x e l ∗ G T ( 4 ) + Y l i n e ∗ G T ( 5 )

    Xge0 X g e 0 , Yge0 Y g e 0 )表示对应于图上坐标( Xpixel X p i x e l , Yline Y l i n e )的实际地理坐标。对一个上北下南的图像,GT(2)和GT(4)等于0, GT(1)是像元的宽度, GT(5)是像元的高度的相反数。(GT(0),GT(3))坐标对表示左上角像元的左上角坐标。

    通过这个放射变换,我们可以得到图上所有像元对应的地理坐标。

    好了,所以我们需要计算对于上面的六参数,我们主要需要计算重投影以后图像左上角的坐标(最小的X坐标值和最大的Y坐标值),这个转换我们可以通过osr.CoordinateTransformation类进行,下面给出实现代码:

    from osgeo import gdal
    from osgeo import osr
    
    # root_ds = gdal.Open('/Users/tanzhenyu/Resources/DataWare/MODIS/MOD09A1.A2017361.h28v06.006.2018005034659.hdf')
    # # 返回结果是一个list,list中的每个元素是一个tuple,每个tuple中包含了对数据集的路径,元数据等的描述信息
    # # tuple中的第一个元素描述的是数据子集的全路径
    # ds_list = root_ds.GetSubDatasets()
    #
    # # 取出第1个数据子集(MODIS反射率产品的第一个波段)进行转换
    # # 第一个参数是输出数据,第二个参数是输入数据,后面可以跟多个可选项
    # gdal.Warp('reprojection.tif', ds_list[0][0], dstSRS='EPSG:32649')
    #
    # # 关闭数据集
    # root_ds = None
    
    
    def reproject(src_file, dst_file, p_width, p_height, epsg_to):
        """
        :param src_file: 输入文件
        :param dst_file: 输出文件
        :param p_width: 输出图像像素宽度
        :param p_height: 输出图像像素高度
        :param epsg_to: 输出图像EPSG坐标代码
        :return:
        """
        # 首先,读取输入数据,然后获得输入数据的投影,放射变换参数,以及图像宽高等信息
        src_ds = gdal.Open(src_file)
        src_srs = osr.SpatialReference()
        src_srs.ImportFromWkt(src_ds.GetProjection())
    
        srs_trans = src_ds.GetGeoTransform()
        x_size = src_ds.RasterXSize
        y_size = src_ds.RasterYSize
        d_type = src_ds.GetRasterBand(1).DataType
    
        # 获得输出数据的投影,建立两个投影直接的转换关系
        dst_srs = osr.SpatialReference()
        dst_srs.ImportFromEPSG(epsg_to)
        tx = osr.CoordinateTransformation(src_srs, dst_srs)
    
        # 计算输出图像四个角点的坐标
        (ulx, uly, _) = tx.TransformPoint(srs_trans[0], srs_trans[3])
        (urx, ury, _) = tx.TransformPoint(srs_trans[0] + srs_trans[1] * x_size, srs_trans[3])
        (llx, lly, _) = tx.TransformPoint(srs_trans[0], srs_trans[3] + srs_trans[5] * y_size)
        (lrx, lry, _) = tx.TransformPoint(srs_trans[0] + srs_trans[1] * x_size + srs_trans[2] * y_size,
                                          srs_trans[3] + srs_trans[4] * x_size + srs_trans[5] * y_size)
    
        min_x = min(ulx, urx, llx, lrx)
        max_x = max(ulx, urx, llx, lrx)
        min_y = min(uly, ury, lly, lry)
        max_y = max(uly, ury, lly, lry)
    
        # 创建输出图像,需要计算输出图像的尺寸(重投影以后图像的尺寸会发生变化)
        driver = gdal.GetDriverByName('GTiff')
        dst_ds = driver.Create(dst_file,
                               int((max_x - min_x) / p_width),
                               int((max_y - min_y) / p_height),
                               1, d_type)
        dst_trans = (min_x, p_width, srs_trans[2],
                     max_y, srs_trans[4], -p_height)
    
        # 设置GeoTransform和Projection信息
        dst_ds.SetGeoTransform(dst_trans)
        dst_ds.SetProjection(dst_srs.ExportToWkt())
        # 进行投影转换
        gdal.ReprojectImage(src_ds, dst_ds,
                            src_srs.ExportToWkt(), dst_srs.ExportToWkt(),
                            gdal.GRA_Bilinear)
        dst_ds.GetRasterBand(1).SetNoDataValue(0)  # 设置NoData值
        dst_ds.FlushCache()
    
        del src_ds
        del dst_ds
    
    
    if __name__ == '__main__':
        src_file = 'HDF4_EOS:EOS_GRID:"MOD09A1.A2017361.h28v06.006.2018005034659.hdf":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01'
        dst_file = 'reprojection.tif'
        reproject(src_file, dst_file, 450, 450, 32649)
    展开全文
  • 透视投影的坐标转换与数学推导

    千次阅读 2018-10-30 17:56:00
    坐标转换 管线在处理顶点数据时需要经过多个坐标系的转换。 1.通过modelview矩阵先...3.透射投影将顶点转换到了单位立方体中,变换到标准化设备坐标系。 4.而标准化设备坐标系中的物体常常会产生形变,为恢复...

    坐标转换


    管线在处理顶点数据时需要经过多个坐标系的转换。

    1.通过model与view矩阵先将其变换到世界坐标系中,再将其变换到观察坐标系中。从而方便后续处理。

    2.经透射投影,将3D坐标映射到视平面,完成从3D到2D的变换,从而在后续可转换为在屏幕显示的2D坐标。

    3.透射投影将顶点转换到了单位立方体中,变换到标准化设备坐标系。

    4.而标准化设备坐标系中的物体常常会产生形变,为恢复物体原比例,将其进行视口变换,转换到最后显示的屏幕坐标系中。

     

    以上为顶点在坐标转换中的大致流程,其中包含的很多细节部分会在下面结合透射投影进行详细分析。

     

    逻辑分析


     

                     图 1                                            图 2                                            图 3                                                         图4

     

    1.转换:我们输入在程序中的坐标都是3D,而最终显示在屏幕上的只能是2D。所以可想而知,需要一个过程来将3D坐标映射到我们的2D视平面上。只有这样我们才能通过显示器去观看。

     

    2.观察:如上述图1,我们在观察面上所观察到的正是观察点E与物体实际坐标的连线在视平面P上的交点,其比例也可形成近大远小的现实感。图2可直观的感受。

     

    3.问题:根据学过的知识,显然这一交点可以很轻松的通过相似三角形来求得。但问题是,我们需要获得多大范围内的物体?即限定观察空间。不在观察空间内的物体又该怎样去处理?即裁剪。

     

    4.观察空间:设定近平面,远平面,和观察角度,来得到一个视椎体(frustum)。而我们也常常将近平面作为视平面。如图3。

     

    5.裁剪:裁剪掉视锥体以外的顶点,只保留以内的。若直接对视锥体裁剪,因平面的倾斜也会导致计算量过大,所以采用将其映射到单位立方体中(CVV),此时都为垂直平面,并且在映射过程中,顶点与视锥体的关系不变,在视锥体外的仍然在视锥体外。裁剪时只需判定一个坐标即可,大大减少了计算量。如图4。

     

    6.恢复:但是在映射到CVV的时候,顶点间的比例产生了变化,物体常常产生形变,所以在映射到屏幕坐标系时,必须对其进行比例恢复。采用的方法是:将近平面(即视平面)的宽高比例一开始就调整为与视口宽高比相同。所以在恢复比例时,只需要再与比例进行乘除等操作即可。(视口并非电脑屏幕大小,视口指的是我们要渲染的尺寸。视口是<=电脑屏幕大小的,我们通常在程序中将视口与窗口大小保持一致,所以视口也表现为程序运行后的窗口尺寸。)

     

    数学推导(2种表达式)

     

     

    透视投影的过程分为4部分进行。

    1、将视锥体中的顶点投影到视平面 

    2、将视锥体映射到CVV(规则观察体)中进行裁剪 

    3、进行透视除法。

    4、进行视口变换,恢复其比例,最终将其映射到屏幕坐标系中。

     

     

    公式推导1(线性插值法)

     

    参数介绍:np:近截面  fp:远截面   p':投影后的顶点  p:物体顶点  N:近截面到观察点的距离   F: eye:观察点   (取近截面为投影面

     

    根据p与p'的相似三角形性质可得:

    同理,在z轴和y轴上可得:

    从而得到投影后得坐标:

     

    到此,已把顶点投影到视平面,此时所有点的z值均为-N。若以此为最终坐标,则无法在后续中进行深度测试等操作,所以需要保留原本顶点在空间中的z值。

    由此可把坐标最终定为:

    此坐标的x与y是相对于视平面的,而z是相对于原先的空间坐标。可以想象若用此坐标来构成原物体,则已经产生形变。

     

     

    在逻辑分析中已解释了在CVV中裁剪的原因,接下来就是把视锥体映射到CVV中,即把视锥体内的物体都变换到CVV中。

    1. 涉及到变换,则需要构造出一个变换矩阵,从而用此矩阵来完成映射。而矩阵中的各种变换就要求我们把顶点坐标转换为齐次坐标。

    2.被忽视的重要环节:任何顶点最终显示在屏幕上,都要经过插值和光栅化的过程,然后再渲染光栅化得到的片段,最终显示在屏幕的像素中。而透视投影的插值也是使用线性插值,但是x',y'与z并非线性关系,经过推导可得出x',y'与1/z是成线性关系的。若能把p'中的z变换为带有1/z的形式,则会更加方便后期操作。

    3.由于要把p'转换为齐次坐标,而齐次坐标在变换为普通坐标时,需要同除w,x'和y'中均为x和y同除以了-z。所以也需要把p'中的z变换为带有-1/z的形式,从而更方便齐次坐标到普通坐标的转换。

    4.要把z方向上的值都映射到CVV中,即区间 [-1,1]。可以选取合适的系数,当z=-N时,取-1。当z=-F时,取1。

     

    综上多点,我们对p'的z进行形式转化,最终得到:

    对于这种形式,就可以很轻松的进行矩阵和齐次变换了

    其中

     

    先在z轴上进行从视锥体到CVV的映射

     

    再对x轴和y轴进行映射(此时是把投影后的x'和y'坐标,映射到CVV中)

     

    我们知道-Nx / z的有效范围是投影平面的左边界值(left)和右边界值(right),即区间[left, right],而y的边界值-Ny / z则映射到区间[bottom, top]。而现在我们想把-Nx / z属于[left, right]映射到-Nx / z属于[-1, 1]中,-Ny / z属于[bottom, top]映射到-Ny / z属于[-1, 1]中。采用的方法是线性插值

     

     

    综上可得p'的坐标为:

    这个坐标是齐次坐标经透视除法后得到的,我们做透视除法的逆操作,同乘1/w得:

    变换为矩阵形式:

    从而求得我们想要得最终变换矩阵M表达式:

     以上就是透视投影矩阵的数学推导过程。但是到此仍然有很多问题,因为只知道公式是无法具体结合应用的,而且在很多程序中给出的都是类似于aspect的公式,好像并没有用到上述推导的公式。还有上述的left和right边界在视口和屏幕坐标系中又是何种意义?在以下,我将进行说明。

     

     

    公式推导2(等比变换法)

     

    在程序种,我们很少用到以上形式,因为参数过多,较为复杂,更为普遍的是以下形式。2种形式的区别在于,对x和y进行映射时采用的方法不同

    公式1采用线性插值,公式2采用等比的方式。

    Z方向的推导过程不变,最终得到的仍为:

    由于此时,x'=-Nx/z与y‘=-Ny/z均在视平面上,只需要将视平面变换到CVV的单位截面即可。设最终变换到CVV中的坐标为x'' , y'',投影平面高为H。因为CVV长宽高均为2

    则由等比变换可得:y'/y''= H/2

    解得:y''=-(y/z)*cot(fovy/2)

    同理可得:x''=-(x/z)*cot(fovy/2)*(1/aspect)

    求矩阵时对其做齐次坐标逆处理,同乘-z。从而综上可得矩阵M为

     

    对比2种推导过程,公式1采用线性插值更复杂,所以常用的是公式2中的等比变换。得到的矩阵也更为简单。所以在程序中若打开底层矩阵的定义,我们可发现采用的均为公式2。

     

     

    参数说明

     

     

    一、投影平面的上下左右边界与代码的关系。

    1.在正常的程序中用到的代码为 glm::perspective( glm::radians(fov) , (float)SCR_WIDTH / (float)SCR_HEIGHT ,  0.1f , 100.0f ) ;以此构造出透视投影矩阵,第一个参数为fov,即y方向上视野角度的大小。第二个参数为我们设定运行后窗口的宽高比(即aspect)。第三、四个分别为近截面和远截面距观察点的距离。

    2.令fovy=radians(fov) 。则tan(fovy/2)=height/0.1f,从而可求出height。再由height乘上aspect可得出width。从而确定了视锥体前后上下左右6个面的坐标。由此就可求出之前推导出的矩阵M中的所有参数了。

    3.直接看glm::perspective的参数,会产生迷惑,因为这个函数隐藏了利用fov和aspect求各个面的过程。

    由此图也可看出,函数中另有一个函数glFrustum。

     

    二、投影平面的边界与管线流程的关系

    1.因为我们通常将视口渲染尺寸设为窗口尺寸,所以以上所有推导的过程中都包含着视口尺寸=窗口尺寸的隐藏条件。

    2.在 glm::perspective( glm::radians(fov) , (float)SCR_WIDTH / (float)SCR_HEIGHT ,  0.1f , 100.0f ) 中,把投影面的宽高比设定成了和视口宽高比相同。这是因为在讲投影面上的图形映射到CVV中时,CVV的宽高比为1,这就导致图形产生失真现象。解决这一问题,就是在视口变换中按照投影面原先设定的比例变换到视口中,又因为原先设定的比例就与视口相同,所以最终渲染出的比例当然也就是恢复原样了。而视口尺寸比例又与窗口尺寸比例相同,则最终在窗口中看到的图形为正常投影后的保真图形。

     

    转载于:https://www.cnblogs.com/jingrui/p/9874194.html

    展开全文
  • (个人也在学习中) 1.椭球体、基准面及地图投影GIS中的坐标系定义是GIS系统的...GIS中的坐标系定义由基准面和地图投影两组参数确定,而基准面的定义则由特定椭球体及其对应的转换参数确定,因此欲正确定义GIS系统...
  • 地理坐标系和投影坐标系之间的关系 基本概念 地理坐标系:为球面坐标。 参考平面地是椭球面,坐标单位:经纬度; 投影坐标系:为平面坐标。参考平面地是水平面,坐标单位:米、千米等; 地理坐标转换投影坐标...
  • 立方体投影算法

    2019-04-13 20:37:49
    投影,数学含义是两个面上的点点的对应关系;在一个面上的点,另一个面上只有唯一的点之对应。每种地图投影都有一组必须定义的参数参数用于指定原点以及为感兴趣区域自定义投影。角度参数使用地理坐标系单位,...
  • 地理坐标系和投影坐标系之间的关系

    万次阅读 多人点赞 2017-12-01 11:26:42
    参考平面地是椭球面,坐标单位:经纬度;...地理坐标转换投影坐标的过程可理解为投影。(投影:将不规则的地球曲面转换为平面)地理坐标系地球的三级逼近大地水准面地球的自然表面有高山也有洼地,
  • 点的中心投影:点的中心投影一般仍是一个点,但当投影线像面平行时,投影点将位于无穷远处。 线段的中心投影:线段的中心投影一般仍是一条线段,但也可能为一个点或者一条射线。 相交线段的中心投影:一
  • 在GIS数据处理过程中,投影是其中一个重要的环节,很多人对定义投影投影转换容易混淆,下面对定义投影投影转换进行介绍: 定义投影简单来说就是对数据定义投影信息,一般也可以分为三种情况: 一是给没有...
  • 本文为基于arcgis转换经纬度和投影坐标教程以及关于它们之间关系的一点理解
  • (一)人脸识别技术之人脸识别过程及识别算法简介

    万次阅读 多人点赞 2018-11-04 23:19:40
    人脸特征提取是将人脸图像信息数字化,将一张人脸图像转换为一串数字(一般称为特征向量)。如对一个人脸,找到他的眼睛左边、嘴唇右边、鼻子、下巴等位置,利用特征点间的欧氏距离、曲率和角度等提取出特征分量,最终...
  • 所谓动态投影指,ArcMap中的Data 的空间参考或是说坐标系统是默认为第一加载到当前工作区的那个文件的坐标系统,后加入的数据,如果和当前工作区坐标系统不相同,则ArcMap会自动做投影变换,把后加入的数据投影变换...
  • matlab人脸识别论文

    万次阅读 多人点赞 2019-10-11 17:41:51
    摘 要 本文设计了一种基于BP神经网络的人脸识别系统,并对其进行了性能分析。该系统首先利用离散小波变换获取包含人脸图像大部分原始信息的低频分量,对图像数据进行降维;...通过系统仿真实验分析发现:人脸特征的提...
  • 【数据库学习】数据库总结

    万次阅读 多人点赞 2018-07-26 13:26:41
    规范化的方法:将关系模式投影分解成两个或两个以上的关系模式。 2,依赖和范式 1)依赖 ①部分函数依赖 设X,Y是关系R的两个属性集合,存在X→Y,若X’是X的真子集,存在X’→Y,则称Y部分函数依赖于X。 举个例子:...
  • 前端面试题

    万次阅读 多人点赞 2019-08-08 11:49:01
    前端面试题汇总 ... 你做的页面在哪些流览器测试过?这些浏览器的内核分别是什么? 21 ... 21 Quirks模式是什么?它和Standards模式有什么区别 21 ...img的alttitle有何异同? strongem的异同? 22 你能...
  • 目录 1.Intro ...(2)高斯投影正反算 4.Environment and Configuration 5.Conclusion 1.Intro 搞了一段时间的WebGIS和接口,突然被领导喊来研究松散型切片(PNG)的切片方法,看到某人写的源码...
  • TensorFlow入门

    千次阅读 多人点赞 2019-04-23 10:09:29
    在实现上, TensorFlow 将图形定义转换成分布式执行的操作, 以充分利用可用的计算资源(如 CPU或 GPU). 一般你不需要显式指定使用 CPU 还是 GPU, TensorFlow 能自动检测. 如果检测到 GPU, TensorFlow 会尽可能地利用...
  • 投影分度带的划分 第一步:查看你所下载(或者要套合的范围)的图像的经纬度中【经度】所在的范围,如下图: 比如,你要图中这个区域的地图来套合你的矢量数据,点击红色箭头指向的地方,如下图: 从上图中我们...
  • OpenCV_估算图像之间的投影关系

    千次阅读 2017-12-18 10:56:16
    相机校准就是设置相机各种参数的过程,就是用相机拍摄特定的图案并分析得到的图像,然后再优化的过程中确定最佳的参数值。 OpenCV推荐使用国际象棋棋盘的图案生成用于校准的三维场景的集合,并且由于图案视平面的,...
  • 通过mapgis转换经纬度与投影坐标,arcgis里的教程参加我的其他博客,可以在分类ArcGIS里查到,详细介绍了投影坐标系经纬度(地理坐标系)的关系。。而且arcgis同样可以进行坐标转换
  • 经纬度对象的投影转换

    千次阅读 2019-07-17 11:15:03
    EPSG:4326(又名WGS84,未投影)是一个地理的非项目坐标系。它是lat,longs GPS显示器。它的单位是十进制度。EPSG:4326无法在平面地图上以有意义的方式显示。 EPSG:3857(又名Pseudo-Mercator,球形墨...
  • 正交投影变换透视投影

    万次阅读 2015-07-15 14:29:59
    相机投影模型 三维计算机图形学的基本问题之一就是三维观察问题:即如何把三维场景投影到要显示的二维图像。大多数经典的解决投影变换方法有两种:正交投影变换和透视投影变化。  正交投影变换用一个长方体来...
  • 地理坐标系与投影坐标系

    千次阅读 2018-08-13 16:38:58
    地理坐标系与投影坐标系 1.基本概念 地理坐标系:为球面坐标。 参考平面地是椭球面,坐标单位:经纬度; 投影坐标系:为平面坐标。参考平面地是水平面,坐标单位:米、千米等; 地理坐标转换投影坐标的过程可...
  • 地理坐标系与投影坐标系辨析

    千次阅读 2019-06-16 16:39:48
    特别地,等距圆柱投影以图像像素点地理位置简单的对应关系已经成为了全球栅格数据集的标准,比如Celestia (一套开放原始码的 3D 天文软件)和NASA的World Wind。 WOW 很小的事情中竟包含如此多的的知识 ...
  • 基本概念 首先简单介绍一下地理坐标系、大地坐标系以及地图投影的概念: 地理坐标系:为球面坐标。 参考平面地是椭球面,... 地理坐标转换投影坐标的过程可理解为投影。(投影:将不规则的地球曲面转换为平面)...
  • 对于坐标系和坐标系之间转换的记录是油生已久的想法,总的来说,坐标系这个概念从毕业以后,参与的两份工作中都离不开它。 由于相关内容较多,重点按照链接在后面介绍中一一做记录。 坐标系是用于表示地理要素、...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 21,949
精华内容 8,779
关键字:

参数转换与投影的关系