原文链接：https://mgimond.github.io/Spatial/coordinate-systems-in-r.html。
译文分上、下两篇，这里为上篇。
1 「关于PROJ环境更新的说明」
新版本的sf工具包使用的是PROJ 6.0C库或更高版本。注意不要把PROJ的版本与R里的proj4工具包弄混淆——proj4和sf都是依赖于PROJ C库构建的R工具包，而PROJ C库是独立于R环境的。关于PROJ的更多信息可以查看https://proj.org/。
PROJ 6.0相比于之前的版本有明显的变化。它对R的空间分析生态造成了很大的影响。如果你使用的proj4和sf工具包是依赖PROJ 6.0之前的版本构建的话，那么接下来本文的一些输入/输出结果可能会与你得到的不太一样。
2 「本练习的示例数据」
通过如下代码可以将示例数据加载到R环境中：
load(url("https://github.com/mgimond/Spatial/raw/main/Data/Sample1.RData"))本文仅会使用其中的两个数据：缅因州的县级矢量图层s.sf；海拔栅格数据elev.r。前者是sf格式的，后者是raster格式的。其他数据可以通过下面的代码删除：
rm(list=c("inter.sf", "p.sf", "rail.sf"))3 「加载sf工具包」
library(sf)注意GEOS、GDAL和PROJ等工具包的版本会与sf关联。这些工具包的版本不同可能会使输出结果与本文不同。你可以通过下面代码查看相关包的版本号：
sf_extSoftVersion()[1:3]
## GEOS GDAL proj.4
## "3.9.1" "3.2.1" "7.2.1"❝❞
译者注：上述版本号为译者运行的结果。原文
GEOS的版本号是3.9.0，其余两个与这里相同。下文若无特别说明，输出内容都是译者运行的结果。
4 「查看一个坐标系统」
提取一个sf对象的坐标系统（coordinate system，CS）所使用的函数是st_crs()：
st_crs(s.sf)
## Coordinate Reference System:
## User input: EPSG:26919
## wkt:
## PROJCRS["NAD83 / UTM zone 19N",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["UTM zone 19N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-69,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["North America - between 72<U+00C2>°W and 66<U+00C2>°W - onshore and offshore. Canada - Labrador; New Brunswick; Nova Scotia; Nunavut; Quebec. Puerto Rico. United States (USA) - Connecticut; Maine; Massachusetts; New Hampshire; New York (Long Island); Rhode Island; Vermont."],
## BBOX[14.92,-72,84,-66]],
## ID["EPSG",26919]]在新版本的PROJ C库中，坐标系统的定义格式是Well Known Text (WTK/WTK2)，它是由一系列的[...]标签组成的。对于投影坐标系统，WTK通常以PROJCRS[...]开头；对于地理坐标系统则以GEOGCRS[...]标签开头。
CRS输出结果还可以包括用户定义的EPSG代码（如本例），或者由大地基准面（datum）和投影类型（projection type）组成的字符串。
提取栅格对象的坐标系统所使用的是raster包中的crs()函数：
library(raster)
crs(elev.r)
## CRS arguments:
## +proj=utm +zone=19 +datum=NAD83 +units=m +no_defs +ellps=GRS80
## +towgs84=0,0,0注意，在该版本的raster工具包中，crs()函数不会输出格式化的WKT字符串。
截至目前，共有两种方式来定义坐标系统：一种是使用EPSG代码（https://spatialreference.org/ref/epsg/）；另一种是使用PROJ4格式化的字符串（https://proj.org/apps/proj.html）。这两种方式在sf和raster工具包中都可以使用。
在新版本的PROJ C库中，你也可以使用Well Known Text (WTK/WTK2)去定义sf对象的坐标系统。这种格式语法更加复杂（正如上文st_crs()函数的输出结果），并且对于手动定义来说，它也不是一个简单的方法。如有可能，请使用权威机构定义过的EPSG代码。不过，如果要定制一个坐标系统，使用PROJ4语法会更加容易。
5 「理解PROJ4的语法」
PROJ4语法是由一系列参数组成的，每个参数均以+为前缀。比如，elev.r的坐标系统使用的是通用墨卡托投影（UTM projection，+proj=utm）的第十九个带（zone19，+zone=19），大地基准面为NAD 1983（+datum=mad83）。
下面列举了一些PROJ4的参数，完整参数请查看网址https://proj.org/usage/projections.html。
| 参数 | 描述 |
|---|---|
| +a | 椭球体的半长轴 |
| +b | 椭球体的半短轴 |
| +datum | 大地基准面名称 |
| +ellps | 椭球体名称 |
| +lat_0 | 原点（origin）的纬度坐标 |
| +lat_1 | 第一个标准纬度平行线 |
| +lat_2 | 第二个标准纬度平行线 |
| +lat_ts | 真尺度纬度（Latitude of true scale） |
| +lon_0 | 中央子午线 |
| +over | 允许经度值在-180-180范围外 |
| +proj | 投影名称 |
| +south | 表示UTM投影南半球的带（zone） |
| +units | 距离单位，可以为米、英尺等 |
| +x_0 | 东移假定值（False easting） |
| +y_0 | 北移假定值（False northing） |
| +zone | UTM投影的带 |
投影类型（即参数+proj）可在网址https://proj.org/operations/projections/查看。
6 「定义一个坐标系统」
定义过的坐标系统能够赋予空间对象。它可以赋给没有坐标参考系的对象，也可以覆盖已存在的坐标系统（后者只有在能确定对象原有的坐标系统是错误的情况下才会使用）。注意，本节的操作不会改变对象背后的坐标值（这个过程留在下节介绍）。
我们假设s.sf原本不存在坐标参考系，然后使用st_set_crs()函数进行手动指定。在下面的例子中，我们使用proj4语法定义坐标系统：
s.sf <- st_set_crs(s.sf, "+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83")现在让我们查看对象的坐标系统：
st_crs(s.sf)
## Coordinate Reference System:
## User input: +proj=utm +zone=19 +ellps=GRS80 +datum=NAD83
## wkt:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6269]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 19N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-69,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16019]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]你可以看到，在User input字段显示的是调用st_set_crs()函数时的proj4字符串。但是，你也会注意到，在WKT字符串中，一些参数比如PROJCRS、BASEGEOGCRS显示是未知的（"unknown"）。这并不是一个太大的问题，因为关键的大地基准面和投影信息是存在的。尽管如此，在可能的情况下，使用EPSG代码定义坐标系统会是个不错的选择：
UTM NAD83 Zone 19N的EPSG代码是26919：
s.sf <- st_set_crs(s.sf, 26919)让我们再来看看坐标系统：
st_crs(s.sf)
## Coordinate Reference System:
## User input: EPSG:26919
## wkt:
## PROJCRS["NAD83 / UTM zone 19N",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["UTM zone 19N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-69,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["North America - between 72掳W and 66掳W - onshore and offshore. Canada - Labrador; New Brunswick; Nova Scotia; Nunavut; Quebec. Puerto Rico. United States (USA) - Connecticut; Maine; Massachusetts; New Hampshire; New York (Long Island); Rhode Island; Vermont."],
## BBOX[14.92,-72,84,-66]],
## ID["EPSG",26919]]关键的投影参数仍然是一致的。但是在WKT字符串开头还额外添加了一些信息。
你同样可以使用前面定义s.sf坐标系统的PROJ4字符串去定义栅格数据的投影，使用的函数是crs()（这里我们假定空间对象缺乏坐标信息或坐标信息不准确）。
crs(elev.r) <- "+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83"注意，我们不需要定义所有的投影参数，只要我们知道未被定义的参数的默认值是正确的就可以了。另外，我们也不需要指定半球，因为NAD83大地基准面仅适用于北美。
让我们来查看栅格对象的坐标系统：
crs(elev.r)
## CRS arguments:
## +proj=utm +zone=19 +datum=NAD83 +units=m +no_defs❝❞
译者注：此处输出结果为作者原文。
使用EPSG代码定义栅格坐标系统的方法如下：
crs(elev.r) <- "+init=EPSG:26919"
crs(elev.r)
## CRS arguments:
## +proj=utm +zone=19 +datum=NAD83 +units=m +no_defs请注意EPSG代码是如何转换成proj4字符串的。
为了重现一些软件如ArcGIS中已定义坐标系统的参数，最好的方法是先提取它们的WKID/EPSG代码，然后将数字放到网站http://spatialreference.org/ref/上进行检索。例如，在ArcGIS中，WKID编码可以在坐标系统属性（coordinate system properties）中查看。
这个编码可以放到网址http://spatialreference.org/ref/的搜索框里，然后查看proj4参数（注意必须在语法选项清单中选择Proj4）。
下面是一些常见的投影：
| 投影 | WKID | 来源 | 语法 |
|---|---|---|---|
| UTM NAD 83 Zone 19N | 26919 | EPSG | +proj=utm +zone=19 +ellps=GRS80 +datum=NAD83 +units=m +no_defs |
| USA Contiguous albers equal area | 102003 | ESRI | +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs |
| Alaska albers equal area | 3338 | EPSG | +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs |
| World Robinson | 54030 | ESRI | +proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs |
扫一扫
6440
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
