将 sf 对象裁剪到 R 中另一个 sf 多边形的范围的问题

problem cropping sf object to extent of another sf polygon in R

我经常下载大尺度的气候数据,并将其裁剪到某个空间对象的范围内。我的代码仍然可以通过将海面温度数据裁剪为 "bathymetric mask"(一种落在边界框内并且比海洋中的某个深度浅的形状)来执行此操作。我从 here ("sst.mon.mean.nc") 下载了一个新的温度数据集,但我不知道如何正确裁剪它。我尝试过的一切都会产生一个空的 sf 对象。

这里是重现问题的代码,下载了上面链接的数据集:

library(tidyverse)
library(raster)
library(sf)
library(oceanmap)

# set bounding box
latrange <- c(35, 45)
lonrange <- c(-78, -66) 

# download COBE data and put in working directory, then run: 
cobe.raw <- raster::stack("sst.mon.mean.nc")
cobe.tidy <- raster::as.data.frame(cobe.raw, xy=TRUE)  %>% 
  pivot_longer(cols=c(-x, -y), names_to="date", values_to="sst") %>%
  filter(!is.na(sst)) 

# get bathymetry with bounding box
bathy <- get.bathy(lon = lonrange, lat = latrange, visualize = F, res = 15) 

# get CRS for reprojecting COBE
bathy.crs <- bathy %>% 
  as("SpatialPolygonsDataFrame") %>% 
  st_as_sf() %>% 
  st_crs() 

# make bathymetry mask
bathy.mask <- bathy %>% 
  as("SpatialPolygonsDataFrame") %>% 
  st_as_sf() %>% # retains CRS 
  dplyr::filter(layer <= 300) # get rid of points over 300 m deep

# also tried making a polygon union 
bathy.mask.union <- st_union(bathy.mask)

# make COBE sf object
cobe.sf <- cobe.tidy %>% 
  mutate(x=x-180) %>% # change lon to -180/+180 from +360
  st_as_sf(coords=c("x","y"), crs = bathy.crs) 

我尝试了 cobe.sf %>% st_intersect(bathy.mask)cobe.sf %>% st_intersect(bathy.mask.union)cobe.sf %>% st_filter(bathy.mask)cobe.sf %>% st_join(bathy.mask, left=FALSE)(它在原始代码中并且仍然适用于所有其他数据文件)。它们都生成一个空对象,这让我怀疑我遗漏了一些关于原始 COBE 数据的愚蠢信息。 COBE 是一个非常大的数据集,但我绘制了 cobe.sf 并且它具有全局范围并且看起来像法线图。有任何想法吗?以下是每个对象的摘要:

> cobe.sf
Simple feature collection with 67932249 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -179.5 ymin: -84.5 xmax: 179.5 ymax: 89.5
epsg (SRID):    4326
proj4string:    +proj=longlat +ellps=WGS84 +no_defs
# A tibble: 67,932,249 x 3
   date          sst      geometry
   <chr>       <dbl>   <POINT [°]>
 1 X1891.01.01 -1.80 (-179.5 89.5)
 2 X1891.02.01 -1.80 (-179.5 89.5)
 3 X1891.03.01 -1.80 (-179.5 89.5)
 4 X1891.04.01 -1.80 (-179.5 89.5)
 5 X1891.05.01 -1.80 (-179.5 89.5)
 6 X1891.06.01 -1.80 (-179.5 89.5)
 7 X1891.07.01 -1.80 (-179.5 89.5)
 8 X1891.08.01 -1.80 (-179.5 89.5)
 9 X1891.09.01 -1.80 (-179.5 89.5)
10 X1891.10.01 -1.80 (-179.5 89.5)
# ... with 67,932,239 more rows
> bathy.mask
Simple feature collection with 497 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -78 ymin: 35 xmax: -66 ymax: 45
epsg (SRID):    4326
proj4string:    +proj=longlat +ellps=WGS84 +no_defs
First 10 features:
   layer                       geometry
