目录html
今天(2019年5月30日)去编译最新版本的GDAL,发现其对Proj.4的依赖已经要求为6.x版本了。因而去https://github.com/OSGeo/proj.4看了一下最新的代码,又去https://proj4.org/看了一下文档,感受5.x和6.x的更新挺大的,有必要测试一下,看工做中的项目是否是要升级过来。git
我没有仔细去看5.x版本的代码,仅看了一下最新的Proj.4
版本6的代码,与早前使用的4.9.3版本简单对比了一下,感受区别仍是挺大的,这里列出几点我关注的地方的对比。github
一、新版本改用C++编写,相比4.9版本代码量增长了很多,功能也多了很多。代码层次结构清晰了许多,好比各类转换算法都在src/transformations目录下能够找到,各类投影方法相关的算法都在src/projections目录能够找到。算法
二、支持了从WKT/WKT2
字符串和EPSG代码直接建立坐标系对象,也支持导出WKT字符串。老版本中记录EPSG坐标系定义的的nad/epsg被弃用,改用SQLite数据库来记录(在data/sql目录下保存着用于生成proj.db
文件的SQL脚本),不过新版本须要依赖SQLite3。sql
三、新版的实现使用了缓存机制,在建立操做坐标系对象及搜索查找等都有用到。代码可见 src/iso19111/factory.cpp、src/iso19111/crs.cpp、src/iso19111/coordinateoperation.cpp、数据库
src/iso19111/coordinateoperation.cpp 等文件。编程
四、新版添加了proj_math.h、math.cpp,添加了pj_hypot
等函数,这解决了一些编译问题(由于以前版本projects.h中声明了hypot
函数,但这个函数在非_WIN32环境中也多是存在math.h中的)。api
如下主要翻译自:PROJ.4 News缓存
此版本的 PROJ 对系统的大地测量功能 (主要是) 引入了一些重要的扩展和改进。架构
引入新功能的主要驱动因素是动态参考框架的出现、高精度全球导航卫星系统的使用日益增长以及对精确坐标变换的相关需求的增长。虽然旧版本的 PROJ 包含一些大地测量功能, 但新框架为将 PROJ 转变为通用地理空间坐标转换引擎奠基了基础。
内部架构也有了许多变化和改进。到目前为止,这些改进都遵循现有的编程接口。可是这个过程已经显示出须要简化和减小代码库,以支持持续的主动开发。
新的主要版本号使该项目在名称上留下了一些难题。在产品的大部分使用寿命中,它被称为PROJ.4,但因为咱们如今已达到版本5,所以名称再也不与版本号对齐。
所以,咱们决定将名称与版本号和该版本分离,而后将产品简称为PROJ。为了表彰软件的历史,咱们将PROJ.4做为组织项目的名称。同一个项目团队也会生成datum-grid 包。
综上所述:
推出新的API在proj.h
proj_
命名空间(名称前缀)PJ_
命名空间(名称前缀)引入“转换管道”(transformation pipelines)的概念,能够经过 菊花链 的方式简化坐标操做,能够对坐标进行复杂的大地转换。
采用 OGC/ISO-19100 地理空间标准系列术语。关键定义是:
新操做
pipeline
)helmert
)horner
)hgridshift
)vgridshift
)molodensky
)deformation
)unitconvert
)axisswap
)ccon
)新的“自由格式”(free format)选项,运行在指定key/value键值对时候经过空白字符进行分割标记,例如proj = merc lat_0 = 45
添加到init-files的元数据,能够经过proj.h
中新的API函数proj_init_info()
读取它们。
添加了具备ITRF转换参数的ITRF2000,ITRF2008和ITRF2014初始化文件,包括板块运动模型参数。
添加椭球参数到GSK2011,PZ90和"danish"。后者相似于已经支持的andrae椭球体,但长半轴略微不一样。
添加哥本哈根本初子午线.
将EPSG数据库更新至9.2.0版。
Geodesic库已更新至1.49.2-c版。
对分析性偏导数的支持已被移除。
改善了Winkel Tripel和Aitoff的表现。
将pj_has_inverse()
函数引入proj_api.h
,使用它检测一个操做是否可反转,而不是检测P->inv
是否存在。
ABI版本号更新为13:0:0。
删除了对Windows CE的支持。
删除了VB6 COM接口。
proj_errno_string()
到 proj.h
--skip-lines
选项。NaN
处返回Nan
值。src/org_proj4_Projections.h
文件。+t_epoch
为时间输入。-multiplier
选项到vgridshift。+transpose
选项替换为+convention
。从如今开始应当显示指定convertion使用,使用+transpose选项将返回错误。PROJ6 进行了普遍的更改, 以增长其功能范围, 从具备所谓 "早期绑定" 大地测量基准转换功能的制图投影引擎, 到更完整的库, 支持坐标变换和坐标参考系统。
做为其余加强功能的基础, PROJ 如今包括由 iso-19111:2019 标准/OGC 抽象规范主题 2: "按坐标引用" 的模型的 C++ 实现, 用于大地测量参照框 (基准), 坐标参考系统和协调操做。这些大地测量对象的构造和查询可经过新的 C++ API 进行, 而且在很大程度上能够从 C API 中的绑定中访问。
这些大地测量对象能够从 OGC 已知文本格式 (WKT) 以不一样的变体导入和导出: ESRI WKT、GDAL WKT 一、WKT2:2015 (ISO 191.2: 2015) 和 WKT2:2018 (ISO 19162: 2018)。还支持从 PROJ 字符串导入和导出 crs 对象。此功能之前在 GDAL 软件库中可用 (WKT2 支持除外, 这是一项新功能), 如今是 PROJ 不可或缺的一部分。
如今, sqlite3 数据库文件 proj.db 中提供了一个统一的大地测量对象数据库、坐标参考系统及其元数据以及这些 CRS 之间的坐标操做。这包括从 IOMP EPSG 数据集 (v9.6.0 版本)、法国国家测绘机构大地测量登记册和 ESRI 投影引擎数据库中导入的定义。PROJ 如今是此 CRS 和坐标操做数据库的 "OSGeo C stack" 中的参考软件, 而之前此功能分布在 PROJ、GDAL 和歌词地理, 并使用 CSV 或其余基于特定文本的格式。
添加了考虑到元数据 (如使用区域和准确性) 的后期绑定坐标操做功能。这能够在许多状况下避免过去使用 WGS84 做为中间系统的要求, 这可能会致使没必要要的精度损失, 或者在没法转换到 WGS84 的状况下有时是没法实现的。这些后期绑定功能如今由 proj_create_crs_to_crs() 函数和 cs2cs 实用程序使用。
添加了一个新的命令行实用程序 projinfo 来查询有关数据库的大地测量对象的信息, 从 WKT 字符串和 PROJ 字符串导入和导出大地测量对象, 并显示两个 CRS(坐标参考系统) 之间可用的坐标操做。
ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
宏接口才可用。 (#836)proj_def.dat
默认文件的支持。 (#201)proj.h
。(#1175)nad2bin
程序。如今能够在proj-datumgrid仓库中查看。proj_geocentric_latitude()
APIetmerc
)。旧的实现能够经过添加 +approx
参数(#404)PROJ_LIB
环境变量中使用多个目录 (#1218)+t_obs
参数,从helmert 和deformation中 (#1264)+dt
参数到deformation中,用于替换移除的 +t_obs
(#1264)-k ellipsoid
选项到 projinfo (#1338)proj_normalize_for_visualization()
尝试应用大多数 GIS 应用和 PROJ <6 使用的轴序 (#1387)PROJ_LIB
搜索路径 (#1398)如下内容主要来自Version 4 to 5/6 API Migration。
这是但愿将代码迁移到使用PROJ 5的开发人员的过渡指南。
原文太长,这里简单抽取一部分。
一、以前老的API两个坐标参考系统之间的任何转换都必须经过未明肯定义的WGS84框架进行中转,新的API取消了对PROJ中转换的限制。虽然任然能够进行这种类型的转换,可是在多数状况下,有更好的替代方案。
二、若是你只关心到米级精度,那么旧的API是够用的。可是WGS84并不是真实世界的基础,其它一切均可以经过WGS84进行转换的观点是有缺陷的。而且这里说的WGS84是6个实现中的哪个呢?
三、对许多坐标系统而言,转换到WGS84,在旧系统之间可能变换精度在厘米级。
四、hub(这里说的就是把全部坐标之间的转换都经过WGS84作中转)框架("基准")概念自己没有大问题。但世界本质是4维的,为了得到大地测量精确变换,可能须要在四个维度下计算,新的API容许这样作。
旧的API下,坐标从坐标系A转换到坐标系B的过程,须要从A变换到WGS84(反算),再从WGS84变换到B(正算)
$ echo 300000 6100000 | cs2cs +proj=utm +zone=33 +ellps=GRS80 +to +proj=utm +zone=32 +ellps=GRS80 683687.87 6099299.66 0.00
新的命令行工具cct
使用新的API,因此一样的转换可能结果不同
$ echo 300000 6100000 0 0 | cct +proj=pipeline +step +inv +proj=utm +zone=33 +ellps=GRS80 +step +proj=utm +zone=32 +ellps=GRS80 683687.8667 6099299.6624 0.0000 0.0000
这里显示了旧API和新API之间的区别,并举了几个例子。 下面咱们用两个不一样的API实现相同的程序。 程序从命令行读取输入经度和纬度,并使用墨卡托投影将它们转换为投影坐标。
咱们首先编写PROJ v.4的程序:
#include <proj_api.h> main(int argc, char **argv) { projPJ pj_merc, pj_longlat; double x, y; // 建立坐标参考系对象 if (!(pj_longlat = pj_init_plus("+proj=longlat +ellps=clrk66")) ) return 1; if (!(pj_merc = pj_init_plus("+proj=merc +ellps=clrk66 +lat_ts=33")) ) return 1; // PROJ.4 API 默认的经纬度都是使用弧度值 while (scanf("%lf %lf", &x, &y) == 2) { x *= DEG_TO_RAD; /* 经度 */ y *= DEG_TO_RAD; /* 纬度 */ // 进行坐标转换 p = pj_transform(pj_longlat, pj_merc, 1, 1, &x, &y, NULL ); printf("%.2f\t%.2f\n", x, y); } pj_free(pj_longlat); pj_free(pj_merc); return 0; }
使用PROJ v.5实现的相同程序:
// 使用新的头文件 #include <proj.h> main(int argc, char **argv) { PJ *P; PJ_COORD c; // 建立墨卡托投影坐标系 P = proj_create(PJ_DEFAULT_CTX, "+proj=merc +ellps=clrk66 +lat_ts=33"); if (P==0) return 1; while (scanf("%lf %lf", &c.lp.lam, &c.lp.phi) == 2) { // 度转弧度 c.lp.lam = proj_torad(c.lp.lam); c.lp.phi = proj_torad(c.lp.phi); // 进行坐标转换(正算) c = proj_trans(P, PJ_FWD, c); printf("%.2f\t%.2f\n", c.xy.x, c.xy.y); } proj_destroy(P); }
旧版 API 函数 | 新版 API 函数 |
---|---|
pj_fwd | proj_trans() |
pj_inv | proj_trans() |
pj_fwd3 | proj_trans() |
pj_inv3 | proj_trans() |
pj_transform | proj_trans_array or proj_trans_generic |
pj_init | proj_create() |
pj_init_plus | proj_create() |
pj_free | proj_destroy() |
pj_is_latlong | proj_angular_output() |
pj_is_geocent | proj_angular_output() |
pj_get_def | proj_pj_info() |
pj_latlong_from_proj | No equivalent |
pj_set_finder | No equivalent |
pj_set_searchpath | No equivalent |
pj_deallocate_grids | No equivalent |
pj_strerrno | No equivalent |
pj_get_errno_ref | proj_errno() |
pj_get_release | proj_info() |
这是但愿迁移代码以使用PROJ 6的开发人员的过渡指南。
这里显示了旧API和新API之间的区别,并举了几个例子。 下面咱们用两个不一样的API实现相同的程序。 程序从命令行读取输入经度和纬度,并使用墨卡托投影将它们转换为投影坐标。
咱们首先编写PROJ v.4的程序:
#include <proj_api.h> main(int argc, char **argv) { projPJ pj_merc, pj_longlat; double x, y; // 建立坐标参考系对象 if (!(pj_longlat = pj_init_plus("+proj=longlat +ellps=clrk66")) ) return 1; if (!(pj_merc = pj_init_plus("+proj=merc +ellps=clrk66 +lat_ts=33")) ) return 1; // PROJ.4 API 默认的经纬度都是使用弧度值 while (scanf("%lf %lf", &x, &y) == 2) { x *= DEG_TO_RAD; /* 经度 */ y *= DEG_TO_RAD; /* 纬度 */ // 进行坐标转换 p = pj_transform(pj_longlat, pj_merc, 1, 1, &x, &y, NULL ); printf("%.2f\t%.2f\n", x, y); } pj_free(pj_longlat); pj_free(pj_merc); return 0; }
使用PROJ v.6实现的相同程序:
#include <proj.h> main(int argc, char **argv) { PJ *P; PJ_COORD c; /* 注意: 在PROJ 6中强烈建议不要使用 PROJ 格式字符串来定义 CRS ,由于 PROJ 格式字符串*/ /* 不是描述 CRS 的一种好方式,准确说来是其对大地基准的描述不够详细。 */ /* 使用权威机构提供的坐标系代码(例如"EPSG:4326", etc...)或者WKT字符串来建立,将可以 */ /* 充分利用PROJ的“transformation engine”来肯定两个CRS直接的最佳转换方式 */ P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "+proj=longlat +ellps=clrs66", "+proj=merc +ellps=clrk66 +lat_ts=33", NULL); if (P==0) return 1; { /* 对于特定的使用状况下(转换先后坐标系已知),这是没有必要的。 */ /* proj_normalize_for_visualization() 确保预期坐标顺序和由 proj_trans() */ /* 返回的顺序是 大地坐标系下先经度后纬度,投影坐标系下先东向后北向 */ /* 若是不是使用上面的PROJ字符串,而是使用 "EPSG:XXXX" 代码,这多是必要的。 */ PJ* P_for_GIS = proj_normalize_for_visualization(C, P); if( 0 == P_for_GIS ) { proj_destroy(P); return 1; } proj_destroy(P); P = P_for_GIS; } while (scanf("%lf %lf", &c.lp.lam, &c.lp.phi) == 2) { /* 不须要转换到弧度 */ c = proj_trans(P, PJ_FWD, c); printf("%.2f\t%.2f\n", c.xy.x, c.xy.y); } proj_destroy(P); }
旧版 API 函数 | 新版 API 函数 |
---|---|
pj_fwd | proj_trans() |
pj_inv | proj_trans() |
pj_fwd3 | proj_trans() |
pj_inv3 | proj_trans() |
pj_transform | proj_create_crs_to_crs() + (proj_normalize_for_visualization() +) proj_trans() , proj_trans_array() or proj_trans_generic() |
pj_init | proj_create() / proj_create_crs_to_crs() |
pj_init | proj_create() / proj_create_crs_to_crs() |
pj_free | proj_destroy() |
pj_is_latlong | proj_get_type() |
pj_is_geocent | proj_get_type() |
pj_get_def | proj_pj_info() |
pj_latlong_from_proj | No direct equivalent, but can be accomplished by chaining proj_create() , proj_crs_get_horizontal_datum() and proj_create_geographic_crs_from_datum() |
pj_set_finder | proj_context_set_file_finder() |
pj_set_searchpath | proj_context_set_search_paths() |
pj_deallocate_grids | No equivalent |
pj_strerrno | No equivalent |
pj_get_errno_ref | proj_errno() |
pj_get_release | proj_info() |