R语言中的地理/投影坐标系统（上）[翻译]

最新推荐文章于 2022-11-21 12:16:23 发布
最新推荐文章于 2022-11-21 12:16:23 发布
阅读量2.3k 收藏 17
点赞数 1
文章标签： python java 编程语言 大数据 mysql
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_54000907/article/details/124977348
版权

原文链接：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）
+zoneUTM投影的带

投影类型（即参数+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）中查看。

2ec584d93cbc19c3f9c3f85526148e05.png

这个编码可以放到网址http://spatialreference.org/ref/的搜索框里，然后查看proj4参数（注意必须在语法选项清单中选择Proj4）。

30acedb21e1539b6c39ec3194ab4f6cf.png

下面是一些常见的投影：

投影WKID来源语法
UTM NAD 83 Zone 19N26919EPSG+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
USA Contiguous albers equal area102003ESRI+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 area3338EPSG+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 Robinson54030ESRI+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
R语言学堂
关注
  • 1
    点赞
  • 踩
  • 17
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
【R-terra包|第二节】：栅格重分类、投影变换、分区统计
地理空间科学的博客
05-06 5201
基于R语言terra的空间运算
地理、投影坐标转换工具
04-28
坐标转换工具，能够进行地理坐标和各种投影之间的良好转换，可以一个一个转化，并支持批量转换。
R语言绘制两种不同形式的投影方法效果图
临风暖阳的博客
11-20 449
笔者用R语言绘制两种不同形式的投影方法效果图
raster | R语言中的空间栅格对象及其基本处理方法（Ⅱ）：投影、属性提取
R语言学堂
07-25 7103
library(raster) 投影设置投影在使用raster()函数创建栅格对象时，调用crs参数可以定义其投影。投影参数具有默认值：r<-raster() crs(r) ##...
Openlayers使用Proj4js坐标转换cgcs2000分带（zone）投影的问题
GIS小虫的博客
05-12 3512
问题源于使用4528投影按照官网示例定义view时地图一片空白，于是开始了对Proj4js和投影转换的抽丝剥茧，这里记录下了解的知识和最后的结论。 官网示例和Proj4js的使用 OL示例见wms-image-custom-proj 这里简单介绍下proj4，我认知的proj4是基于一些约定参数组成的用于定义坐标参考的一组参数，和WKT相比更简约，省略了WKT里的完整参数定义。比如4528的wkt PROJCRS["CGCS2000 / 3-degree Gauss-Kruger zone 40", .
GIS坐标系、投影和转换
08-05
GIS（Geographic Information System，地理信息系统）是处理和分析地理空间数据的重要工具，其中坐标系、投影和转换是GIS中的核心概念。理解这些概念对于进行有效的地理数据分析至关重要。 坐标系是GIS的基础，它...
vb.rar_vb坐标正算_高斯_高斯投影 VB
09-19
这个“vb.rar_vb坐标正算_高斯_高斯投影 VB”压缩包内容显然涉及到使用VB（Visual Basic）编程语言实现高斯投影坐标的正算算法。下面我们将深入探讨高斯投影、坐标正算以及如何用VB来实现这一过程。 一、高斯投影 ...
经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法.rar
01-11
在地理信息系统和测绘领域，高斯投影是一种广泛使用的地图投影方法，它将地球表面的球面坐标（经度和纬度）转换为平面上的直角坐标（X和Y），以便于地图制作和分析。本资源包含的是经纬度BL转换到高斯平面直角坐标XY...
学习利用开源GIS完成数据录入、服务发布与展示（二）
weixin_42512595的博客
03-20 294
一、安装QGis，试验导入shape file。 QGIS导入数据的地方比较特别，在数据库/数据管理器中有导入选项。 在导入shape file时报错：试验多次后，发现用qgis打开查看属性显示无效投影，在坐标参考系部分可以看到是值是“Unknown CRS: PROJCRS["NAD83 / North Carolina (ftUS)”。可以用qgis导出，设置空间参考就可以顺利导入了。 ...
R语言中的地理/投影坐标系统（下）[翻译]
R语言学堂
05-28 2088
原文链接：https://mgimond.github.io/Spatial/coordinate-systems-in-r.html。译文分上、下两篇，这里为下篇。❝「译者注」：在原文的本部分出现了许多世界地图，并且地图涉及到了我国领土边界的标识，但是这种标识是错误的。为了避免这一错误并尽可能地保留原文面貌，在确实需要标记区域边界时，删去了亚洲国家和地区；没有必要标识边...
R语言：栅格数据处理
hu397313168的博客
11-21 2075
>％来自dplyr包的管道函数，其作用是将前一步的结果直接传参给下一步的函数，从而省略了中间的赋值步骤，可以大量减少内存中的对象，节省内存。然后把它转化为空间信息数据框，这里需要声明经纬所在列，以及坐标系定义为4326（WGS 84）本帖子会简单介绍如何读入并处理栅格数据。从世界气温图中裁剪出肯尼亚的部分，使用crop函数。插值之前让我们先统一坐标系，然后转化为sp对象。空间插值之前，先合并气温的时空分布数据。，用到的是澳洲的地图。裁剪出真正属于肯尼亚的部分。一次读入所有栅格数据。
学习mysql4.1的空间数据库特性(1)
飞羽飞の猪の博客
10-16 9172
06年10月16日需求:公司现在gis系统有大量的经纬度坐标，都是车辆行驶几个月采集过来的坐标。现在在做区域搜索功能（就是在上个月某个时间内车辆在哪里行驶）。因为现在数据库单个车辆的经纬度坐标太多，所以搜索起来非常慢（让人忍受不了的速度），随着时间的推移，可能搜索起来会越来越慢，所以公司决定考虑使用mysql的空间数据库特性。把经纬度坐标都转移到成空间数据库坐标. 公司使用mysql4
R语言做地理分析：数据导入，合并，绘图，导出
weixin_38980061的博客
09-28 6380
学习笔记，原文链接????：https://andrewmaclachlan.github.io/CASA0005repo/geographic-information.htmlhttps://andrewmaclachlan.github.io/CASA0005repo/geographic-information.html 0 数据准备 0.1 伦敦边界数据 (.shp) 0.2 伦敦2008-2020垃圾乱丢事件统计数据 (.csv) https://pan.baidu.com/s/1mBri
R语言空间数据处理笔记—矢量-栅格数据交互
热门推荐
万里写入胸怀间
06-27 1万+
矢量-栅格数据交互 根据矢量数据范围裁剪栅格数据 统一投影 vector = st_transform(vector,projection(raster)) 设定矢量数据范围外栅格数据的值为NA raster_masked= mask(raster,as(vector,'Spatial') 本Markdown编辑器使用StackEdit修改而来，用它写博客，将会带来全新的体...
基于R语言进行栅格数据统计及Raster包简介
撼沧
05-25 1万+
基于R语言进行栅格数据统计及Raster包简介 写在前面：日常数据处理中常需要对单个栅格图层进行统计，包括计算其均值、最值和总和等等。使用R语言可方面快速地实现该需求，主要借助的package为Raster包。 library(raster) x = raster('data/data.tif') cellStats(x,stat='mean') x表示待计算的栅格； stat表示需要统计的内容，包括sum, mean, min, max, sd（标准差）, skew（偏度） 和 rms（均方根）；
R语言空间数据分析学习笔记3——坐标系设置
GeoSuper的博客
11-30 7583
地理空间数据一般都带有坐标系，最常用的是WGS 84，如果我们已经读入了一个地理数据，比如保存在data中，那么可以查看它的空间坐标系： library(pacman) p_load(sf) data = st_read("G:/Rdata/China/fujian1.shp") # Reading layer `fujian1' from data source `G:\Rdata\China\fujian1.shp' using driver `ESRI Shapefile' # Simple feat
R绘图系统中的坐标系
hsdcc217的博客
12-03 2571
在R语言中的坐标系
QGIS WKT自定义坐标系研究及其与ArcGIS的对比
QGISClass的博客
01-08 4346
我国全国版图常用的投影为双标准纬线圆锥投影，如等积Albers投影或者等角Lambert投影，所使用的地理坐标参照系则依据数据生产时间而有所不同，如早期的北京54坐标系和西安80坐标系，以及近年来所普遍采用的国家2000大地坐标系。本文以获取自中国科学院资源环境科学与数据中心的2015年省级行政区划数据为示范，该数据为SHP格式，采用的就是基于北京54坐标系的Albers投影。 这种双标准纬线、根据制图区域选定中央经线、地理坐标系、投影方式的组合，在大多数GIS软件的预定义坐标系中普遍不存在，需要自.
R实现地理位置与经纬度相互转换
weixin_34367845的博客
01-11 960
本实例要实现目标通过输入城市名或者地名，然后找出其经度纬度值，以及通过可视化展现其线路流向以及周边地图展示 address_list数据： 山西省太原市小店区亲贤北街77号 贵州省贵阳市云岩区书香门第B栋3单元 北京市通州区神树商业街168号 贵州省贵阳市南明区兴关路51 北京市东城区长巷二条乙5号 山西省太原市杏花岭区北肖墙12号 北京 北京市通州区 北京市的的的的的 天津...
地理坐标系转投影坐标
最新发布
07-27
地理坐标系转投影坐标是将地球表面上的地理坐标（经度和纬度）转换为平面上的投影坐标（通常是直角坐标）。这个过程称为地理坐标系到投影坐标系的转换。 在进行地理坐标系到投影坐标系的转换时，通常需要选择一个特定的投影方法或坐标系统。常见的投影方法包括经纬度坐标系统、UTM（通用横轴墨卡托投影）、高斯-克吕格投影等。 具体的转换方法和工具可能会因使用的地理信息系统（GIS）软件或编程语言而有所不同。一般而言，可以使用专门的地理信息系统软件，如ArcGIS、QGIS等，来进行地理坐标系和投影坐标系之间的转换。此外，也可以使用编程语言中的库或函数来实现这一转换，如Python中的pyproj库、JavaScript中的proj4库等。 在进行地理坐标系到投影坐标系的转换时，需要提供相关的参数，如选择的投影方法、中央经线、假东原点等。这些参数的选择取决于具体应用场景和需求。 总之，地理坐标系到投影坐标系的转换是地理信息处理中常见且重要的步骤之一，可以将地球上的地理位置信息转换为能够在平面上进行分析和可视化的投影坐标。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值