Proj.4 升级新版本5.x和6.x

Proj.4 升级新版本5.x和6.x

0、缘起

今天(2019年5月30日)去编译最新版本的GDAL,发现其对Proj.4的依赖已经要求为6.x版本了。因而去https://github.com/OSGeo/proj.4看了一下最新的代码,又去https://proj4.org/看了一下文档,感受5.x和6.x的更新挺大的,有必要测试一下,看工做中的项目是否是要升级过来。git

一、5.x和6.x更新状况简述

我没有仔细去看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.cppsrc/iso19111/crs.cppsrc/iso19111/coordinateoperation.cpp数据库

src/iso19111/coordinateoperation.cpp 等文件。编程

四、新版添加了proj_math.hmath.cpp,添加了pj_hypot等函数,这解决了一些编译问题(由于以前版本projects.h中声明了hypot函数,但这个函数在非_WIN32环境中也多是存在math.h中的)。api

如下主要翻译自:PROJ.4 News缓存

PROJ 5.x 更新

此版本的 PROJ 对系统的大地测量功能 (主要是) 引入了一些重要的扩展和改进。架构

引入新功能的主要驱动因素是动态参考框架的出现、高精度全球导航卫星系统的使用日益增长以及对精确坐标变换的相关需求的增长。虽然旧版本的 PROJ 包含一些大地测量功能, 但新框架为将 PROJ 转变为通用地理空间坐标转换引擎奠基了基础。

内部架构也有了许多变化和改进。到目前为止,这些改进都遵循现有的编程接口。可是这个过程已经显示出须要简化和减小代码库,以支持持续的主动开发。

新的主要版本号使该项目在名称上留下了一些难题。在产品的大部分使用寿命中,它被称为PROJ.4,但因为咱们如今已达到版本5,所以名称再也不与版本号对齐。

所以,咱们决定将名称与版本号和该版本分离,而后将产品简称为PROJ。为了表彰软件的历史,咱们将PROJ.4做为组织项目的名称。同一个项目团队也会生成datum-grid 包。

综上所述:

  • PROJ.4项目提供产品PROJ,如今版本为5.0.0。
  • PROJ的基础组件是库libproj。
  • 其余PROJ组件包括应用程序proj,它为libproj提供命令行界面。
  • PROJ.4项目还分发了基准网格(datum-grid)包,在编写本文时,它是1.6.0版本。

