Shapefile 投影文件

Posted 2023-03-23

【中文标题】Shapefile 投影文件

【英文标题】：Shapefile projection files

【发布时间】：2014-10-07 05:30:45

【问题描述】：

我正在处理一些涉及 Shapefile 的开放数据，并且我设法将数据转换为 GeoJSON 并提取点的坐标（以及其他信息）。但是坐标似乎不是纬度/经度坐标：

[ 368667.2455000002, 5009510.4067 ]
[ 367885.47140000015, 5019804.3237 ]
[ 395852.80260000005, 5027699.9354 ]
[ 379358.1364000002, 5036798.5747 ]
[ 351968.9621000001, 5017404.8727 ]
[ 375123.64269999973, 5033338.467499999 ]
[ 378133.7736999998, 5032617.237500001 ]
[ 385791.0351, 5010557.9144 ]
[ 349796.77770000044, 5013571.1559999995 ]
[ 367271.0566999996, 5030212.897399999 ]
[ 378808.0292999996, 5013921.017100001 ]
[ 336650.69820000045, 5039983.886399999 ]
[ 364178.05599999987, 5015957.9625 ]
[ 362715.9132000003, 5023112.371200001 ]
[ 351321.0865000002, 5013763.298699999 ]
[ 373254.0789000001, 5026533.7453000005 ]
[ 355235.6211000001, 5016957.0644000005 ]
[ 327797.9938000003, 5036758.195699999 ]
[ 362836.86930000037, 5015895.811000001 ]
[ 375479.41530000046, 5001245.2017 ]

经过进一步研究，我发现我必须根据 PRJ 文件中的数据“投影”这些坐标：

PROJCS["City of Ottawa",
    GEOGCS["GCS_North_American_1983",
        DATUM["D_North_American_1983",
            SPHEROID["GRS_1980", 6378137.0, 298.257222101]
        ],
        PRIMEM["Greenwich", 0.0],
        UNIT["Degree", 0.0174532925199433]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting", 304800.0],
    PARAMETER["False_Northing", 0.0],
    PARAMETER["Central_Meridian", -76.5],
    PARAMETER["Scale_Factor", 0.9999],
    PARAMETER["Latitude_Of_Origin", 0.0],
    UNIT["Meter", 1.0]
]

我已经读过What does an esri projection file represent? 和随后的https://gis.stackexchange.com/questions/2383/does-re-projecting-an-esri-shapefile-only-update-the-content-of-the-prj-file/2386#2386。尽管它们描述了投影文件的作用/是什么，但没有讨论如何完成转换。

与 Shapefile 不同，我无法找到有关 PRJ 文件的规范。是否有用于将坐标转换为纬度/经度坐标的算法？

我也愿意使用库，但它们需要与 android（GeoTools 不适合我）或 javascript（Node）兼容。

【问题讨论】：

这么多年过去了，你有更新吗？如何使用 nodeJS 或 JS 使用该投影文件将这些坐标转换为 GPS 坐标（wgs84）？

【参考方案1】：

如何实际完成转换可能非常复杂，尤其是在基准面之间进行转换时（例如 NAD83 到 WGS84 等）。有很多工具可以做到这一点，但我不是熟悉移动设备，但我知道有些人已经让 GDAL 和 OGR 在 Android 上构建。见http://trac.osgeo.org/gdal/wiki/BuildingForAndroid。

您需要找到数据的 EPSG 代码。根据您发布的内容，这些坐标位于 EPSG:32189（MTM 区域 9）中。为了找到这个，我使用了网站http://prj2epsg.org，您可以在其中上传 .prj 文件，它会为您找到 EPSG 代码。您可以在 spatialreference.org 查看有关该投影的更多信息。该站点有几个有用的功能，例如每个投影的可下载 prj 文件，以及用于上述 proj4 库的 proj4 文件。

希望这会有所帮助。

【讨论】：

通常，NAD83 和 WGS84 之间的距离在一米以内。出于实际目的，它们可以被认为是相同的。

【参考方案2】：

根据最新的 GeoJSON 规范 (RFC7946)，CRS 必须是 WGS84(lat/lng)。 GDAL 套件等工具假定您的输出的 CRS 位于 WGS84 中。

也就是说，重投影步骤必须在从 shapefile 转换为 GeoJSON 时或之前完成。

例如，来自GDAL的ogr2ogr，并且假设PRJ伴随文件与shapefile位于同一目录中（将从PRJ文件中读取坐标系，无需指定输入EPSG代码），

ogr2ogr -f GeoJSON output.json input.shp

请注意，在最近的 GDAL 版本中，不需要强制输出坐标系为 WGS84，如第一段所述。以前的语法是，

ogr2ogr -t_srs EPSG:4326 -f GeoJSON output.json input.shp

您可以为ogr2ogr 使用 NodeJS 包装器

例如，

https://github.com/wavded/ogr2ogr

对应的npm包，

https://www.npmjs.com/package/ogr2ogr

编辑 1：

或者，可以使用 SWIG 为 GDAL 生成的 NodeJS 绑定，

https://www.npmjs.com/package/gdal

编辑 2：

尽管如此，如果 OP 真的想将仍在 EPSG:32198 中的 GeoJSON 转换为 WGS84 中的 GeoJSON，这将完成这项工作：

ogr2ogr -t_srs EPSG:4326 -s_srs EPSG:32198 -f GeoJSON output.json input.json

请注意，GDAL 与 EPSG 数据库一起分发，将从中读取以下字符串：

+proj=tmerc +lat_0=0 +lon_0=-76.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

对应问题中的 PRJ 文件。

【参考方案3】：

使用 Javascript，我发现并实现的最佳解决方案（它是 working）依赖于使用库 proj4。但我知道这个库也适用于 Python。

在 Javascript 中是：

var firstProjection = 'PROJCS["NAD83 / Massachusetts Mainland",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",42.68333333333333],PARAMETER["standard_parallel_2",41.71666666666667],PARAMETER["latitude_of_origin",41],PARAMETER["central_meridian",-71.5],PARAMETER["false_easting",200000],PARAMETER["false_northing",750000],AUTHORITY["EPSG","26986"],AXIS["X",EAST],AXIS["Y",NORTH]]';
var secondProjection = "+proj=gnom +lat_0=90 +lon_0=0 +x_0=6300000 +y_0=6300000 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
//I'm not going to redefine those two in latter examples.
proj4(firstProjection,secondProjection,[2,5]);
// [-2690666.2977344505, 3662659.885459918]

如果只给出 1 个投影，则假定它是从 WGS84（GPS 标准）投影的。

proj4(firstProjection,[-71,41]);
// [242075.00535055372, 750123.32090043]

【讨论】：

GDAL 在木头下使用 PROJ。有几个 SWIG 生成的 PROJ 绑定，其中之一是 PyProj（python 绑定）。 Proj4js 是用 JavaScript 重写的部分 PROJ。