如何转换地图 Window 和数据点的 CRS 以匹配 SF 对象？

Question

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

我正在使用 rnaturalearth 的 countries 数据作为我的 spatialpolygondataframe。我认为这个数据的原始投影 (EPSG: 4326 proj4string: "+proj=longlat +datum=WGS84 +no_defs") 看起来不太好，所以我将 countries 数据转换为 Lambert正形圆锥投影 (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')。我不确定这是最好的投影，但我认为这是一个开始的地方，我愿意接受建议。

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

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

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())

Answer 1

恐怕这是经典之作。你的纬度和经度是错误的方式。纬度和经度是人们所说的通常顺序 -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)

现在应该可以了！

如何转换地图 Window 和数据点的 CRS 以匹配 SF 对象？

How to Transform CRS of Map Window and Datapoints to Match SF Object?

r

ggplot2

sf

proj