1      5 POLYGON ((-67.25 45, -67 45...
2     44 POLYGON ((-67 45, -66.75 45...
3    127 POLYGON ((-66.75 45, -66.5 ...
4    122 POLYGON ((-66.5 45, -66.25 ...
5     95 POLYGON ((-66.25 45, -66 45...
6      5 POLYGON ((-67.75 44.75, -67...
7     47 POLYGON ((-67.5 44.75, -67....
8     77 POLYGON ((-67.25 44.75, -67...
9     50 POLYGON ((-67 44.75, -66.75...
10   159 POLYGON ((-66.75 44.75, -66...
> bathy.mask.union
Geometry set for 1 feature 
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -78 ymin: 35 xmax: -66 ymax: 45
epsg (SRID):    4326
proj4string:    +proj=longlat +ellps=WGS84 +no_defs
MULTIPOLYGON (((-75.5 39.5, -75.25 39.5, -75.25...

您从 0-360 度经度到 -180-180 度经度的转换已关闭。例如,请参见此处:https://gis.stackexchange.com/questions/201789/verifying-formula-that-will-convert-longitude-0-360-to-180-to-180

数据: 请注意,我采样了 10,000 行。

# cobe.tidy
set.seed(123)
cobe.tidy <- raster::as.data.frame(cobe.raw, xy=TRUE)  %>% 
  pivot_longer(cols=c(-x, -y), names_to="date", values_to="sst") %>%
  filter(!is.na(sst))

cobe.sample <- cobe.tidy %>% 
  sample_n(10000)

原始版本: cobe 投影不正确,您的测深结果在 China/Mongolia 中(或者更确切地说,China/Mongolia 最终在 New英格兰)。

# make COBE sf object
cobe.sf <- cobe.sample %>% 
  mutate(x=x-180) %>% # change lon to -180/+180 from +360
  st_as_sf(coords=c("x","y"), crs = bathy.crs) 

ggplot() +
  geom_sf(data = cobe.sf) +
  geom_sf(data = bathy.mask.union, fill = "red")

新版本: 新英格兰附近的水深测量正确显示。

# make COBE sf object
cobe.sf <- cobe.sample %>% 
  mutate(x= ifelse(x >= 180, x-360, x)) %>% # change lon to -180/+180 from +360
  st_as_sf(coords=c("x","y"), crs = bathy.crs) 

ggplot() +
  geom_sf(data = cobe.sf) +
  geom_sf(data = bathy.mask.union, fill = "red")

要剪辑,我们使用 st_intersection()。在这种情况下,您还可以执行 cobe.sf[bathy.mask.union,].

output <- cobe.sf %>% 
  st_intersection(bathy.mask)

ggplot() +
  geom_sf(data = output) +
  geom_sf(data = bathy.mask.union, fill = "red", alpha = 0.5)

output

Simple feature collection with 10 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -77.5 ymin: 36.5 xmax: -67.5 ymax: 43.5
CRS:            +proj=longlat +ellps=WGS84
# A tibble: 10 x 3
   date          sst     geometry
   <chr>       <dbl>  <POINT [°]>
 1 X1992.11.01  8.42 (-77.5 43.5)
 2 X1896.10.01 15.4  (-67.5 40.5)
 3 X1913.09.01 19.1  (-70.5 40.5)
 4 X1911.01.01  8.03 (-67.5 40.5)
 5 X1932.12.01  7.34 (-69.5 42.5)
 6 X1928.03.01  2.67 (-70.5 42.5)
 7 X1974.07.01 20.4  (-73.5 40.5)
 8 X2013.02.01  4.32 (-67.5 43.5)
 9 X1899.11.01 11.4  (-71.5 41.5)
10 X1941.06.01 22.1  (-75.5 36.5)

看来你需要做的是

cobe.raw <- raster::stack("sst.mon.mean.nc")
cobe.rawr <- rotate(cobe.raw)

然后从那里拿走