5.0.0 更新

  • 推出新的API在proj.h

    • 新版API增长了4D空间坐标转换功能
    • 新API中的函数使用proj_命名空间(名称前缀)
    • 新API中的数据类型使用PJ_命名空间(名称前缀)
  • 引入“转换管道”(transformation pipelines)的概念,能够经过 菊花链 的方式简化坐标操做,能够对坐标进行复杂的大地转换。

  • 采用 OGC/ISO-19100 地理空间标准系列术语。关键定义是:

    • 在通用层面上,坐标操做是基于从一个坐标参考系统到另外一个坐标参考系统的一对一关系的坐标变换。
    • 变换(transformation )是一种坐标操做,其中两个坐标参考系统基于不一样的基准,例如,从全局参考框架改变到区域框架。
    • 转换(conversion )是一种坐标操做,其中两个坐标参考系统都基于相同的数据,例如,坐标单位的变化。
    • 投影是从椭球坐标系到平面坐标系的坐标转换。虽然投影只是根据标准进行的转换, 但它们在 PROJ 中被视为单独的实体, 由于它们占库中绝大多数操做。
  • 新操做

  • 新的“自由格式”(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接口。

5.1.0 更新

  • 添加函数proj_errno_string()proj.h
  • 验证管道步骤之间的单位,并确保转换的完整性。
  • 在没有参数的状况下调用cctgie程序时打印帮助。
  • CITATION文件添加到源码发行。
  • 添加了Web墨卡托操做。
  • 提升了赤道附近(4326)向spherical Mercator(3857)前向转换的数值精度。
  • 为cct添加了--skip-lines选项。
  • 始终对输入NaN 处返回Nan值。
  • 移除没有使用的src/org_proj4_Projections.h文件。
  • Java Native Interface绑定已更新
  • 水平和垂直网格移位操做扩展到时域(temporal domain)。

5.2.0 更新

  • 在unitconvert中增长了对deg,rad和grad的支持。
  • 当没有另外指定时,假设+t_epoch为时间输入。
  • 添加了逆拉格朗日投影。
  • 添加-multiplier选项到vgridshift。
  • 添加等地球投影。(等地球地图投影是用于世界地图的新的等面积假圆柱投影)
  • 为gie添加了“require_grid”选项。
  • 将 Helmert 变换的+transpose选项替换为+convention。从如今开始应当显示指定convertion使用,使用+transpose选项将返回错误。
  • 改进的逆spherical Mercator投影数值精度
  • 当前cct会将前向输入的坐标,转发到输出流。 (#1111)

PROJ 6.x 更新

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(坐标参考系统) 之间可用的坐标操做。

6.0.0 更新

  • 公开接口中移除projects.h(当前仅做为内部接口)(#835)
  • 不推荐使用proj_api.h接口。头文件仍然可用,但将在下一个主要版本的PROJ中删除。如今须要定义ACCEPT_USE_OF_DEPRECATED_PROJ_API_H宏接口才可用。 (#836)
  • 删除了对nmake构建系统的支持。(#838)
  • 删除对proj_def.dat 默认文件的支持。 (#201)
  • 构建PROJ须要支持C++11编译器。
  • 构建添加对SQLite 3.7的依赖。(#1175)
  • 添加 projinfo 命令行程序。
  • 添加一些用于处理ISO19111功能函数的到proj.h。(#1175)
  • 更新了cs2cs以使用后期绑定特性。(#1182)
  • 移除 nad2bin 程序。如今能够在proj-datumgrid仓库中查看。
  • 在proj中删除了对Chebyshev多项式的支持
  • 从proj.h中移除proj_geocentric_latitude() API
  • 更改proj行为:如今只容许投影的初始化 (#1162)
  • 更改tmerc行为:如今默认扩展横轴墨卡托算法 (etmerc)。旧的实现能够经过添加 +approx参数(#404)
  • 更改行为:默认椭球如今设置为GRS80 (以前是WGS84) (#1210)
  • 容许在 PROJ_LIB 环境变量中使用多个目录 (#1218)
  • 添加Lambert Conic Conformal (2SP Michigan) 投影 (#1142)
  • 添加 Bertin1953 投影 (#1133)
  • 添加Tobler-Mercator 投影 (#1153)
  • 添加Molodensky-Badekas 变换 (#1160)
  • 添加pushpop 坐标操做 (#1250)
  • 添加+t_obs 参数,从helmert 和deformation中 (#1264)
  • 添加+dt 参数到deformation中,用于替换移除的 +t_obs (#1264)

6.1.0 更新

  • 包含 QGIS 中定义的椭球 (#1137)
  • 添加 -k ellipsoid 选项到 projinfo (#1338)
  • cs2cs支持4D坐标 (#1355)
  • WKT2 解析:更新到OGC 18-010r6 (#1360 #1366))
  • 更新googletest内部版本到v1.8.1 (#1361)
  • 数据库更新:EPSG v9.6.2 (#1462), IGNF v3.0.3, ESRI 10.7.0 并添加operation_version列(#1368)
  • 添加proj_normalize_for_visualization() 尝试应用大多数 GIS 应用和 PROJ <6 使用的轴序 (#1387)
  • 添加 noop(空)操做 (#1391)
  • 用户设置的路径优先于 PROJ_LIB 搜索路径 (#1398)
  • 缩小数据库大小 (#1438)
  • 添加对compoundCRS 和concatenatedOperation组件命名的支持 (#1441)

二、从PROJ.4向新版本迁移

如下内容主要来自Version 4 to 5/6 API Migration

迁移到5.x版本

这是但愿将代码迁移到使用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 函数 新版 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()

迁移到6.x版本

这是但愿迁移代码以使用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 函数 新版 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()
相关文章
相关标签/搜索