时间戳

首页 > 解决方案 > 如何转换地图窗口和数据点的 CRS 以匹配 SF 对象？

r - 如何转换地图窗口和数据点的 CRS 以匹配 SF 对象？

问题描述

我正在尝试创建一张显示站点位置的北美太平洋海岸（加拿大南部、美国、下加利福尼亚州）的地图。理想情况下，地图投影将保留海岸线的形状/方向。

我使用rnaturalearth'scountries数据作为我的spatialpolygondataframe. 我认为这个数据的原始投影（EPSG: 4326 proj4string: "+proj=longlat +datum=WGS84 +no_defs"）看起来不太好，所以我将countries数据转换为 Lambert Conformal Conic projection (proj4string:+proj= lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs')。我不确定这是最好的投影，但我认为这是一个开始的地方，我愿意接受建议。

我正在努力绘制经过改造的科幻小说来展示北美的太平洋海岸线。我知道我需要转换我的窗口坐标以匹配 LCC 北美投影以及我的站点坐标，但是 Rstudio 给了我一个空白绘图窗口和错误。

我尝试按照 [ https://www.r-bloggers.com/zooming-in-on-maps-with-sf-and-ggplot2/][1]中列出的步骤进行操作，但经过多次尝试都没有成功并尝试其他方法。我在转换我的地图窗口时哪里出错了？这是我第一次尝试映射变换投影。任何建议将不胜感激！提前致谢！

library(tidyverse)
library(ggplot2)
library(sf)
library(rnaturalearth) # map data source
library(rgdal)

## Download sf from rnaturalearth
country <- ne_download(scale = 10, type = 'countries', category = 'cultural', returnclass = 'sf')
st_crs(country) # show CRS of data
st_proj_info() # lcc is available

# Lambert Conformal Conic projection proj4 according to https://epsg.io/102009
LCC_North_America <- '+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
# transform country data to LCC N. America
country_LCC <- st_transform(country, crs = LCC_North_America)

## site locations, need to transform to LCC N.America and plot
lat <- c(48.98604, 47.72322, 37.96389, 33.61817)
long <- c(-122.7593, -122.6559, -122.4881, -117.9052)

disp_wind_4326 <- st_sfc(st_point(c(-127,49)), st_point(c(-112, 29)), crs = 4326)
disp_wind_LCC <- st_transform(disp_wind_4326, crs = LCC_North_America)
disp_window <- st_coordinates(disp_wind_LCC)

(NOOC_coast_fig <- ggplot() +
   geom_sf(data = country_LCC) +
   coord_sf(xlim = disp_window[,'X'], ylim = disp_window[,'Y'],
            datum = LCC_North_America, expand = FALSE) +
   theme_bw())

标签： rggplot2sfproj

解决方案

很经典，我怕。你的纬度和经度是错误的方式。纬度和经度是人们所说的通常顺序-Y然后X，这无济于事。

disp_wind_4326 <- st_sfc(st_point(c(-127, 49)), st_point(c(-112, 29)), crs = 4326)

对于自助，值得检查对象的结构以查看它们是否符合您的预期。当我运行您的代码时，我查看disp_wind_LCC并看到了这个：

> disp_wind_LCC
Geometry set for 2 features  (with 2 geometries empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    102009
proj4string:    +proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
POINT EMPTY
POINT EMPTY

此时很明显，输入出现了问题，因为它们没有创建有效的几何图形。您不能在 (-90, 90) 之外有纬度，因此它给出了一个空点。

根据您在下面的评论，下一个问题非常模糊。您正在将全球数据转换为您的 LCC 投影。如果您要绘制整个country_LCC数据集，您会发现存在大量失真。

特别是，南极洲根本无法应对。这很常见：如果您提供的本地投影地理坐标远远超出其使用的域，则可能会出现严重错误。

> st_bbox(country_LCC[country_LCC$ADMIN == 'Antarctica',])
      xmin       ymin       xmax       ymax 
-335167885 -226033419   42202752   22626164 
> st_bbox(country_LCC[country_LCC$ADMIN == 'United States of America',])
    xmin     ymin     xmax     ymax 
-6653437 -1523147  2124650  4338858

所以你的情节是正确的，但是有一个错误的南极洲多边形位于其他一切之上。因此，在这里和一般情况下，在进行任何重新投影之前，将全局数据减少到感兴趣的区域：

country <- subset(country, ADM0_A3 %in% c("USA", "CAN", "MEX"))
country_LCC <- st_transform(country, crs = LCC_North_America)

现在应该可以了！

推荐阅读