将 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)
然后从那里拿走
我经常下载大尺度的气候数据，并将其裁剪到某个空间对象的范围内。我的代码仍然可以通过将海面温度数据裁剪为 "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)
然后从那里拿走