【GIS开发】osgEarth依赖库PROJ（Python）

PROJ(Cartographic Projections and Coordinate Transformations Library)，
地图投影的程序库，遵循X/MIT开源协议，用于坐标转换和坐标参考系统转换。

1、OSGeo/PROJ（C++）

https://proj.org/
https://github.com/OSGeo/PROJ
https://download.osgeo.org/proj

PROJ - Cartographic Projections and Coordinate Transformations Library

PROJ 是一款通用坐标转换软件，可将地理空间坐标从一个坐标参考系 (CRS) 转换为另一个坐标参考系 (CRS)。这包括制图投影和大地测量变换。PROJ 在X/MIT 开源许可下发布

PROJ 包括命令行应用程序，用于轻松转换来自文本文件或直接来自用户输入的坐标。除了命令行实用程序之外，PROJ 还公开了一个 应用程序编程接口，简称 API。API 允许开发人员在他们自己的软件中使用 PROJ 的功能，而无需自己实现类似的功能。

PROJ 最初只是作为一个制图应用程序，让用户可以使用许多不同的制图投影将大地坐标转换为投影坐标。多年来，随着需求变得明显，对基准转换的支持也慢慢进入 PROJ。今天，PROJ 支持一百多种不同的地图投影，并且可以使用除了最晦涩的大地测量技术之外的所有方法来转换基准面之间的坐标。

PROJ库的编译依赖库需求如下：

• C99 compiler
• C++11 compiler
• CMake >= 3.9
• SQLite3 >= 3.11: headers and library for target architecture, and sqlite3 executable for build architecture.
• libtiff >= 4.0 (optional but recommended)
• curl >= 7.29.0 (optional but recommended)

1.1 编译sqlite3

https://sqlite.org/download.html
这里静态编译SQLite库的64位版本。

将下载好的两个文件sqlite-amalgamation-3380100.zip、和sqlite-dll-win64-x64-3380100.zip解压到SQLite文件夹内,如下所示：
sqlite-amalgamation-3380100： shell.c 、sqlite3.c、sqlite3.h、sqlite3ext.h
sqlite-dll-win64-x64-3380100.zip：sqlite3.def、sqlite3.dll


打开VS2017，新建一个静态库工程，如下：

将下载的sqlite3的相关，添加到工程里，然后设置“不使用预编译头”，编译工程如下：

C/C++ --> 预处理器 --> 预处理器定义：添加预定义处理如下：

_USRDLL
SQLITE_ENABLE_RTREE
SQLITE_ENABLE_COLUMN_METADATA
SQLITE_ENABLE_FTS5
SQLITE_ENABLE_UNLOCK_NOTIFY

设置模块定义文件,链接器 --> 输入 --> 模块定义文件：

sqlite3.def

修改模块定义文件，在它内容最后面添加：sqlite3_unlock_notify
然后配置类型改为静态库lib。

1.2 编译libtiff

http://download.osgeo.org/libtiff/

1.3 编译openssl

https://www.openssl.org/
https://github.com/openssl/openssl
http://slproweb.com/products/Win32OpenSSL.html

openssl源码编译需要perl环境，配置比较麻烦。这里直接采用网上编译好的文件（http://slproweb.com/products/Win32OpenSSL.html），如下图所示：

下载安装之后，文件夹如下：

1.4 编译curl

https://curl.se/download.html
https://github.com/curl/curl

1.5 编译PROJ9

官网地址：
http://download.osgeo.org/proj/
从官网上下载PROJ最新代码：proj-9.0.0.tar

解压文件proj-9.0.0.tar后，可以看见cmake的脚本文件CMakeLists.txt，如下图：

于是我们使用cmake构建PROJ库的VS2017的工程文件，进而编译代码。

2、pyproj（python库）

2.1 概述

Python interface to PROJ (cartographic projections and coordinate transformations library).

• Minimum supported PROJ version is 8.0
• Minimum supported Python version is 3.8

2.2 安装

https://github.com/pyproj4/pyproj

pip install pyproj

有时候通过pip官方直接安装第三方库会失败，使用国内镜像相对稳定些。
pipy国内镜像目前有：

http://pypi.douban.com/  豆瓣
http://pypi.hustunique.com/  华中理工大学
http://pypi.sdutlinux.org/  山东理工大学
http://pypi.mirrors.ustc.edu.cn/  中国科学技术大学

阿里云：

pip install robotframework -i http://mirrors.aliyun.com/pypi/simple --trusted-host mirrors.aliyun.com 阿里云

豆瓣：

pip install pyproj -i http://pypi.douban.com/simple/ --trusted-host pypi.douban.com

pip list

2.3 代码测试

• Initializing CRS

from pyproj import CRS
crs = CRS.from_epsg(4326)
crs = CRS.from_string("epsg:4326")
crs = CRS.from_proj4("+proj=latlon")
crs = CRS.from_user_input(4326)

• Transformations from CRS to CRS

from pyproj import CRS
crs_4326 = CRS.from_epsg(4326)
crs_26917 = CRS.from_epsg(26917)
print(crs_26917)

from pyproj import Transformer
transformer = Transformer.from_crs(crs_4326, crs_26917)
transformer = Transformer.from_crs(4326, 26917)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:26917")
print(transformer)

• Geodesic line length

from pyproj import Geod
lats = [-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
        -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7]
lons = [-74, -102, -102, -131, -163, 163, 172, 140, 113,
        88, 59, 25, -4, -14, -33, -46, -61]
geod = Geod(ellps="WGS84")
total_length = geod.line_length(lons, lats)
print(f"{total_length:.3f}")

• Geodesic area

from pyproj import Geod
geod = Geod('+a=6378137 +f=0.0033528106647475126')
lats = [-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
        -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7]
lons = [-74, -102, -102, -131, -163, 163, 172, 140, 113,
        88, 59, 25, -4, -14, -33, -46, -61]
poly_area, poly_perimeter = geod.polygon_area_perimeter(lons, lats)
print(f"{poly_area:.3f} {poly_perimeter:.3f}")

• WGS84坐标系下，计算该经纬度对应的（x,y）
  https://pyproj4.github.io/pyproj/stable/api/proj.html

from pyproj import Proj
p = Proj(proj='utm',zone=10,ellps='WGS84', preserve_units=False)
x,y = p(-120.108, 34.36116666)
# lonlat to xy
print('x=%9.3f y=%11.3f' % (x,y))
# xy to lonlat
print('lon=%8.3f lat=%5.3f' % p(x,y,inverse=True))

lons = (-119.72,-118.40,-122.38)
lats = (36.77, 33.93, 37.62 )
x,y = p(lons, lats)
print('x: %9.3f %9.3f %9.3f' % x)
print('y: %9.3f %9.3f %9.3f' % y)

lons, lats = p(x, y, inverse=True) # inverse transform
print('lons: %8.3f %8.3f %8.3f' % lons)
print('lats: %8.3f %8.3f %8.3f' % lats)

p2 = Proj('+proj=utm +zone=10 +ellps=WGS84', preserve_units=False)
x,y = p2(-120.108, 34.36116666)
print('x=%9.3f y=%11.3f' % (x,y))

p = Proj("epsg:32667", preserve_units=False)
print('x=%12.3f y=%12.3f (meters)' % p(-114.057222, 51.045))

后记

如果你觉得该方法或代码有一点点用处，可以给作者点个赞、赏杯咖啡；╮(￣▽￣)╭
如果你感觉方法或代码不咋地//(ㄒoㄒ)//，就在评论处留言，作者继续改进。o_O???
谢谢各位童鞋们啦( ´ ▽´ )ﾉ ( ´ ▽´)っ！！！